Valence Cristiani

Maxime Lecomte

Philippe Maurine

Revisiting Mutual Information Analysis: Multidimensionality, Neural Estimation and Optimality Proofs

Keywords: Side-channel analysis, Mutual information, Deep learning, Multidimensionality, MINE

Recent works showed how Mutual Information Neural Estimation (MINE) could be applied to side-channel analysis in order to evaluate the amount of leakage of an electronic device. One of the main advantages of MINE over classical estimation techniques is to enable the computation between high dimensional traces and a secret, which is relevant for leakage assessment. However, optimally exploiting this information in an attack context in order to retrieve a secret remains a non-trivial task especially when a profiling phase of the target is not allowed. Within this context, the purpose of this paper is to address this problem based on a simple idea: there are multiple leakage sources in side-channel traces and optimal attacks should necessarily exploit most/all of them. To this aim, a new mathematical framework, designed to bridge classical Mutual Information Analysis (MIA) and the multidimensional aspect of neural-based estimators, is proposed. One of the goals is to provide rigorous proofs consolidating the mathematical basis behind MIA, thus alleviating inconsistencies found in the state of the art. This framework allows to derive a new attack called Neural Estimated Mutual Information Analysis (NEMIA). To the best of our knowledge, it is the first unsupervised attack able to benefit from both the power of deep learning techniques and the valuable theoretical properties of MI. Simulations and experiments show that NEMIA outperforms classical side-channel attacks, especially in low-information contexts.

Introduction

Context

Side-Channel Analysis (SCA) could be defined as the process of gaining information on a secret hold by a system through leakage that comes from its practical implementation. In the most famous examples, an adversary exploits physical leakages of an electronic device such as its power consumption [START_REF] Kocher | Differential power analysis[END_REF] or Electromagnetic (EM) emanations [START_REF] Quisquater | Electromagnetic analysis: Measures and counter-measures for smart cards[END_REF] to recover a cryptographic key. Many other side-channels have been pointed out in the literature such as timing attacks [START_REF] Billy | Remote timing attacks are still practical[END_REF], cache monitoring [START_REF] Percival | Cache missing for fun and profit[END_REF] or even network packets length analysis [SSH + 14]. In any case, the problem can be reduced to the following form: an adversary is able to learn realizations of a leakage variable L, often called a trace, and aims at using it to infer information about another related secret variable S.

From an information theory point of view, the maximum amount of information one could extract from a side-channel trace is bounded by the Mutual Information I(S, L). This quantity is, indeed, central in the side-channel domain. The goals of the different actors could be summarized as follows:

• Designers aim at implementing countermeasures to decrease as much as possible I(S, L), under computational and efficiency constraints.

• Evaluators aim at estimating I(S, L) as closely as possible to assess leakages in a worst-case scenario.

• Attackers aim at developing strategies to partially or fully exploit I(S, L) in order to recover a secret.

The main problem with this paradigm is that I(S, L) is famously hard to estimate from drawn samples when the variables live in a high dimensional space, which is generally the case of L (i.e. power traces often consist of thousands of time samples). Classical MI estimators suffer from the so called "curse of dimensionality" and require an exponential (w.r.t. the dimension) amount of data to produce reliable results. This explains why, despite its valuable theoretical properties, I(S, L) is not directly used for side-channel analysis. Instead, one often compute max i I(S, L[i]) where L[i] stands for the i-st sample of the trace, but this may not represent the true available information when multiple samples leak or when there exist some dependencies between these samples. However, in a recent work [START_REF] Cristiani | Leakage Assessment through Neural Estimation of the Mutual Information[END_REF], authors took advantage of a new deep learning technique called Mutual Information Neural Estimation (MINE) [BBR + 18] to develop a side-channel tool able to reliably estimate the MI between the secret and full traces, drastically reducing the impact of high dimensionality on the estimation reliability. This tool allows one to get an absolute leakage quantification from raw traces which is helpful for designers or evaluators to perform leakage assessment. However, knowing the amount of potentially usable information is not the same as actually exploiting it to retrieve a secret, and authors left open questions regarding this tool from the attacker's point of view. Is an adversary also able to use the inherent multidimensional properties of MINE to exploit at the same time all the potential leakage sources ? And if so, what is the optimal way to do it ? This paper aims at answering these questions.

Side-channel attacks are mainly divided into two categories: supervised SCA, where the adversary can first perform a characterization of the target, and unsupervised SCA in which this profiling step is not possible. For profiled SCA, one is theoretically able to exploit all the information I(S, L) by perfectly learning the target's leakage model during the characterization phase. Deep learning attacks have been shown to effectively extract all the available information when using the negative log likelihood as loss function [START_REF] Masure | A comprehensive study of deep learning for side-channel analysis[END_REF]. Therefore the problem is closed, at least in theory, for profiled SCA. However, this is not the case for unsupervised attacks, where the true leakage model of the target is unknown to the adversary. In this situation, only a fraction of I(S, L), which value depends on the correctness of one's a priori on the leakage model, can be exploited. For example, the Correlation Power Analysis (CPA) [START_REF] Brier | Correlation power analysis with a leakage model[END_REF] is efficient for linear dependencies between the leakage and a certain function of the intermediate variable (often being the Hamming weight function). The Linear Regression Analysis (LRA) [?] also assumes a linear dependency but can handle different weights for each bit of the intermediate variable.

Mutual Information Analysis (MIA), however, has been introduced as a generic strategy able to capture any kind of dependencies. Papers addressing the theoretical background behind MIA [GBTP08, PR09, VCS09, BGP + 11] all present MIA as SCA distinguisher able to recover the correct key without any knowledge on the target nor on its leakage model. However, this leakage model free strategy only works to target non-bijective intermediate variables which makes it well suited for the DES (as the DES S-boxes are not bijectives) but not for more recent algorithms such as the AES. This explains why MIA has not often been used in practice.

A second version of the MIA allowing to target any intermediate variables (and is therefore applicable in many more contexts) has also been developed. These two versions are not separated in the literature but we decided to do so in this paper to clarify the relationship between MIA and leakage model a priori. Indeed, this second version is not leakage model free i.e. it requires some assumptions on the leakage model to work. However, one of the main advantages of this attack over CPA or LRA is that it is not limited to linear leakage model. This gain in genericity comes at the cost of efficiency: CPA has always been proved to work better than MIA in classical attack scenarios since leakage models are often linear. Therefore MIA is more seen as a great tool in theory that does not offer much in practice.

However, one of the main advantages of MIA is that it generalizes well to higher dimension variables and offers a way to potentially use a bigger part of the information contained in a side-channel trace. This has not really been used in the literature (except to extend MIA for masked implementation [PR09,BGP + 11]) due to MI estimators limitations. But recent breakthrough regarding neural estimation encourages to revisit classical MIA in order to make it highly multidimensional, to get closer to an optimal attack regarding the amount of information being used from the traces.

Even if neural estimation techniques can be applied in the leakage model free version of the MIA, we are more interested in the second version of MIA since it does not impose restrictions on the targeted algorithm. However, we argue that the mathematical framework behind this version (developed in [GBTP08, PR09, VCS09, BGP + 11]) is not complete or even wrong at some points and rely too much on intuition instead of proofs. As a result, it is difficult to derive the best way to use the new MI estimators, especially in the context of high dimensional variables, where intuition quickly falls short. That is why rebuilding a mathematical framework along with rigorous proofs on how to conduct an optimal multidimensional MIA is one of the contributions of this paper.

Contributions

1. Clarifying the State Of The Art (SOTA) around the MIA.

We explicitly split MIA into two different versions (which is not done in the SOTA), which helps understanding the need or not of an a priori on the leakage model (2.2). We then highlight inconsistencies with the second version mainly related to the fact that MIA relies on a distinguisher computing a score for each key hypotheses, but the wrong hypotheses scores are not taken into account in the analysis (2.3). This leads us to define a new generic version of MIA which objective is related to a maximization problem that includes the impact of the wrong hypotheses scores (2.4).

2. Providing rigorous proofs to analytically solve the mathematical problems emerging from our new version of MIA.

One of the main contributions of this paper, given by theorem 1 (2.5), is to solve the optimization problem defined in (2.4). Then, theorem 2 provides an extension of the analysis in the context of masking (3). Both theorems are designed to take into account the potential multidimensionality of the leakage and therefore are suited to support the use of the new neural MI estimators.

Presenting a new unsupervised multidimensional attack: the Neural Estimated Mutual Information Analysis (NEMIA).

Mathematical results are then combined with recent breakthroughs regarding neural MI estimation in high dimension. This allows to derive, to the best of our knowledge, the first unsupervised side-channel attack able to benefit from both deep learning techniques (highly multidimensional, no pre-processing of the data...) and the valuable theoretical properties of MI (4).

Providing

Simulations and experiments to support the analysis.

Simulations are provided both to empirically validate the mathematical analysis as well as to gain intuition about their meaning and about which situations are best suited for the use of NEMIA (5). Eventually, practical experiments on the ASCAD database (both on raw traces and on artificially noised traces) are conducted and shows that this new attack outperforms classical SCA strategies in terms of number of traces needed and noise resiliency (6).

Mutual Information Analysis

Background

Notations. Random variables are represented as upper-case letters such as X. They take their values in the corresponding set X depicted with a calligraphic letter. Lower case letters such as x stand for elements of X . Probability density function associated to the variable X is denoted by p X Information theory. The entropy H(X) [START_REF] Shannon | A mathematical theory of communication[END_REF] of a random variable is a fundamental quantity in information theory which indicates how much information one would gain, in average, by learning a particular realization x of X. It is defined as the expectation of the self-information log 2 (1/p X ). In a discrete context:

H(X) = x∈X p X (x) • log 2 ( 1 p X (x) ) (1)
In a side-channel environment where L represents the acquired data, one is not interested in the absolute information provided by X but rather in the amount of information revealed about a second variable such as a secret S. This is exactly what is measured by the mutual information I(S, L). It is defined as:

I(S, L) = H(S) -H(S | L) = H(L) -H(L | S) (2)
where H(A | B) stands for the conditional entropy of A knowing B:

H(A | B) = b∈B p B (b) • H(A | B = b) (3)
Another useful way to characterise I(S, L) is to express it as the Kullback-Leibler (KL) divergence between the joint distribution and the product of the marginals:

I(S, L) = D KL (p S,L || p S ⊗ p L ) = s∈S l∈L P (s, l) • log P (s, l) P (s) • P (l) (4)
Unsupervised attacks. Suppose an adversary wants to recover the secret key used by the physical implementation of a cryptographic algorithm. He has access to a set of measurements (traces) (L i ) 1≤i≤n labeled with the plaintext P i used for the encryption. The general idea of an unsupervised side-channel attack is to make a series of hypotheses k i , on a sub-part of the key and to use a distinguisher D(k) allowing to rank the different candidates. Distinguishers use statistical dependencies between traces and an intermediate variable Z k * = g(P, k * ) that depends on the plaintext and the correct key k * through a deterministic function g : P × K → Z related to the underlying algorithm. For simplicity, g(P, k) is denoted g k (P ) in the rest of the paper.

Common distinguishers such as Pearson's coefficient or coefficient of determination in a linear regression exploit some a priori on the leakage model. A common intuition about mutual information used as a distinguisher [START_REF] Gierlichs | Mutual information analysis[END_REF] is that it has been introduced precisely to reduce the need to have an a priori. It is often found in the literature (e.g. [BGP + 11]) that it aims at generality, leading to successful attacks without requiring specific knowledge or assumptions on the target. While this is true in some sense, this assertion is mitigated hereafter.

State of the art

This section presents the state of the art of MIA [GBTP08, PR09, VCS09, BGP + 11] and is organized to discuss and clarify the importance of the adversary's leakage model a priori. MIA uses a distinguisher D which takes the following form:

D(k) = I f (Z k ), L 1 (5)
with f being a function transforming the guessed intermediate variables Z k . This function is one of the main concerns of this paper. It is often called the "model" of the adversary.

The requirement of a model may seem contradictory with the claims of genericity of the MIA. Actually, this model can be replaced by the identity function making the MIA independent of any a priori on the leakage model. This version of the MIA is presented hereafter.

MIA version 1. (Leakage model free)

In its most basic form, MIA uses I(Z k = g k (P ), L) as a distinguisher, making hypotheses on k. With ϕ : Z → R n representing the leakage model of the target, L can be written as L = ϕ(Z k * ) + N , with N being a random variable independent of Z k for all k, and representing the noise. With these notations, the distinguisher becomes:

D(k) = I Z k , ϕ(Z k * ) + N (6)
Proposition 1. This distinguisher is maximized for the correct key hypothesis k * .

Proof. Using equation 2:

D(k * ) -D(k) = H(L) -H(L | Z k * ) -H(L) -H(L | Z k ) = H ϕ(Z k * ) + N | Z k -H ϕ(Z k * ) + N | Z k * (7)
Since adding knowledge can only decrease entropy:

D(k * ) -D(k) ≥ H ϕ(Z k * ) + N | Z k , Z k * -H ϕ(Z k * ) + N | Z k * (8)
Now using the independence of N and the fact that ϕ(Z k * ) is entirely determined by Z k * :

D(k * ) -D(k) ≥ H(N ) -H(N ) ≥ 0 (9)
which concludes the proof.

This strategy does not require any assumption on the leakage model of the target. However, it only works if the correct key hypothesis is distinguishable from the other ones, or, in other words, if D(k) < D(k * ), ∀k = k * , which is not guaranteed by proposition 1. An important property of the MI is that it is preserved by injective transformations of the input variables. So if different key hypotheses yield Z k variables differing from each other only by a permutation (for example if the g k functions are bijective), I(Z k , L) would be constant for all k and the distinguisher could not discriminate key candidates. Therefore, g k has to be non-injective. For example, one could target the output of the first round DES S-box.

While this form of MIA is effectively leakage model free, it comes with a huge constraint since in many interesting cases g k is bijective. In the AES case, this means that one cannot target the output of the first S-box since Sbox[k * ⊕ P ] is bijective with P . In [START_REF] Prouff | Theoretical and practical aspects of mutual information based side channel analysis[END_REF], authors suggest to target the bitwise addition between two S-box outputs during the MixColumns operation. This requires making hypotheses on 16 bits of the key (leading to 2 16 MI computations). Moreover, it is only feasible if this operation leaks enough information which may not be the case in practice. Indeed, for hardware implementations, this step is usually fully combinatorial and does not use any register. This explains why most of the MIA experiments in the literature have been performed on the DES.

MIA version 2. (Leakage model dependant)

It is still possible to target Z k * = g k * (P ) for bijective g k functions. The idea is to apply a non-injective function f to Z k and use I(f (Z k ), L) as distinguisher. The application of f create a partition of Z so f will be called the "partition function" in the rest of this paper. Since no data transformation can create information (this is the so called data processing inequality [START_REF] Beaudry | An intuitive proof of the data processing inequality[END_REF]), the application of f inevitably decreases the initial information: ∀f, ∀k, I(f (Z k ), L) ≤ I(Z k , L). The goal is then to find a function that decreases more I(Z k , L) than I(Z k * , L) and therefore, enhance the discriminating power of the analysis.

For example, assuming that bits leak independently, [START_REF] Gierlichs | Mutual information analysis[END_REF] proposes to drop one bit of Z. This is equivalent to redefine the intermediate variable as a restrictive number of bits of g k * (P ), and apply MIA version 1 with no partition function. Another idea is to use a guessed version φ of the leakage model ϕ. Indeed, I ϕ(Z k ), ϕ(Z k * ) + N is clearly maximized for k = k * . Therefore, if φ is not too far from ϕ, I φ(Z k ), ϕ(Z k * ) + N may still be maximized for k = k * . It is shown in [START_REF] Veyrat-Charvillon | Mutual information analysis: How, when and why?[END_REF] that error in the approximation of ϕ may be less penalizing than for other attacks.

In addition, MIA is more flexible in the sense that it is not limited to exploit linear dependencies and gives an option to mount a successful attack with any leakage model. However, it should be emphasized that, for this version, the adversary must have a good enough a priori on the leakage, otherwise, the attack is unsuccessful. A suitable choice for the partition function necessarily uses assumptions on ϕ.

While we think this point needed to be clarified, we do not see this as a criticism of MIA. Actually, proposition 2 hereafter proves that hopes of finding a leakage model free strategy able to target a bijective intermediate variable Z k * are vain, even outside the context of MIA. Proposition 2. Let g k be a bijective map for all k. For any strategy S which takes as input a set of traces L = ϕ(g k * (P i )) 1≤i≤n and outputs a ranking of the different key hypotheses, there exists a leakage model φ that would rank k * in the last position such that the attack completely fails.

Proof. First, apply S on traces obtained through any leakage model ϕ 0 and denote by k the last key returned by S. Now, for all P , define φ0 (g k * (P )) = ϕ 0 (g k(P )), which completely defines φ0 as gk is bijective. Applying S on traces obtained through φ0 would now rank k * in the last position.

This proposition shows that there does not exist any generic distinguisher, that would both:

1) Exploit statistical dependencies between traces and an intermediate variable bijectively related to the plaintext.

2) Work whatever the leakage model of the target.

Since MIA version 2, with a fixed partition function, verifies 1), it necessarily fails for some leakage models or, in other words, has to use an assumption on the leakage model to succeed. Even though it requires an analysis on what partition function should be used, the rest of this paper is more focused on MIA version 2 since it is more generic in the sense that it can be applied in many more situations.

About the distinguishability

Both versions of the MIA guarantee that D(k) is maximized for k = k * in theory. However, we argue that this is not enough to guarantee a successful attack in practice, when noise comes into play. What is really important is the difference between D(k * ) and the others, or in other words, the distinguishability of the correct hypothesis through the distinguisher D. The idea found in the literature is that for a wrong key hypothesis: «false predictions will form a partition corresponding to a random sampling of [L] and therefore simply give scaled images of the global side-channel probability density function. Hence, the estimated mutual information will be equal (or close) to zero in this case.» [BGP + 11].

We do not agree with this fact since the wrong hypotheses scores totally depend on the partition function f and on g k . As explained in the previous section, if the g k 's are bijective, all the scores would be equal if f is also bijective. This fact is well noted in all the papers about MIA but we would like to emphasize that even for non-bijective f the wrong hypotheses score depends on the "degree of bijectiveness" of f . Intuitively, the more compact f is (in the sense of more collisions through f ) the smaller the wrong hypotheses scores would be. But the same is true for the correct score which means that there is a trade-off between how much one wants to decrease I f (Z k ), L for the wrong k and keep I f (Z k * ), L high, to enhance the distinguishability.

Towards an optimal partition function f

In the SOTA, the partition function is not seen as a parameter on which a maximization research could be done. Therefore, no research on finding the optimal function f has been conducted. It is generally fixed to one or two constant choices, except in [START_REF] Prouff | Theoretical and practical aspects of mutual information based side channel analysis[END_REF] where authors proposed that f could be a generic function. However, it is stated that the adversary:

«does not need a good linear approximation of ϕ but only a function [f ] such that the mutual information

I f (Z k * ), ϕ(Z k * ) is non-negligible » [PR09].
Again, this condition is necessary but not sufficient. Even if bijective functions are excluded one can create the following f 0 function such that:

f 0 (x) = 0, if x ∈ {0, 1} x, else (10)
Being almost the identity function, f 0 is such that I f 0 (Z k * ), ϕ(Z k * ) is high but would have a very low discriminating power. This shows that the wrong hypotheses scores can not be left out of the analysis. One typically wants to find the f function maximizing the distinguishability of the correct hypothesis. Therefore, let us define the optimal set of functions F opt as:

F opt = arg max f : Z→R n I f (Z k * ), L -max k =k * I f (Z k ), L 2 (11)
F opt is a set since the maximum is reached by an infinite amount of functions. Indeed, if f opt ∈ F opt , for any bijection b, b • f opt is also in F opt since bijections do not affect MI.

Analytical resolution

Being consistent with proposition 2, F opt depends on L and therefore on the leakage model. Since knowledge on ϕ is required anyway, this section assumes a full knowledge on ϕ in order to conduct an analytical analysis to find the optimal f function. Traces are also supposed to be acquired in an ideal set-up, without noise, so that, at least for a significant sub-part of the trace, L = ϕ(Z k * ). Bounds taking into account imperfect knowledge on ϕ as well as noise will be given in section 2.7.

A natural choice for the partition function would be to take f = ϕ because it maximizes the left term in (11):

I f (Z k * ), ϕ(Z k * ) .
But it is a priori possible to find a function that would maximize the global objective without maximizing the left term of (11). Th. 1 actually proves that it is not possible and that whatever the leakage model, ϕ ∈ F opt . The main demonstration requires the use of a helper which is introduced in the form of a lemma hereafter.

Lemma 1. Let f : Z → R n be any function. For any leakage model ϕ: Z → R n there exists a decomposition of f into f = f 2 • f 1 , with f 1 : Z → N, f 2 : N → R n , satisfying the two following properties: 1) ∃ f 3 : Im f 1 → R n such that f 3 • f 1 = ϕ 2) ∀z ∈ Z, f 2 | f1 ϕ -1 ({ϕ(z)}) is bijective of reciprocal f -1 2 | f2•f1 ϕ -1 ({ϕ(z)})
Proof. The proof is given in appendix A.

Theorem 1. Let P follow a uniform distribution. Let Z k represent the hypothetical intermediate variables such that Z k = g k (P ) with bijective g k 's. Let ϕ: Z → R n be the leakage model of the target so that L = ϕ(Z k * ). Then, ϕ ∈ F opt . Proof. Let S f = I f (Z k * ), L -max k =k * I f (Z k ),
L represent the distinguishability score for a given function f such that:

F opt = arg max f : Z→R n {S f }
Since all the Z k follow a uniform distribution (P follows a uniform distribution and the g k functions are bijective), the entropy H(f (Z k )) is equal for all k. Then using

I(A, B) = H(A) -H(A | B): S f = -H f (Z k * ) | L + min k =k * H f (Z k ) | L (12)
Symmetrically, using

I(A, B) = H(B) -H(B | A): S f = -H L | f (Z k * ) + min k =k * H L | f (Z k ) (13)
Let f : Z → R n be any function. Applying lemma 1, there exist f 1 and f 2 satisfying the two properties given in lemma 1, such that

f = f 2 • f 1 . The goal is to show that S f ≤ S ϕ .
The proof is divided into two phases: first show that S f ≤ S f1 using (12), then show that S f1 ≤ S ϕ using (13). Let us start with (12):

S f = -H f 2 • f 1 (Z k * ) | L + min k =k * H f 2 • f 1 (Z k ) | L ≤ -H f 2 • f 1 (Z k * ) | L + min k =k * H f 1 (Z k ) | L ( 14
)
since applying f 2 in the second term can only decrease entropy (see lemma 2). The goal is now to remove f 2 in the first term:

-H f 2 • f 1 (Z k * ) | L = l∈L f2∈Im f2 P (l) • P ( f2 | l) • log(P ( f2 | l)) (15) P ( f2 | l) = P f 2 • f 1 (Z k * ) = f2 | ϕ(Z k * ) = l = P f 1 (Z k * ) ∈ f -1 2 ( f2 ) | ϕ(Z k * ) = l (16) Knowing that ϕ(Z k * ) = l implies that Z k * ∈ ϕ -1 ({l}) and also that f 1 (Z k * ) ∈ f 1 (ϕ -1 ({l})).
Let A l denotes f 1 (ϕ -1 ({l})) to avoid heavy notations. Then:

ϕ(Z k * ) = l =⇒ f 1 (Z k * ) ∈ A l =⇒ f 1 (Z k * ) ∈ f -1 2 (f 2 (A l )) (17)
which means that:

P ( f2 | l) = P f 1 (Z k * ) ∈ f -1 2 | f2(A l ) ( f2 ) | l if f2 ∈ f 2 (A l ) 0 else (18) Lemma 1 states that f 2 | A l is bijective of reciprocal f -1 2 | f2(A l ) , so if f2 ∈ f 2 (A l ): P ( f2 | l) = P f 1 (Z k * ) = f -1 2 | f2(A l ) ( f2 ) | l (19)
Let us plug this result back into (15):

-H f 2 • f 1 (Z k * ) | L = l∈L f2∈f2(Al) P (l) • P f 1 (Z k * ) = f -1 2 | f2(A l ) ( f2 ) | l • log P f 1 (Z k * ) = f -1 2 | f2(A l ) ( f2 ) | l (20)
Now, one can apply the following change of variable in the second sum:

f1 = f -1 2 | f2(A l ) ( f2 ): -H f 2 • f 1 (Z k * ) | L = l∈L f1∈Al P (l) • P f 1 (Z k * ) = f1 ) | l • log P f 1 (Z k * ) = f1 ) | l (21) Finally, since P f 1 (Z k * ) = f1 ) | l = 0 when f1 ∈ Im f 1 \ A l ,
one can artificially add some terms equal to 0 in the second sum:

-H f 2 • f 1 (Z k * ) | L = l∈L f1∈Im f1 P (l) • P ( f1 | l) • log P ( f1 | l) = -H f 1 (Z k * ) | L (22)
Applying this result to (14) gives:

S f ≤ -H f 1 (Z k * ) | L + min k =k * H f 1 (Z k ) | L S f ≤ S f1 (23)
which concludes the first step of the demonstration.

iacrtans class documentation

Now the goal is to show that S f1 ≤ S ϕ . Lemma 1 guarantees that there exists f 3 such that f 3 • f 1 = ϕ. Let us use this in (13):

S f1 = -H L | f 1 (Z k * ) + min k =k * H L | f 1 (Z k ) ≤ -H L | f 1 (Z k * ) + min k =k * H L | f 3 • f 1 ϕ (Z k ) (24)
since applying f 3 to the known variable can only increase the global entropy (see lemma 3). Now using L = ϕ(Z k * ):

-H L | f 1 (Z k * ) ≤ 0 -H L | f 1 (Z k * ) ≤ -H ϕ(Z k * )|ϕ(Z k * ) = 0 (25)
Therefore:

S f1 ≤ -H L | ϕ(Z k * ) + min k =k * H L | ϕ(Z k ) S f1 ≤ S ϕ (26)
Finally, using both part of the demonstration:

S f ≤ S f1 ≤ S ϕ (27)
which ensures that ϕ is better than any other functions and so that ϕ ∈ F opt .

Remark 1. Demonstration of Th. 1 would have worked exactly the same if one had first fixed a particular hypothesis k, and tried to maximize

S f,k = I f (Z k * ), L -I f (Z k ), L .
Therefore, for each k, ϕ maximizes the distance between the score of k * and k which is an even stronger version of the theorem. One could not be sure that such a function would exist a priori, that is why F opt has not been defined with this criterion. However, this shows a posteriori that Th 1 is still valid even if one decides to redefine F opt , for example using the distance with the mean (instead of the maximum) of the wrong hypotheses scores.

Interpretation. This theorem tells that to conduct an optimal MIA, one has to transform the targeted variable Z k by applying the leakage model ϕ (or any bijection of ϕ) and use I(ϕ(Z k ), L) as a distinguisher. Note the multidimensional aspect of this theorem since both ϕ(Z k ) and L can live in high dimensional space. This is a key point in this paper that will be discussed in detail in section 4.2 which bridges this theorem with newest multidimensional MI estimators in order to derive a new attack. Note that this theorem also implies that if the leakage model is itself bijective, MIA is not a valid strategy since the distinguishability score would be bounded by 0.

Selecting leakage model a priori

In a real-life experiment, one might not perfectly know the leakage model ϕ but only an estimation φ. This section provides a procedure to evaluate the correctness of φ, helping to choose from multiple guesses φ1 , . . . , φn . This test relies on the following observation:

Proposition 3. Let L = ϕ(Z k * ) + N , with N an independent random variable representing the noise. Then: ϕ = max f [I f (Z k * ), L ]
Proof. On one hand:

I f (Z k * ), L = H(L) -H L | f (Z k * ) ≤ H(L) -H ϕ(Z k * ) + N | f (Z k * ), ϕ(Z k * ) ≤ H(L) -H(N ) (28)
and on the other hand:

I ϕ(Z k * ), L = H(L) -H L | ϕ(Z k * ) = H(L) -H ϕ(Z k * ) + N | ϕ(Z k * ) = H(L) -H(N ) (29) Then: I ϕ(Z k * ), L ≥ I f (Z k * ), L (30)
which concludes the proof.

The identity function obviously also maximizes: I f (Z k * ), L so combining this with proposition 3:

I Z k * , L = I ϕ(Z k * ), L (31)
or,

I Z k * , L = max k [I ϕ(Z k ), L ] (32)
Then, if k * is known (for example in an evaluation setup) one can use equation 31 and estimate I Z k * , L and I φ(Z k * ), L and compare them. If the latest is too far from I Z k * , L , one may reject φ as being a good approximation of the true underlying leakage model. If k * is unknown, the adversary can still use equation 32 estimating I φ(Z k ), L for all k, and comparing the maximum with I Z k0 , L (k 0 can be chosen randomly since all the Z k variables are just permutation of each other which does not affect MI). Note that this test is only a rejection test since passing the test does not guarantee a good estimation of ϕ: for example, the identity function always passes the test.

Leakage model uncertainty and noise

Let assume that the adversary has chosen a given estimation φ of ϕ. Let also assume that the ideal data L = ϕ(Z k * ), used in theorem 1, are now noisy so that the acquired data takes the following form: L = ϕ(Z k * ) + N , with N an independent random variable. This section aims at complementing theorem 1 by lower bounding the distinguishably score S φ that one would get in practice in such a context:

S φ = I φ(Z k * ), L -max k =k * I φ(Z k ), L (33)
Our goal is to compare S φ with the optimal score S ϕ (from theorem 1) that one would get with the perfect knowledge of ϕ and un-noised data such that:

S ϕ = I ϕ(Z k * ), L -max k =k * I ϕ(Z k ), L (34)
Proposition 4. S φ is lower-bounded by the following inequality:

S φ ≥ S ϕ -H(N ) -H ϕ(Z k * ) | φ(Z k * ) -max k =k * H φ(Z k ) | ϕ(Z k ) ( 35
)
Proof. Using the same argument as in (13) one has:

S φ = -H ϕ(Z k * ) + N | φ(Z k * ) + min k =k * H ϕ(Z k * ) + N | φ(Z k ) (36)
Since removing noise on the right term can only decrease entropy:

S φ ≥ -H ϕ(Z k * ) + N | φ(Z k * ) + min k =k * H ϕ(Z k * ) | φ(Z k ) (37)
Now since H(A + B) ≤ H(A) + H(B) and using the independence of N :

S φ ≥ -H(N ) -H ϕ(Z k * ) | φ(Z k * ) + min k =k * H ϕ(Z k * ) | φ(Z k ) (38) Using H(A | B) ≥ H(A | C) -H(B | C)
which can be shown through information Venn diagram:

S φ ≥ -H(N ) -H ϕ(Z k * ) | φ(Z k * ) + min k =k * H ϕ(Z k * ) | ϕ(Z k ) -H φ(Z k ) | ϕ(Z k ) ≥ -H(N ) -H ϕ(Z k * ) | φ(Z k * ) + min k =k * H ϕ(Z k * ) | ϕ(Z k ) -max k =k * H φ(Z k ) | ϕ(Z k ) (39)
Now let S ϕ appear in the equation:

min k =k * H ϕ(Z k * ) | ϕ(Z k ) = min k =k * H ϕ(Z k * ) | ϕ(Z k ) - 0 H ϕ(Z k * ) | ϕ(Z k * ) = S ϕ (40) So: S φ ≥ S ϕ -H(N ) -H ϕ(Z k * ) | φ(Z k * ) -max k =k * H φ(Z k ) | ϕ(Z k ) (41)
which concludes the proof.

This proposition describes the impact of the noise and leakage model approximation in a quantitative way. Its qualitative interpretation is fairly intuitive. It clearly shows that one has two strategies to get closer to the optimal score: reducing the noise entropy or improving his guess on φ. When H(N ) tends towards 0 and φ gets closer to ϕ , S φ tends towards the optimal score S ϕ . It also captures the fact that bijective errors do not impact the outcome of the attack since if there exists a bijection between φ(Z k ) and ϕ(Z k ), both terms

H ϕ(Z k * ) | φ(Z k * ) and max k =k * [H φ(Z k ) | ϕ(Z k ) ]
would be equal to 0.

MIA against masked implementations

Masking is one of the most widely used countermeasures to protect implementations of block ciphers against side-channel analysis [START_REF] Chari | Towards sound approaches to counteract power-analysis attacks[END_REF]. The idea is to split each sensitive intermediate value Z, into d shares, following the relation:

Z = Z 1 * Z 2 * • • • * Z d ( 42
)
for a group operation * . The d -1 shares Z 2 , ..., Z d are randomly chosen and the last one, Z 1 , is processed such that (42) is satisfied. Assuming the masks are uniformly distributed, the knowledge of d -1 shares does not tell anything about Z. However, partial knowledge on the d shares can be exploited to retrieve information on Z. That is why, to defeat masking, one should use a distinguisher able to combine the leakage of at least d samples of the traces (assuming masks do not leak at the same time). Higher-order correlation attacks [START_REF] Messerges | Using second-order power analysis to attack dpa resistant software[END_REF] exploit a combining function, C : R d → R , which transforms a multidimensional leakage into a single value such that the output of C correlates with Z.

The optimal combining function is unknown but most of the time, the centered product between the shares is used [START_REF] Prouff | Statistical analysis of second order differential power analysis[END_REF].

MIA, a natural choice against masking

Although higher-order CPA attacks lead to successful key recoveries, they are not optimal from an information-theoretic point of view. Indeed, by the data processing inequality [START_REF] Beaudry | An intuitive proof of the data processing inequality[END_REF], the application of the combining function leads to an information loss. Opposed to Pearson's correlation, mutual information can deal with dependencies of multidimensional variables. Therefore, no combining function is required which makes MIA a very natural strategy against masked implementations. An extension of MIA in the context of masking has been proposed in [START_REF] Prouff | Theoretical and practical aspects of mutual information based side channel analysis[END_REF] and [BGP + 11]. The idea is very similar to the non-masked case. Concepts of MIA versions 1 and 2 still apply and one can use I(f (Z k ), L) as a distinguisher.

About the partition function in the presence of masking

Using I(f (Z k ), L) as distinguisher still raises the question of the optimal f function. Th. 1 cannot be applied straightforwardly since, for masked implementation, the leakage cannot be expressed as a deterministic function ϕ(Z k * ) modulo some noise. Instead, with Z i representing the shares, one now has:

L = i ϕ i (Z i ) (43)
for some functions ϕ i : Z → R n . Note that, as for the unmasked case, a noise-free version of the leakage is first considered to simplify the analysis. Noise will be added in section 3.3. Most of the time, the ϕ i supports can be supposed disjoint (i.e. leakages of the shares do not overlap). In that case, the leakage vector could be summarized as:

L = [ϕ 1 (Z 1 ), . . . , ϕ d (Z d )] (44)
with ϕ i taking its values in a subspace of R n . Even with this simplification, we could not solve analytically the problem of finding an optimal partition function, or, in other words, a function f ∈ F opt as defined in (11). However, we still give some useful insights in the common case of boolean masking on a device leaking the Hamming weight (or Hamming distance with a known value) of the shares.

For this specific case, [BGP + 11] tried to use the Hamming weight as well as the identity function for f (they were attacking the output of a DES S-box, therefore a non-injective intermediate variable). The Hamming weight produced better results. Their justification is that the Hamming weight is closer to the underlying leakage model of the circuit. We do not find this justification intuitive since even in the ideal case where the leakage could be expressed as:

L = [HW (Z k * ⊕ M ), HW (M )] (45)
HW (Z k * ) is not directly related to any physical leakage. More generally, there is no proof that if all shares leak with the same leakage model ϕ, taking f = ϕ is the optimal (or even a good) option. However, a closer inspection reveals that, in the specific case of a Hamming weight leakage model, it is possible to give a formal justification for using f = HW . This is the purpose of Th. 2. Note that such a justification can not be found in the SOTA.

Let us introduce F Lef t as the left part of equation 11:

F Lef t = arg max f : Z→R n I f (Z k * ), L (46)
This set does not consider the wrong hypotheses. Therefore it is not hard to find a function f ∈ F Lef t : the identity or any bijective function works. The problem is that with a bijective map, I f (Z k * ), L = I f (Z k ), L for any k. However, a non-injective function f such that f ∈ F Lef t would naturally decrease I f (Z k ), L and create some distinguishability. Such a function is not a priori likely to exist. But the following theorem shows that, while being highly non-injective, HW ∈ F Lef t .

Theorem 2. Let L represent the leakage of a masked variable Z k * with a mask M . Let both shares follow any bijection b 1 and b 2 of a Hamming weight leakage model so that:

L = b 1 HW (Z k * ⊕ M ) , b 2 HW (M ) (47)
Then, HW ∈ F Lef t or in other words:

I HW (Z k * ), L = I Z k * , L .
The following proof may be generalizable to higher-order (see section 5.3 for an empirical validation), but for simplicity, only first-order masking is considered here.

Proof. Since bijective transformations do not impact mutual information, one can consider without loss of generality that:

L = HW (Z k * ⊕ M ), HW (M ) (48)
Now let us evaluate I f (Z k * ), L using eqation 4:

I f (Z k * ), L = f ∈f (Z) l∈L P ( f , l) • log P ( f , l) P ( f ) • P (l) (49)
One can split the first sum by summing on z instead of f :

I f (Z k * ), L = z∈Z l∈L P (z, l) • log P (l | f (z)) P (l) = z∈Z l∈L P (z) • P (l | z) • log P (l | f (z)) P (l) (50)
Since the identity function is bijective and maximizes this quantity, it would be enough to show that P (l | HW (z)) = P (l | z) for any given z and a given l = [HW (z ⊕ m), HW (m)] for a fixed m. Let us start by the latter term:

P (l | z) = P (HW (m)) • P (HW (z ⊕ m) | z, HW (m)) (51)
To compute the right term one can evaluate the cardinal of the set M of all the masks m satisfying the following conditions:

1) HW (m ) = HW (m) 2) HW (z ⊕ m ) = HW (z ⊕ m)
and divide by the number of byte with a Hamming Weight of HW (m) which is 8 HW (m) .

To evaluate this cardinal, we first show an invariance property. For any m ∈ M, let n m denotes the number of bits set to 1 in m such that there is also a bit set to 1 at the same position (0 to 7) in z. Then:

HW (z ⊕ m ) = HW (m ) + HW (z) -2 * n m ⇐⇒ n m = HW (m ) + HW (z) -HW (z ⊕ m ) 2 (52)
Now since m satisfies the above two conditions:

n m = HW (m) + HW (z) -HW (z ⊕ m) 2 (53)
which does not depend on m anymore. As n m has to be a positive integer, the above equation shows that:

HW (m) + HW (z) -HW (z ⊕ m) / ∈ 2N =⇒ M = ∅ (54)
This allows us to define a generic n as:

n = HW (m)+HW (z)-HW (z⊕m) 2 , if HW (m) + HW (z) -HW (z ⊕ m) ∈ 2N -1, otherwise (55)
so that ∀m ∈ M, n m = n.

Reciprocally, one can see that each byte m such that HW (m ) = HW (m) and n m = n is in M. So to form a valid m ∈ M one has to choose first the position of the n '1s' superposing with the '1s' in z, which lead to HW (z) n possibilities. Then, choose the positions of the remaining '1s', which lead to 8-HW (z) HW (m)-n possibilities. Therefore, with the convention l k = 0 when k is strictly negative:

P (HW (z ⊕ m) | z and HW (m)) = HW (z) n • 8 -HW (z) HW (m) -n • 1 8 HW (m) (56)
Injecting this into (51) gives:

P (l | z) = 8 HW (m) 2 8 • HW (z) n • 8 -HW (z) HW (m) -n • 1 8 HW (m) = 1 2 8 • HW (z) n • 8
Now using result from (56):

A = z s.t. HW (z )=HW (z) 1 8 HW (z) • HW (z ) n • 8 -HW (z ) HW (m) -n • 1 8 HW (m) = HW (z) n • 8 -HW (z) HW (m) -n • 1 8 HW (m) (60)
since all the terms are constant in the sum and there are exactly 8 HW (z) of them. Now plugging this into (58) gives:

P (l | HW (z)) = 1 2 8 • HW (z) n • 8 -HW (z) HW (m) -n = P (l | z) (61) Thus, I HW (Z k * ), L = I Z k * , L (62)
which ensures that HW ∈ F lef t and concludes the proof.

Interpretation. This theorem shows that when the shares leak in Hamming weight, it is sound to use f = HW in practice because it creates some distinguishability decreasing the information only for the wrong hypotheses. Since the Hamming distance with a computable value can be rewritten as a Hamming weight, it also works in that case. However, Th. 2 is not generalizable to any leakage model ϕ (for example on 3 bits words, ϕ = 2b 1 + b 2 + b 3 gives a counter-example). Knowing if there exists a generic strategy against masking (depending on ϕ but working for any ϕ) or if one will always be condemned to work on a case-by-case basis is an interesting question and may be handled in future works.

Remark 2. Note that since I Z k , L = I HW (Z k ), L = max k [I HW (Z k ), L ], the procedure described in section 2.6 can also be applied on a masked implementation, to test the validity of the Hamming weight leakage model hypothesis. If the Hamming weight is too far from the true model, a practical alternative is to use only specific bits of the unmasked variable as partition function. An example of this is given in section 6.

Considering the distinguishability score:

S f = I f (Z k * ), L -max k =k * I f (Z k ), L (63)
HW has not been shown to be optimal. However, a partial result can be given introducing the concept of "wider" function.

Definition 1. A function f is said wider than g if there exists another function h such that: h • f = g.

Corollary 1. Let L be defined as in (47). Then, for any function h wider than HW , S HW ≥ Sh.

Proof. The proof is given in appendix B.

Even though we do not conjecture so, a function doing better than the HW may exist. But a straightforward consequence of Th. 2, given by corollary 1, is that HW has a better distinguishability score than any other wider function.

Noise and multidimensionality

The advantage of MINE is to be able to exploit the information contained in multiple samples at the same time. In a Hamming weight leakage scenario, the Hamming weight of a variable is probably not going to leak perfectly on a single sample. Instead, multiple samples may leak a noisy version of it. To ensure that it is sound to use MINE to mount a new attack in the case of masking, one would need a multidimensional version of Th.2. This is exactly the purpose of corollary 2, in which the noise is directly included.

In the context of masking the actual useful part of the leakage could be expressed as:

L = b 1 HW (Z k * ⊕ M ) + N 1 , . . . , b m1 HW (Z k * ⊕ M ) + N m1 ,
b1 HW (M ) + N1 , . . . , bm2 HW (M ) + Nm2 (64) with b i and bj being bijective maps, and N i and Nj being discrete noise variables independent of the shares. The following corollary shows that Th 2 is still valid in that case.

Corollary 2. Let L be defined as in ( 64). Then, one still has HW ∈ F Lef t as defined in (46).

Proof. As for Th. 2, one can drop, without loss of generality, the bijections in L as they do not affect the MI. Let N be the noise vector [N 1 , . . . , N m1 , N1 , . . . , ] and L the noise-free version of the leakage so that L = L + N . As for Th. 2, it is enough to show that P (l | HW (z)) = P (l | z) for any given l and z. Decomposing on all the possible values of the noise one has:

P (l | z) = n∈N P (n) • P (L = l | z and n) = n∈N P (n) • P ( L = l -n | z) ( 65
)
Since L is noise free, it consists of the repetition of the same two variables: HW (Z k * ⊕ M ) (m 1 times) and HW (M ) (m 2 times). So for the probability P ( L = l -n | z) to be non-zero, the vector l -n should be constant on its first m 1 coordinates, and constant on its m 2 last one. Let N c be the subset of N verifying the precedent property. If n / ∈ N c , then:

P ( L = l -n | z) = P ( L = l -n | HW (z)) = 0 (66) Else, if n ∈ N c , then, with a n = (l -n)[1], b n = (l -n)[m 1 + m 2 ] and L = [HW (Z k * ⊕ M ), HW (M )]: P ( L = l -n | z) = P ( L = [a n , b n ] | z) (67)
So (65) can be rewritten as:

P (l | z) = n∈Nc P (n) • P ( L = [a n , b n ] | z) (68) Since, Th. 2 tells that P ( L = [a n , b n ] | z) = P ( L = [a n , b n ] | HW (z)): P (l | z) = n∈Nc P (n) • P ( L = [a n , b n ] | HW (z)) P (l | z) = n∈Nc P (n) • P ( L = l -n | HW (z)) P (l | z) = P (l | HW (z)) ( 69
)
which concludes the proof.

This corollary shows that it is sound to use I(HW (Z k ), L) as distinguisher even when considering a noisy multidimensional leakage vector. Th. 2 still applies and MINE may benefit from the different leakage sources resulting in an attack (presented in the next section) exploiting a bigger amount of the available information.

Neural Estimated Mutual Information Analysis (NEMIA)

This section aims at formally describing the new attack proposed in this paper. Note that throughout this work, a tool able to compute I(Z, L) with high dimensional variables has been assumed to exist. Obviously, this research has been driven by recent progress regarding neural estimation techniques. But this work is not absolutely related to MINE. It would stay sound with any MI estimator able to work in high dimension. In particular, any progress in the field, which is likely to happen since it is a very active domain, would instantly impact the attack efficiency. In this work, the most basic version of MINE is used. It should be seen as a proof of concept with almost no hyper-parameters tuning and without considering recent optimizations nor improvements in the technique (non-exhaustively: [CL20, LSN + 19, CABH + 19]). A study focused on deep learning optimizations would be interesting but is out of the scope of this paper. Basic principles of MINE are recalled hereafter.

Mutual Information Neural Estimation

Technical details about the utilization of MINE in a side-channel context can be found in [START_REF] Cristiani | Leakage Assessment through Neural Estimation of the Mutual Information[END_REF]. However, a high-level picture is still given in this section. The general idea is to express I(Z, L) as the Kullback-Leibler divergence between the joint distribution and the product of the marginals:

I(Z, L) = D KL (p Z,L || p Z ⊗ p L ).
Then, to exploit the Donkser-Varadhan variational formulation of the KL-divergence that states that if p and q are two densities defined over a compact set Ω ∈ R d :

D KL (p || q) = sup T : Ω→R [E p [T ] -log(E q [e T ])] (70)
This allows to express MI as a supremum. Then, the following loss function can be defined:

L(θ) = E p Z,L [T θ ] -log(E p Z ⊗p L ) [e T θ ]) (71)
and deep learning techniques can be applied to maximize this loss over all the functions T θ parametrized by a neural network with parameters θ ∈ Θ. The objective function should converge towards the supremum so that its final value constitutes the MI estimation. Formally:

Definition 2. (MINE) Let A = {(z 1 , l 1 ), . . . , (z n , l n )} and B = {( ∼ z 1 , ∼ l 1 ), . . . , ( ∼ z n , ∼
l n )} be two sets of n empirical samples respectively from p Z,L and p Z ⊗ p L . Let F = {T θ } θ∈Θ be the set of functions parametrized by a neural network. MINE is defined as follows:

I(S, X) n = sup T ∈F E A [T ] -log(E B [e T ]) (72)
where E X [•] stands for the expectation empirically estimated over the set X .

In practice one only has samples from the joint distribution: A = {(z 1 , l 1 ), . . . (z n , l n )} of the labeled traces. Samples from the product of the marginals can be artificially generated by shuffling the variable L using a random permutation ρ: B = {(z 1 , l ρ(1) ), . . . , (z n , l ρ(n) )}. [START_REF] Cristiani | Leakage Assessment through Neural Estimation of the Mutual Information[END_REF] is the overfitting. Indeed, the loss function may overestimate the true MI. Therefore, one can introduce a validation loss function to detect overfitting and to produce a more reliable estimation. The idea is to split A and B into training datasets A t and B t and validation datasets A v and B v . Then, only the training datasets are used for back-propagation so that the loss function evaluated on the validation datasets cannot overestimate the MI. That is why only validation loss functions are considered/plotted in this paper. For robustness, the MI estimation is not set to be the supremum of the validation loss, but instead, the supremum of a moving average along the epochs with a window size of w which depends on the variability between epochs (w = 10 in this paper).

Validation loss function. One of the main problems of MINE pointed out in

Architecture. The network's input layer consists of a concatenation of both Z and L variables. Authors in [START_REF] Cristiani | Leakage Assessment through Neural Estimation of the Mutual Information[END_REF] have shown that the representation of Z is important and that one should use the One-Hot Encoding (OHE) or a binary encoding of Z (unless otherwise specified we used the OHE in this paper). The output layer is a single neuron as the function T output has to be a real value. Other layers are not specified and should be adapted to the underlying problem (e.g. convolutional layers to counter jitter or traces misalignment).

For our experiments, unless otherwise specified, a simple Multi-Layer Perceptron (MLP) is used with two hidden layers of 20 neurons each. A batch normalization layer is added after the first layer and dropout layers with p = 0.25 are inserted after each hidden layer in order to mitigate overfitting. The activation function is set to the Exponential Linear Unit (ELU) and the batch size to 1000. The validation dataset represents 20 percent of the full dataset.

Multidimensional paradigm

MINE is by essence a tool that estimates MI in a multidimensional way, enabling to compute the MI between f (Z k ) and full (or at least large part of) traces. This was not possible with classical MI estimators which do not scale with high dimensional variables. Until now, MIA was only performed with the following distinguisher:

D old (k) = max i I(f (Z k ), L[i]) (73)
where L[i] represents the i-th sample of the trace. This way, trace dimension is kept low, allowing methods such as the histogram or the kernel density estimation [START_REF] Prouff | Theoretical and practical aspects of mutual information based side channel analysis[END_REF] to produce reliable results. However, this comes at the cost of sacrificing some, and maybe a large part, of the available information. MINE allows to directly use:

D new (k) = I(f (Z k ), L) (74)
as a distinguisher. This comes with two main advantages:

• Intermediate variables often leak at multiple instants in the trace. MINE allows to exploit all these leakage sources at the same time.

• Other intermediate variables, statistically dependent from the first one, can also leak information. For example, there could be some useful information about an AES key, before and after the application of the first S-box. In this context, MINE could exploit leakage from both intermediate variables at the same time.

Theorem 1 states that the optimal distinguisher is I(ϕ(Z k ), L) with ϕ being the leakage model. It is important to note that ϕ(Z k ) itself can be multidimensional. Therefore, an optimal MI attack should exploit this multidimensionality of the leakage model to increase the distinguishability of the correct hypothesis. However, it is frequent that multiple samples leak with the same underlying model: for example, a noisy version of the Hamming weight of Sbox[k * ⊕ P ] can leak multiple times in the trace. In such a context, the deterministic parts of the leakage of all these samples are all bijectively related. As adding bijection of the same variables multiple times would not change the MI, one can keep only one version of each different sub-leakage model. For example, if the target leaks (maybe multiple times) the Hamming weight of the first S-box of an AES and the Hamming distance between the S-box and k ⊕ P , Z k could be defined as k ⊕ P and one could replace ϕ(Z k ) by the two-dimensional vector:

HW Sbox[Z k ] , HW Sbox[Z k ] ⊕ Z k (75)
Remark 3. In practice, one may deliberately drop some intermediates variables for not being enough discriminating for wrong key candidates making them less tolerant regarding errors in the estimation of ϕ. For example, it is theoretically possible to use leakage on a xor: HW (k ⊕ P ) (assuming a Hamming weight a priori) but it is preferable to use intermediate variables where each bit depends on multiple bits of k such as the output of an S-box. Indeed, these variables are more discriminating since single bit errors on k are diffused to the whole variable which prevents from rewarding wrong hypotheses with several correct bits.

Scalability with masking order. In the context of masking, another advantage of multidimensionality emerges. In a classical d-order attack one often does not know the exact leakage time of each share, and therefore, has to compute the value of the distinguisher for each possible tuple (i 1 , . . . , i d ) and select the maximum. In the case of MIA the old distinguisher takes the following form:

D old (k) = max i1,...,i d {I(f (Z k ), L[i 1 , . . . , i d ])} (76)
For long traces, this can become a huge constraint since the total number of tuples grows exponentially with the masking order. Our version of the MIA which uses I(f (Z k ), L) as distinguisher, does not suffer from this since the complexity of the attack stays constant whatever the masking order. Note that it does not mean that masking is useless: it still decreases exponentially the information contains in side-channel traces [?] and an attack may require more traces to succeed. However, for a fixed number of traces, the computational effort required to mount a NEMIA is constant whatever the order of the attack.

Attack description

A step-by-step description of the NEMIA is given hereafter. It takes as input a set of traces and outputs a ranking of the key hypotheses.

1. Define an a priori φ on the leakage model. It can be multidimensional if multiple intermediate variables related to the key leak information. Also, a single intermediate variable can have different leakage models at different times. The test described in section 2.6 can be used to detect wrong a priori. Even if MIA is tolerant regarding estimation errors on ϕ, better a priori lead to more efficient attacks.

2. Compute, for all k, the hypothesis vectors:

H k = φ(Z k ).
3. Remove as much uninformative parts from traces as possible, using knowledge on the underlying algorithm (approximate leakage time windows). This step is not mathematically mandatory since uninformative components do not affect MI. However, a more compact leakage variable L, helps for the machine learning part of the attack mitigating the overfitting problem.

4. Compute I(H k , L), for all k, with MINE. Each estimation is the supremum of a moving average along the epochs of the validation loss function.

5. Rank the key hypotheses.

For masked implementation, the only step that changes is the construction of H k . If the shares have a Hamming weight leakage model, Th.2 proves that it is sound to use the Hamming weight of the corresponding unmasked intermediate variable in H k (one may do this for multiple intermediate variables). For a generic leakage model of the shares, the best strategy to adopt remains an open question. It appears that, in some cases, it is efficient to keep a restrictive number of bits of the unmasked variable as partition function, for example in a situation where some bits of the shares leak much more information than the others (an example of this is given in section 6).

Simulation experiments

In order to gain confidence in the mathematical results presented in this paper, as well as to gain intuition about their implications, this section presents simulation experiments.

The importance of the a priori

The main message of Th.1 is that, to maximize the distinguishability of the correct hypothesis, one should use the leakage model ϕ to create the hypothesis vectors H k . In a classical side-channel scenario, with no other specific information, one may often guess a Hamming weight leakage of the intermediate variables. This is justified by electronic arguments. However, it has been shown that bits may have different leakage behaviours, such as leakage weighting or even a sign inversion [START_REF] Cristiani | A Bit-Level Approach to Side Channel Based Disassembling[END_REF]. To illustrate Th.1, 10k synthetic traces leaking a slightly modified version ϕ 0 of the Hamming weight have been generated. They consist of a single sample leaking the Hamming weight of Z k * = Sbox(k * ⊕ P ) but with a flipped sign for bit 0 so that:

ϕ 0 (z) = -z 0 + 7 i=1 z i ( 77
)
with z i representing the i-th bit of z. Some Gaussian noise has been added to the traces so that L = ϕ 0 (Z) + N (0, 1). Fig. 1 shows the results of a NEMIA with k * = 0, both with HW and ϕ 0 as partition function. As predicted by Th.1, the distinguishability score:

S f = I f (Z k * ), L -max k =k * I f (Z k ), L (78)
is higher for f = ϕ 0 than for f = HW . Obviously, an attacker may not know ϕ 0 and an attack with the Hamming weight still succeeds in that case. However, this shows that, if by any means, an adversary knows the particularity of bit 0 of such a target, he can perform more efficient attacks. Semi-supervised attacks. This opens the idea of semi-supervised attacks. One of the main problems of profiling attacks is the portability [START_REF] Elaabid | Portability of templates[END_REF]. Indeed, during the characterization phase, the adversary learns a perfect representation of the leakage model which may overfit on the particular target which is profiled. It has been shown that portability to other targets is not trivial. Therefore NEMIA could be turned into a semi-supervised attack where the purpose of the characterization phase is only to learn general leakage characteristics, such as the sign or weighting of each bit, and use them as an improved a priori for a NEMIA. Since NEMIA is agnostic towards bijective errors in the leakage model estimation, it has a better chance of being portable on many other targets similar to the one used for profiling.

The potential of multidimensionality

One of the main advantages of NEMIA is its potential to exploit at the same time, multiple leakage sources. It is possible that multiple intermediate variables leak information on the key and each particular variable may leak multiple times in the traces.

To highlight this multidimensional potential, 10k synthetic traces have again been generated. Let suppose that a target running an AES both leaks A k * = HW (Sbox[k * ⊕ P ]) and B k * = HW (Sbox[k * ⊕ P ] ⊕ (k * ⊕ P )). One could imagine that the bus leaks the Hamming weight of the S-box data and that the update of the state register leaks the Hamming distance with its precedent value (e.g. [MEP + 08]). Let also suppose that the target leaks a noisy version of these models on 5 samples for each, so that the leakage variable L = ϕ(k * ⊕ P ) is equal to:

L = [A k * + N 1 , . . . , A k * + N 5 , B k * + N 6 , . . . , B k * + N 10 ] (79)
with N i = N (0, 50 2 ). An important amount of noise (σ = 50) has been added because the optimization brought by NEMIA is really highlighted in a context where the available information is very low. Fig. 2 shows the evolution of the average rank of k * for three attacks (each repeated 10 times with random sampling of the traces).

1. A classical CPA. It appears that for low numbers of traces, CPA is doing better than NEMIA but this is not very meaningful since attacks with guessing entropies greater than 50 (on a single key byte) are not really exploitable for a full key recovery. On the other hand, and as predicted by the theory, N EM IA F ull converges faster towards a ranking of 0 than N EM IA P artial , and both converge faster than CPA. To the best of our knowledge, classical MI-based attacks always performed worse than CPA in the literature, when considering the Hamming weight model, i.e. in the CPA playground. This experiment shows that in a low-information scenario (noisy traces), it may be worth considering NEMIA and its multidimensional advantages over classical attacks which treat samples one by one.

N EM IA

Empirical validation of theorem 2

Th.2 may seem very counterintuitive since it basically says that: when shares of a boolean masking leak in a Hamming weight model, one has:

I HW (Z k * ), L = I Z k * , L (80)
which is surprising since HW is highly non-injective and should at first glance, decrease the information. Corollary 2 says that this is even true when multiple samples leak a noised version of the Hamming weight of the shares. To verify this claim, 100k synthetic traces have been generated considering the following leakage:

L = HW (Z k * ⊕ M ) + N 1 , . . . , HW (Z k * ⊕ M ) + N 10 , HW (M ) + N 11 , . . . , HW (M ) + N 20 (81)
with Z k * = Sbox(k * ⊕ P ), N i = N (0, 1) and M being uniformly distributed in Z/256Z.

Fig. 3a shows the evolution of the loss function for both the HW and the identity function for the correct key hypothesis. As predicted, both converge towards the same value which confirms experimentally that the application of the HW does not alter information. The HW function is even doing a little better which can be explained with practical machine learning considerations. Indeed, the information being constant, it is easier for the network to learn with a 9-classes variable than with a 256 classes variable (note that in this experiment, id(Z k * ) has been encoded in binary rather than in OHE, because it produced slightly better results). Also, since overfitting was not really a problem in this experiment, the dropout parameter has been set to p = 0.1. Fig. 3b shows the result of the same experiment leaded on a second-order masking, with three shares and 10 leakage samples for each. Noise has been a bit decreased (σ = 0.5 instead of 1) to keep comparable level of information. The result sustains that Th.2 may be generalized to higher-order and that MINE is able to extract information even with a second-order masking. Since it is a masked implementation, the test described in remark 2 has first been conducted. Results are presented in Fig. 4a. I Z k * , L is more than four times greater than I HW (Z k * ), L which indicates that the underlying leakage of the shares is far from a pure Hamming weight model. In parallel to this, authors in [START_REF] Timon | Non-profiled deep learning-based side-channel attacks with sensitivity analysis[END_REF] presented a new methodology to use deep learning in the context of unsupervised attacks. Their method requires a partition function (which also has to be non-injective) and they reported that, for the ASCAD database, only keeping the value of the Least Significant Bit (LSB) produced better results than the Hamming weight without giving further explanations.

In a real attack scenario, an adversary mounting a NEMIA could obviously try to use every single bit of the unmasked variable as partition function. But in order to gain some intuition, and since the masks values are given in the database, we first performed a linear regression on both shares, assuming bits leak independently so that the actual leakage of share s is: 4b and4c show the evolution of the α i coefficients, on a leakage window for both shares. Since the implementation is protected by a boolean masking, a mono-bit leakage is exploitable only if it is present on the same bit of both shares. Out of the 8 bits, bit 0 (LSB) is clearly the one that leaks the most information since its coefficients are among the greatest ones in both shares. Thus, we computed with MINE I(Z k * [0], L) where Z k * [0] represents the LSB of Sbox(k * [3] ⊕ P [3]). It returned 0.09 bit, which is two times more than the information left with the Hamming weight (see Fig. 4a). This indicates that the LSB may be a good partition function since it is highly non-injective and still keep a decent amount of information for the correct hypothesis. We also tried with other bits but the information, while being non-zero, was significantly lower. Even though attacks with the Hamming weight were successful, we decided to use the LSB as partition function for the rest of our analysis. The attacks presented in this section uses traces that corresponds to a concatenation of two leakage window of size 100 per share. First window corresponds to samples between 140 and 240 while second window corresponds to samples 500 to 600. Also, we used a slightly different network architecture for NEMIA with convolutional layers applied to the traces (see Appendix D) which helps the network with feature extraction even though traces are not jittered.

7 i=0 α i s i + β. Figs.
Fig. 5a shows the guessing entropies for both a second order CPA and NEMIA. Attacks has been repeated 10 times to average the position of k * . For a number of traces greater than 800 the correct key k * was always the one predicted by NEMIA while the rank of k * decreases much slower for CPA. This Note that due to the low number of traces, the batch size used by the network has been set to 100 instead of 1000 for this experiment.

Pushing towards the limit. In order to evaluate the potential of NEMIA to exploit leakage even in very low information context, the dataset has been artificially degraded adding Gaussian noise N (0, σ 2 ) to each sample. Both NEMIA and second order CPA have been performed with σ going from 0 to 20, using the whole 50k traces. For each level of noise, the attacks have been repeated 10 times in order to compute the average rank of the correct hypothesis. Results are presented in Fig. 5b. They confirm that NEMIA is able to succeed even in situations where a classical higher order CPA would not.

Conclusion and perspectives

This paper first proposes a clarification of the state of the art around the MIA. It provides rigorous proofs whose goal is to derive the optimal MI-based attack working with high-dimensional traces. Combined with recent breakthroughs on neural MI estimation techniques, this allows to mount a new attack: the NEMIA, which benefits from both the strength of deep learning and information theory. Being able to exploit at the same time multiple leakage sources, it pushes the amount of effectively used information (depending on the strength of the attacker a priori) closer to the actual existing information between traces and secret. Simulations and real case experiments are presented to support the mathematical theory developed in this paper. They also show that NEMIA outperforms classical uni/bi-variate side-channel attacks and that this strategy may be worth to consider in low-information/high-noise situations, where all (or a large part of) the the available information contained in traces need to be used to mount a successful attack.

Several lines of research emerge from this paper. The mathematical analysis could be further extended, especially in the context of masking, in order to develop strategies for generic leakage model of the shares or for other masking schemes such as arithmetic masking. On the practical side, integrating the latest optimization on neural estimation techniques, as well as deep learning research on optimal networks architecture and hyperparameters would allow to mount more efficient attacks, taking as input larger portion of the traces, leading to better/easier attacks.

A Proof of lemma 1

Lemma 1. Let f : Z → R n be any function. For any leakage model ϕ: Z → R n there exists a decomposition of f into f = f 2 • f 1 , with f 1 : Z → N, f 2 : N → R n , satisfying the two following properties:

1) ∃ f 3 : Im f 1 → R n such that f 3 • f 1 = ϕ 2) ∀z ∈ Z, f 2 | f1 ϕ -1 ({ϕ(z)}) is bijective of reciprocal f -1 2 | f2•f1 ϕ -1 ({ϕ(z)})
Proof. Let us create a partition of Z = n i=1 P i where two elements z 1 , z 2 ∈ Z are in the same P i if and only if:

• ϕ(z 1 ) = ϕ(z 2 ) • f (z 1 ) = f (z 2 )
Then, one may define f 1 as f 1 (z) = i, ∀z ∈ P i . Since f 1 only collides for z that already collides through ϕ, there exists f 3 such that f 3 • f 1 = ϕ. As f is constant on P i , let us denote by v i its output on elements of P i . Then f 2 can be defined as f 2 (i) = v i so that f 2 • f 1 = f . Now let us prove 2). Let z ∈ Z and a, b ∈ f 1 (ϕ -1 ({ϕ(z)})) such that f 2 (a) = f 2 (b). There exists z a and z b such that a = f 1 (z a ) and b = f 1 (z b ) with ϕ(z a ) = ϕ(z b ) = ϕ(z). So:

• ϕ(z a ) = ϕ(z b ) • f 2 (f 1 (z a )) = f 2 (f 1 (z b )) ⇐⇒ f (z a ) = f (z b )
which means that z a and z b are in the same P i and thus collides through f 1 . So a = b which proves that f 2 | f1(ϕ -1 ({ϕ(z)})) is injective. Then, considering its set of destination being its image, one can say that this function is bijective with reciprocal function: f -1 2 | f2•f1(ϕ -1 ({ϕ(z)})) .

B Proof of corollary 1

Definition 1. A function f is said wider-than g if there exists another function h such that: h • f = g.

Corollary 1. Let L be defined as in (47). Then, for any function h wider than HW , S HW ≥ Sh.

Proof. There exists h such that h • h = HW . So:

S HW = I HW (Z k * ), L -max k =k * I HW (Z k ), L = I h • h(Z k * ), L -max k =k * I h • h(Z k ), L (82)
Since removing h in the second term can only increase the information:

S HW ≥ I h • h(Z k * ), L -max k =k * I h(Z k ), L (83)
By Th.2, HW maximizes over g the quantity: I g(Z k * ), L , so removing h in the first term cannot increase the information: Proof. The data processing inequality [START_REF] Beaudry | An intuitive proof of the data processing inequality[END_REF] ensures that applying f to any variables can not increase its mutual information with another variable so:

S HW ≥ I h(Z k * ), L -max k =k * I h(Z k ), L S HW ≥ Sh

D Network arhitecture for Ascad

Fig. 6 shows the network architecture used for the last ASCAD experiment. Note that when using convolutional layers with MINE, the convolutional layers should only be applied to the trace variable and not to f (Z k ) which would not make sense.

  -HW (z) HW (m) -n (57) Now let us evaluate P (l | HW (z)): P (l | HW (z)) = P (HW (m)) • A P (HW (z ⊕ m) | HW (z) and HW (m)) t. HW (z )=HW (z) P (z | HW (z)) • P (HW (z ⊕ m) | z and HW (m))

Figure 1 :

 1 Figure 1: I f (Z k L in terms of k, with k * = 0

Figure 2 :

 2 Figure 2: Guessing entropies for CPA, N EM IA P artial and N EM IA F ull

Figure 3 :

 3 Figure 3: Comparison of I Z k * , L and I HW (Z k * ), L on masked synthetic traces

  Sbox(k * [3] ⊕ P [3]) ⊕ r[3]

Figure 4 :Figure 5 :

 45 Figure 4: Analysis of the ASCAD leakage model: a) Test from remark 2 -b) & c) Coeficients of a linear regression on the given variable

  Let A and B be a two discrete random variables. Let f : A → R n be any function. Then:H f (A) | B ≤ H(A | B) (85)

I

  f (A), f (A) | B ≤ I(A, A | B) H f (A) | B ≤ H(A | B) (86)Lemma 3. Let A and B be a two discrete random variables. Let f : A → R n be any function. Then:H A | f (B) ≥ H(A | B) (87)Proof. Again, using the data processing inequality[START_REF] Beaudry | An intuitive proof of the data processing inequality[END_REF]:I A | f (B) ≤ I(A, B) H(A) -H A | f (B) ≤ H(A) -H(A | B) H A | f (B) ≥ H(A | B) (88)

Figure 6 :

 6 Figure 6: Network architecture for ASCAD experiment

Due to MI estimator limitation D(k) is often replaced in practice by maxi I(f (Z k ), L[i]), where L[i]represents the i-th sample of the trace. This does not affect the theory described in this section so we decided to keep it as described in eq. 5 for the sake of simplicity. More details are provided in section 4.2.

Note that other criteria such as the distance with the mean of the wrong hypotheses can also be used as discussed in remark 1.


Nils Laurent
email: nils.laurent1@univ-grenoble-alpes.fr

Sylvain Meignen
email: sylvain.meignen@univ-grenoble-alpes.fr

Julie Fontecave-Jallon

Bertrand Rivet
email: bertrand.rivet@gipsa-lab.grenoble-inp.fr

A Novel Algorithm for Heart Rate Estimation Based on Synchrosqueezing Transform

Keywords: HR estimation, non-negative matrix factorization, Wasserstein distance, synchrosqueezing transform

To estimate the heart rate (HR) from electrocardiogram (ECG), time-frequency representations such as the short-time Fourier transform (STFT) is often used. As the STFT is constrained by the choice of a specific analysis window used for its definition, we alternatively propose to estimate HR from a synchrosqueezed STFT. More precisely, we build a novel algorithm inspired by non-negative matrix factorization to estimate HR by determining the minimal Wasserstein distance between a synchrosqueezed STFT and columns of a specific dictionary matrix. Throughout numerical simulations carried out on both synthetic and real ECGs, we show in what way to use a synchrosqueezed STFT rather than STFT improves HR estimation.

I. INTRODUCTION

HR estimation from electrocardiogram (ECG) recordings has been a very active research topic in the past decades, and can be either signal-based, as for instance in techniques using R-peak detection [START_REF] Pan | A real-time qrs detection algorithm[END_REF], [START_REF] Daniel T Kaplan | Simultaneous qrs detection and feature extraction using simple matched filter basis functions[END_REF], or time-frequency (TF) based [START_REF] Huang | Ecg arrhythmia classification using stft-based spectrogram and convolutional neural network[END_REF], [START_REF] Shah Shaik | A method for qrs delineation based on stft using adaptive threshold[END_REF]. TF analysis has been used on ECG signals in many different circumstances as for instance to detect abnormal behaviors like sleep apnea [START_REF] Af Quiceno-Manrique | Detection of obstructive sleep apnea in ecg recordings using time-frequency distributions and dynamic features[END_REF], to study arrhythmia [START_REF] Sultan | Ecg arrhythmia classification using time frequency distribution techniques[END_REF], to identify coronary artery deseases [START_REF] Sharma | A new method to identify coronary artery disease with ecg signals and time-frequency concentrated antisymmetric biorthogonal wavelet filter bank[END_REF] or to separate fetal ECG from that of the mother [START_REF] Su | Extract fetal ecg from single-lead abdominal ecg by de-shape short time fourier transform and nonlocal median[END_REF]. Regarding this last application, there exist many other alternative techniques to extract fetal ECG among which one may cite blind source separation [START_REF] De Lathauwer | Fetal electrocardiogram extraction by blind source subspace separation[END_REF], [START_REF] Akhbari | Fetal electrocardiogram r-peak detection using robust tensor decomposition and extended kalman filtering[END_REF], blind adaptive filtering [START_REF] Graupe | Blind adaptive filtering for non-invasive extraction of the fetal electrocardiogram and its nonstationarities[END_REF], Kalman filter [START_REF] Niknazar | Fetal ecg extraction by extended state kalman filtering based on single-channel recordings[END_REF], to name but a few.

In the present paper, the focus is put on TF-based HR estimation and, more precisely, we investigate the STFT of ECGs. In that context, the basic principle ruling HR estimation is that an ECG corresponds at each time instant to the sum of several harmonics whose frequencies are multiples of the so-called fundamental frequency (FF) which is the sought HR. For its estimation, we first propose to compare STFT at each time instant with the columns of a dictionary matrix, being aware that such an approach will necessarily be constrained by the choice of the window used in STFT definition. To cope with this issue, we carry out the same approach but using the reassigned STFT named synchrosqueezing transform (FSST) instead of STFT. The rationale to consider such a time-frequency representation (TFR) is that we believe the reassignment process, by sharpening the TFR it is based on, should lead to a more accurate HR estimation. Note that the synchrosqueezing transform was first introduced to reassign the continuous wavelet transform [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF] and then extended to STFT reassignment in [START_REF] Oberlin | The Fourier-based synchrosqueezing transform[END_REF], [START_REF] Thakur | The synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications[END_REF] . As FSST reassigns exactly only signals made of pure harmonics, an extension of FSST was introduced to deal with more frequency modulated signals [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF], called FSST2. In this regard, we will therefore investigate whether to use FSST2 instead of FSST also brings some improvements in terms of HR estimation.

The paper is organized as follows: in Section II, we recall some basic definitions on STFT, synchrosqueezing transforms, and non-negative matrix factorization (NMF) and then detail our novel algorithm, inspired by NMF, that uses either STFT or its reassigned versions to perform HR estimation, in Section III. Numerical results on both simulated and real ECG signals conclude the paper.

II. DEFINITIONS A. Time frequency Analysis

Considering a signal f ∈ L 1 (R) ∩ L 2 (R) and a real window g ∈ L ∞ (R) ∩ L 2 (R), its STFT is defined as: V g f (t, η) = R f (τ )g(τ -t)e -2iπ(τ -t)η dτ. (1)
This work was supported in part by the University Grenoble Alpes under IRS Grant "AMUSETE" and the ANR ASCETE project with grant number ANR-19-CE48-0001-01.

In this paper, we assume ECG is modeled by a specific kind of multicomponent signals (MCSs):

f (t) = P p=1 A p (t)e 2iπpφ (t) , (2)
where φ is the fundamental frequency (FF) and A p the instantaneous amplitude (IA).

B. Synchrosqueezing Transforms

While STFT is very commonly used to analyze signals of type (2), it is constrained by the Heisenberg uncertainty principle that stipulates that STFT cannot be perfectly localized in time and frequency at the same time. Therefore, the choice of the analysis window g is central in that context.

To improve the TF localization of STFT, different reassignment techniques were introduced among which FSST [START_REF] Oberlin | The Fourier-based synchrosqueezing transform[END_REF], [START_REF] Thakur | The synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications[END_REF] is one of the most popular. In a nutshell, using the local instantaneous frequency estimator given by a reassignment operator

ω(t, η) = {ω(t, η)} with ω(t, η) = η - 1 2iπ V g f (t, η) V g f (t, η)
,

where {X} is the real part of complex number X. FSST reassigns STFT coefficients with modulus larger than a threshold λ through:

T λ f (t, ω) = |V g f (t,η)|>λ V g f (t, η)δ(ω -ω(t, η))dη. (3)
A discrete time and frequency version of the synchrosqueezing transform given by (3) can then be derived, corresponding to a K × N matrix, N being the length of the signal and K the number of frequency bins. Since FSST is unable to properly reassign MCSs containing frequency modulated modes, an extension was proposed to deal with this issue through the so-called second order synchrosqueezing transform (FSST2) [START_REF] Oberlin | The Fourier-based synchrosqueezing transform[END_REF], [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF]. The latter is defined by means of the second order local complex IF estimate of f [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF]:

ω(2) (t, η) = ω(t, η) + q(t, η)(t -t(t, η)) if ∂ t t(t, η) = 0 ω(t, η) otherwise, in which: t(t, η) = t + V tg f (t,η) V g f (t,η) , and q(t, η) = 1 2iπ V g f (t,η)V g f (t,η)-(V g f (t,η)) 2 V tg f (t,η)V g f (t,η)-V tg f (t,η)V g f (t,η)
. Finally, ω (2) (t, η) = {ω (2) (t, η)} corresponds to the desired local IF estimate. FSST2 is then defined by replacing ω by ω (2) in (3). In the same way as for FSST, a discrete time and frequency version of FSST2 can be derived.

C. NMF Basics

Non-negative matrix factorization (NMF) decomposes a given non-negative data matrix X ∈ R K×N into two non-negative matrices, the dictionary matrix W ∈ R K×Ke and the activation matrix H ∈ R Ke×N such that X ≈ WH, K e being the number of components in the dictionary. The decomposition is based on minimizing the reconstruction error of X through WH, which can be formulated as [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF] min

W,H D(X|WH) subject to W ≥ 0, H ≥ 0. (4)
The most popular cost functions D are the Euclidean distance, Kullback-Leibler (KL) divergence and Itakura-Saito (IS) distance.

To solve the problem defined in (4), a maximization-minimization algorithm based on multiplicative updates [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF] is widely used [START_REF] Leplat | Blind audio source separation with minimum-volume beta-divergence NMF[END_REF].

D. Definition of the Dictionary Matrices

In the context of ECG signals analysis, we consider a fixed dictionary W, in which each column corresponds to a specific fondamental frequency (FF). In the following study, we will consider that X is either STFT (see an illustration in Fig. 1 (a)), FSST or FSST2 moduli (see an illustration in Fig. 1 (b) for FSST2). Taking into account that, for resting patients, FF may vary between 30 and 180 beats per minute (bpms), when X is either FSST or FSST2 moduli, we consider the dictionary matrix given in Fig. 1 (d), which we denote by W 0 . This matrix consists of K e columns corresponding to 151 components each of which corresponds to an integer FF between 30 and 180 bpms. The coefficients in each column are set to one when they correspond to multiples of FF, and zero elsewhere. With such a dictionary matrix, only some of the harmonics are taken into account, since fewer harmonics are present in W 0 when FF increases. Furthermore, the fact that the different harmonics have different amplitude is not considered. We will see that these aspects do not much alter HR estimation performance with the algorithm we propose. Dictionary matrix W 0 used with FSST2-based approach. For the sake of legibility, in (c) and (d) we do not normalize the column of the matrices.

When X corresponds to STFT modulus, one has to take into account the width of the analysis window in the definition of the dictionary matrix. For the sake of simplicity, we here consider a Gaussian window of the form e -π x 2 σ 2 . In the Fourier domain, it corresponds to e -πσ 2 η 2

. So, to build the dictionary in the STFT case, we consider a simple convolution of each column of W 0 with e -πσ 2 η 2 . This results in the dictionary matrix denoted by W 1 (see Fig. 1 (c)).

III. NMF-LIKE DECOMPOSITION BASED ON FSST MAGNITUDE A. Computation of Activation Matrix Using Wasserstein Distance

With a fixed dictionary matrix, the computation of H when D is the Euclidian distance boils down to a simple least-square minimization problem:

H = argmin H0 X -WH 0 2 ,
and thus H satisfies W T WH = W T X, and is not necessarily unique. Alternatively, we impose that H be associated with only one function in the dictionary, and we compute H by minimizing, at each time instant, a specific earth mover distance (EMD) [START_REF] Peleg | Fast and robust earth mover's distances[END_REF]. Let us recalll that EMD is a sliced Wasserstein distance that aims at comparing probability distributions, and that it has already been used in the TF context for instance in [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF], [START_REF] Daubechies | Conceft: concentration of frequency and time via a multitapered synchrosqueezed transform[END_REF]- [START_REF] Duong | High-order synchrosqueezing transform for multicomponent signals analysis-with an application to gravitational-wave signal[END_REF].

We propose to compute the sliced Wasserstein distance d between X at time n, corresponding to the column X n of X (normalized by its sum so that it can be viewed as a probability distribution) and each column of the matrix W, the i th column of W being denoted by W i (normalized by its sum), and then define

i 0 = argmin i d(W i , X n ). (5)
Finally, H n , the n th column of H, is such that H i,n = 0 if i = i 0 , and H i0,n = 1.

B. Algorithm for HR Computation

The approach we propose for HR estimation is based on minimal Wasserstein and consists of several steps. First, one defines a probability of false HR detection on a test signal. To do so, we first estimate the mean HR i M by the median of the component indices associated with the minimal Wasserstein distance over the whole test signal duration. Let us denote by îM this estimate. To measure good detections, we consider the interval Î = [ 3 îM 4 , 3 îM 2 ], and then define B as follows:

B = #{n, argmin i d(W i , X n ) ∈ Î} N . (6)
The rationale supporting the definition of Î is that the largest sub-harmonic is located at îM 2 and the smallest harmonic at 2 îM , and thus Î corresponds to the half intervals between these and FF. The probability of false detection is then defined as

P f = 1 -B.
Assuming P f is known, the first step of the proposed HR estimation for a given ECG signal consists of computing n med an estimate of the minimal time index such that the probability that the median of the detection computed from time indices 1 to n med is exact with probability p 0 . Assuming the locations of false detections are randomly distributed, n med can be computed using the binomial law as follows:

n med = argmin n ( n/2 k=0 C k n P k f (1 -P f ) n-k ≤ p 0 ). (7)
Knowing n med , we compute HR estimates for time indices between 1 and n med , following the scheme described in Algorithm 1. To check a posteriori the validity of the computed n med , we compute the median of the detection between time indices 1 and n med , and then recompute an estimate of B from this set using the same approach as before. If B is lesser than 1 -βP f (β is fixed to 3/4 here) n med is increased until B becomes larger than 1 -βP f . Having defined HR estimation for the first n med time indices, the estimation carries on following the approach described in Algorithm 2, in which γ is some parameter whose influence will be discussed in the Results section.

Algorithm 1: HR estimation step 1 It is worth noting here that to improve the computational efficiency of the proposed algorithm we compute only the Wasserstein distance on a small number of elements of the dictionary.

Input : W,X, n med for n = 1, • • • , n med do ĩ[n] = argmin i∈[1,Ke] d(W i , X n ) end i M = median ĩ[n], n = 1, • • • , n med for n = 1, • • • , n med do î[n] = argmin i∈[ 3 4 i M , 3 2 i M ] d(W i , X n ) end

IV. RESULTS

In this section, we first investigate false detections depending on the TFR used to compute B, then we explain how to fix the parameter γ in Algorithm 2 by studying the behavior of this algorithm on synthetic ECG generated by the MATLAB function Algorithm 2: HR estimation step 2 "ecgsyn" developed in [START_REF] Patrick E Mcsharry | A dynamical model for generating synthetic electrocardiogram signals[END_REF]. Finally, we conclude by showing some results on real ECGs extracted from the SiSEC database [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF].

Input : W,X, n med , î[1, • • • , n med ], γ for n = n med + 1, • • • , L do ∆ n = std î[k], k = 1, • • • , n -1 î[n] = argmin [î-γ∆n,î+γ∆n] d(W i , X n ) end HR Mean B STFT B FSST B FSST2

A. Determination of the Probability of False Detection

To compare the different probability of false detection depending on the type of TFR used to define X, we consider synthetic ECGs generated as in [START_REF] Patrick E Mcsharry | A dynamical model for generating synthetic electrocardiogram signals[END_REF], compute the corresponding STFT, FSST, and FSST2 and then (5) for each time instant when X is associated with one of the just mentioned TFRs. With these synthetic signals, the mean HR component i M is known, and we measure good detection by considering I = [ 3i M 4 , 3i M 2 ] and then

B = #{n, argmin i d(W i , X n ) ∈ I} N , (8)
the probability of false detection being 1 -B. We numerical notice that, with STFT or FSST2 the minimum of the Wasserstein distance is much more likely to be associated with a sub-harmonic than when FSST is used, resulting in a much higher probability of false detection, and this is true whatever the mean HR ( the results are depicted in Table. I). An illustration of this is given in Fig. 2 (a)-(c) when the TFR are either STFT, FSST and FSST2, respectively: the white curves correspond to the indices associated with the minimal Wasserstein distance for each time index n. We notice that false detections usually occur when HR changes rapidly, and, as the signal is quasi-harmonic, to make a linear chirp approximation on the modes, as is done in FFST2, results in error at these time indices. As our plan is to use Wasserstein distance for HR computation, we will make use of FSST rather than STFT or FSST2.

B. Determination of γ on Synthetic ECG The goal of this section is to define an appropriate choice for parameter γ used in Algorithm 2. For that purpose, we investigate how good this algorithm is to estimate the mean and standard deviation of HR, which are input parameters of the function "ecgsyn". In what follows, we denote by HR syn the mean heart rate and by σ syn its standard deviation. To make HR estimation realistic we consider HR syn = 60 bpms (corresponding to component index 31) and σ syn = 2 bpms. Such a signal enables a clear separation between harmonics and enough variability of HR, so that its study is useful for the analysis of real ECGs that comes next. For that type of signals, we compute the mean of î along with its standard deviation, expressed in terms of a number of components, both as a function of γ. The results depicted in Fig. 3 (a) and (b) show that the proposed algorithm enables very accurate estimations of both the mean and standard deviation of HR provided γ is taken around 2.5. Note that the range for γ to retrieve the mean and standard deviation of HR is very similar when either FSST or FSST2 is considered. We here study ECGs extracted from SiSEC database [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF], which correspond to thoracic recordings, which can be contaminated by different types of noises, and in particular impulse noises. EMD based on FSST computed on such a signal, and γ fixed as explained in Section IV-B, is depicted in Figure 4 (a). The presence of impulsive noises tends to create some local minima of the Wasserstein distance at frequencies corresponding to the largest sub-harmonic of the true FF. As FSST somehow reduces the noise while enhances the significant harmonics, the proposed HR estimation algorithm never gets trapped by minima that arise at sub-harmonic locations. We display in Fig. 4 (b) î (white line) along with the TF region corresponding to standard deviation ∆ n , when n varies.

C. Study of Real ECGs

V. CONCLUSION

In this paper, we have introduced a novel procedure for HR estimation from the reassigned time-frequency representation of ECG recordings. The proposed approach has been developed on the time-frequency representation associated with the synchrosqueezing transform of the short-time Fourier transform. The simulations carried out on synthetic ECGs help us finding the most appropriate TFR for HR estimation and the right parameters values. An application to real ECG signals was then proposed, for which promising results were derived. In a near future, we plan to extend this approach to more complex situations principally the monitoring of fetal activity, to study more in detail the impact of noise on the proposed HR estimate, and also to investigate how to use the proposed algorithm for abnormal behavior detection like arrythmia or sleep apnea.

Fig. 1 .

 1 Fig. 1. (a): STFT of an ECG generated following [19]; (b): FSST2 of the same ECG; (c): Dictionary matrix W 1 used with STFT-based approach; (d):Dictionary matrix W 0 used with FSST2-based approach. For the sake of legibility, in (c) and (d) we do not normalize the column of the matrices.

Fig. 2 .

 2 Fig. 2. (a): modulus of STFT, the white curve corresponds to the indices associated with the minimal Wasserstein distance ( B = 39.61%); (b):same as (a) but for FSST ( B = 94.34%) (c):same as (a) but for FSST2 ( B = 75.4%).

Fig. 3 .

 3 Fig. 3. (a): Estimation of HR mean from î. (b): Estimation of the corresponding standard deviation.

Fig. 4 .

 4 Fig. 4. (a) FSST based EMD of a real ECG signal. (b) HR detection on the same signal using γ = 3.

TABLE I COMPUTATION

 I OF B ON SYNTETHIC ECG SIGNALS WITH DIFFERENT HR MEANS WHEN EITHER STFT, FSST OR FSST2 ARE USED TO COMPUTE X