Juliette Legrand
email: juliette.legrand1@inrae.fr

Philippe Naveau

Marco Oesting

Evaluation of binary classifiers for asymptotically dependent and independent extremes

Keywords: Multivariate Extreme Value Theory, Binary Classification, Multivariate Regular Variation, Hidden Regular Variation, River Discharges

Machine learning classification methods usually assume that all possible classes are sufficiently present within the training set. Due to their inherent rarities, extreme events are always under-represented and classifiers tailored for predicting extremes need to be carefully designed to handle this under-representation. In this paper, we address the question of how to assess and compare classifiers with respect to their capacity to capture extreme occurrences. This is also related to the topic of scoring rules used in forecasting literature. In this context, we propose and study a risk function adapted to extremal classifiers. The inferential properties of our empirical risk estimator are derived under the framework of multivariate regular variation and hidden regular variation. A simulation study compares different classifiers and indicates their performance with respect to our risk function. To conclude, we apply our framework to the analysis of extreme river discharges in the Danube river basin. The application compares different predictive algorithms and test their capacity at forecasting river discharges from other river stations. As a byproduct, we study the special class of linear classifiers, show that the optimization of our risk function leads to a consistent solution and we identify the explanatory variables that contribute the most to extremal behavior.

Introduction

In binary classification, one typically considers data of the form (X, Y ) where Y ∈ {-1, 1} represents a binary response to the input X ∈ [0, ∞) d . In this paper, we focus on the case that Y = Y (u) represents the occurrence of an extreme event, Y (u) = 1 indicating that a random quantity H crosses a level u, called the threshold, and Y (u) = -1 otherwise, that is

Y (u) = +1, if H > u, -1, otherwise. (1)
In the following, for simplicity, we focus on the case that H is a non-negative random variable such that P(H > u) > 0 for all u > 0 and its upper end point is infinite.

In extreme value analysis, one is interested in the behavior of Y (u) for high levels, that is for u → ∞, and, therefore, also any classifier g needs to be adapted to the threshold u. Thus, for every u > 0, let g(•; u) be a measurable function from R d to {-1, 1}. In order to evaluate the quality of the classification at a certain level u, we consider a loss function l u that assigns a cost to a classifier g(•; u) and a realization (x, y (u) ). Here, it is important to note that, by definition of rare events, P(Y (u) = 1) is very small and P(Y (u) = -1) is close to one as u gets large. This imbalance can lead to atypical and/or undesirable comparisons of classifiers. For example, the "always optimistic" classifier that never forecasts an extreme can be defined as g(X; u) ≡ -1, almost surely. To see how to handle this naive classifier, the classical risk function defined as the expectation of the indicator 1{g(X; u) = Y (u) } can be written as P g(X; u) = Y (u) = P(H > u, g(X; u) = -1) + P(H ≤ u, g(X; u) = 1).

(2)

If g(X; u) ≡ -1, then P(g(X; u) = Y (u) ) = P(H > u) goes to zero as u gets large. Hence, the classical risk function E(1{g(X; u) = Y (u) }) will systematically favor the always optimistic classifier for extremes. To avoid this undesirable feature, the loss function has to be modified. One natural idea is to re-scale by P(H > u) and introduce the loss function 1{g(X; u) = Y (u) }/ P(H > u). In this case, the risk, i.e., the expected loss E(1{g(X; u) = Y (u) }/ P(H > u)), goes towards one as u gets large.

Another trivial but also interesting case is the "crying wolf" forecaster who always issues g(X; u) ≡ +1, see also the forecaster's dilemma (e.g. [START_REF] Lerch | Forecaster's Dilemma: Extreme Events and Forecast Evaluation[END_REF]. In this case, Equation (2) implies that P g(X; u) = Y (u) = P(H ≤ u) and, consequently, the risk E(1{g(X; u) = Y (u) }/ P(H > u)) goes towards infinity as u gets large. This limiting cost indicates that the "crying wolf" forecaster is much worse than the overly optimistic one. Both of them are unreasonable in practice and there is no reason to strongly favour one over the other one. For this reason, we propose to use a following weighted loss function l u (g; (x, y)) = 1 P(Y (u) = 1 or g(X; u) = 1) 1{g(x; u) = y} and the associated risk R (u) (g) = E(l u (g; (X, Y (u) )) = P(g(X; u) = Y (u) ) P(Y (u) = 1 or g(X; u) = 1) .

(3)

By construction, the event {g(X; u) = Y (u) } implies that {Y (u) = 1} or {g(X; u) = 1} and therefore, necessarily, R (u) (g) ∈ [0, 1]. In particular, the naive classifier g(X; u) ≡ -1 possesses unit risk with R (u) (g) = 1 at each level u > 0. Similarly, the risk of the "crying wolf" classifier g(X; u) ≡ +1 is then equal to R (u) (g) = P(H ≤ u) and converges to one as u → ∞.

Hence, the value of one is reached by the two worst cases scenarios in terms of classifiers. This unit value provides a clear benchmark that can be compared to any other classifier satisfying the existence of the limit

R(g) = lim u→∞ R (u) (g) ∈ [0, 1].
We call such classifiers extremal.

In the weather forecast literature (e.g. [START_REF] Schaefer | The critical success index as an indicator of warning skill[END_REF], the definition of R (u) (g) can be linked to the critical success index, also called the threat score. The critical success index computes the total number of correct event forecasts (hits) divided by the total number of forecasts plus the number of misses (hits + false alarms + misses). Hence, 1 -R (u) (g) can be understood as a critical success index for extremes. In the context of rare events forecasts, [START_REF] Stephenson | The extreme dependency score: a non-vanishing measure for forecasts of rare events[END_REF] highlighted some advantages and drawbacks of various risk functions, including the critical success index. In particular, these authors linked forecast scoring rules with two dependence indices used in EVT (see [START_REF] Coles | Dependence measures for extreme value analyses[END_REF]

χ = lim u→1 P(U > u | V > u) and χ = lim u→1 2 log(P(U > u) log(P(U > u, V > u)) -1 ,
where the two random variables U and V follow the same continuous uniform distribution on [0, 1]. The choice of uniform marginals can be made whenever the forecast can be assumed to be calibrated, i.e. observations and forecasts follow the same marginal distributions and can be transformed into uniforms. Concerning the extremal dependence strength between U and V , if χ > 0, then the variables U and V are said to be asymptotically dependent and χ = 1. If χ = 0, then the variables U and V are said to be asymptotically independent and χ < 1 captures some second order extremal dependence information. [START_REF] Stephenson | The extreme dependency score: a non-vanishing measure for forecasts of rare events[END_REF] advocated the use of χ and called it the extreme dependency score. Later on, [START_REF] Ferro | Extremal dependence indices: Improved verification measures for deterministic forecasts of rare binary events[END_REF] proposed two different scores and studied their properties. But the link with the concept of asymptotic independence was not clear and the convergence results of their estimators were not fully developed.

In contrast to χ, one drawback of χ is that its formula is not easy to explain to practitioners. In comparison, R(g) as a type of the critical success index can be interpreted with ease. Hence, it is of interest to extend this definition to the asymptotic independent case.

In the machine learning literature, there is a vast body of work about imbalanced data (see, e.g. [START_REF] Hilario | Learning from Imbalanced Data Sets[END_REF]. Within this context of imbalanced problems, R(g) can be understood as a function of precision and recall, two metrics well-used to score binary classifiers in learning research (see, e.g. He and Ma, 2013, and the F -score). Still, the review paper by [START_REF] Haixiang | Learning from class-imbalanced data: Review of methods and applications[END_REF] did not mention any imbalanced method based on extreme value theory (see, e.g. [START_REF] Resnick | A Probability Path[END_REF]. To our knowledge, very few theoretical links have been made to bridge the imbalanced learning and multivariate extreme value theory. One noticeable exception is [START_REF] Jalalzai | On Binary Classification in Extreme Regions[END_REF] who worked on binary classifiers for extremes under the regular variation hypothesis. But they did not focus on R (u) (g). Instead, they studied a different setting where the object of interest was P(g(X) = Y | X > u) where X represents a norm with u large. Hence, their conditioning event was { X > u}, while our conditioning depends on Y with the set {Y (u) = 1 or g(X; u) = 1}, see Equation (3). So, their interest was centered on the classifier performance when the norm of the explanatory vector X was large. Our focus is on large values of H in the production of extreme events of the type Y (u) = 1 when H > u, see Equation (1). [START_REF] Jalalzai | On Binary Classification in Extreme Regions[END_REF] provided various theoretical results based on the main assumption that the conditional distribution of X given Y = ±1 was regularly varying with an angular measure that depends on Y = ±1.

In this study, one part of our results is based on the concept of hidden regular variation (see, e.g. [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF][START_REF] Heffernan | Hidden regular variation and the rank transform[END_REF][START_REF] Ferro | A probability model for verifying deterministic forecasts of extreme events[END_REF]. In particular, we take advantage of the model of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] to derive the asymptotic properties of our estimators.

Our paper is organized as follows. In Section 2, we propose and study a risk function that can handle both the asymptotic dependent and independent cases. Estimators are also constructed and their asymptotic properties derived. Section 3 focuses on a simulation example that highlights the difficulty to compare common classifiers in the case of asymptotic independence. In Section 4, we revisit the well studied example of the Danube river application and see how the choice of the metric can change the ranking of classifiers. Note that, besides the proofs of all propositions, the appendix addresses the questions of how to optimize the linear classifier for extremes and how to choose the relevant features, see Section B.

Risk, upper tail equivalence and extremal dependence

The following lemma provides a flexible blueprint to link risk functions with probabilities based on general sets. We will apply it under different setups linked to extreme events. Lemma 1. Let A ε be a sequence of measurable sets of increasing sizes with decreasing ε ∈ [0, 1], in particular A 1 ⊆ A ε ⊆ A 0 . Let B ε be the same type of set sequence such that P(A 1 ∩ B 1 ) > 0. The following ratio R(A ε , B ε ) can be written as

R(A ε , B ε ) := P(A 1 B 1 | A ε ∩ B ε ) P(A 1 ∪ B 1 | A ε ∩ B ε ) = 1 - 1 P(B 1 |A 1 ∩ B ε ) + 1 P(A 1 | A ε ∩ B 1 ) -1 -1 , (4
) where A 1 B 1 denotes the symmetric difference. In addition, we have the three following properties for R:

(a) R(A ε , B ε ) is non-increasing in ε with R(A 1 , B 1 ) = 0.
(b) Let A ε be another sequence of measurable sets of increasing sizes with decreasing ε.

If

A 1 = A 1 and A ε ⊆ A ε for some ε ∈ [0, 1) then R(A ε , B ε ) ≤ R(A ε , B ε ).
(c) If for any ε ∈ [0, 1] and ε ∈ [0, 1], there exist some positive constants a and b and some positive function c ε,ε such that

P(A ε ∩ B ε ) = c ε,ε P a (A ε ) P b (B ε ), ( 5
) then R(A ε , B ε ) = 1 - c ε,1 c 1,1 (P(A 1 | A ε )) -a + c 1,ε c 1,1 (P(B 1 | B ε )) -b -1 -1 . ( 6
)
We deduce from Equation (4) that R(A ε , B ε ) = 0 if and only if

P(B 1 | A 1 ∩ B ε ) = P(A 1 | A ε ∩ B 1 ) = 1.
The second property of this lemma indicates that, for a given ε, the risk function becomes smaller if the set A ε is as small as possible.

Equation ( 5) can be viewed as a mixing condition that leads to a simple expression of R(A ε , B ε ) based on disjoint events.

We want to use the ratio (4) in order to generalize the risk R (u) (g) defined in (3). In the latter definition, adapted to extremes, the set {H > u} was considered and one could set B 1 = {H > u}, for instance, and consequently, B ε = {H > εu} would be a natural sequence if H is regularly varying. The choice of A ε is open and it can play an important role †. In the next paragraphs, we set ε = 0 and discuss the choice of A 1 and A 0 . In Section 2.1, we will work with ε > 0. In that case, we will call the corresponding risk the conditional risk.

In the rest of this paper, we restrict our attention on a particular form of classifiers

g(X; u) = +1, if g(X) > u, -1, otherwise, (7)
for some function g : R d → (0, ∞). The function g(.) does not have to be a norm. It can be understood as any projection/summary of the explanatory variables X onto the positive real line (see, e.g. [START_REF] Aghbalou | Tail inverse regression for dimension reduction with extreme response[END_REF], for projection techniques for extremes). This corresponds to the set A 1 = {g(X) > u} in Lemma 1. A special case of this lemma is to set ε = 0 and when A 0 and B 0 are equal to the full set, i.e. P(A 0 ) = P(B 0 ) = 1, and A 1 = {g(X) > u} and B 1 = {H > u}. In this case, Equation (4) tells us that R (u) (g) defined by (3) satisfies

R (u) (g) = 1 - 1 P(H > u | g(X) > u) + 1 P(g(X) > u | H > u) -1 -1
.

This leads to the following expression of R (u) (g)

R (u) (g) = 1 - P(g(X) > u | H > u) 1 -P(g(X) > u | H > u) + P(g(X) > u)/ P(H > u) . ( 8
)
Within the class defined by Equation ( 7), the effect of the marginal distributions of H and g(X) on R(g) ‡ can be explained. When g(X) has a lighter tail than H, i.e. P(g(X) > u)/ P(H > u) → 0, then we also have P(g(X) > u | H > u) → 0, and consequently R (u) (g) → 1, as u → ∞.

In the case where g(X) possesses a heavier tail than H, i.e. P(g(X) > u)/ P(H > u) → ∞, we can also show that

R (u) (g) → 1, as u → ∞.
This indicates that, whenever g(X) and H are not tail equivalent, the classifier g(X) cannot outperform naive classifiers with respect to R (u) (g) for large u. Hence, the marginal behavior of g(X) has a direct impact of its predicting capacity in terms of R(g). This is a widely known fact in forecast verification. In particular, a paradigm †Although we will apply Lemma 1 to sets A ε that are rare events, this is a not a necessity.

‡the shortcut notation R(g) corresponds to the case where g(X; u) = +1 is built from the event g(X) > u in the associated R(g).

promoted by Gneiting and his co-authors (see, e.g. [START_REF] Gneiting | Probabilistic forecasts, calibration and sharpness[END_REF] is to work only with calibrated forecasts. In our case, calibration means that g(X) and H have the same distributions, and consequently P(g(X) > u)/ P(H > u) = 1 for all u. In practice, it may be difficult to ensure that this constraint holds for extremes (see, e.g. [START_REF] Lerch | Forecaster's Dilemma: Extreme Events and Forecast Evaluation[END_REF]Taillardat et al., 2019). To illustrate this, suppose that P(X > u) = u -1 for all u ≥ 1, and the variable H, independently of the value of X, is either equal to δX or (2 -δ)X with probability .5 and the constant δ ∈ (0, 1). Then, we have

P(H > u) = 1 2 P (δX > u) + 1 2 P ((2 -δ)X > u) =      P(X > u), if u ≥ 2 -δ, 1 2 1 + δ u , if δ < u ≤ 2 -δ, 1, if δ ≥ u.
Hence, X and H are more than tail equivalent, they have identical tail behavior for large u. Concerning classifiers, linear ones of the type g a (X) = aX with a > 0 belong to the class defined by ( 7). They are tail equivalent to H, and R(g a ) < 1. Although X and H have identical tail behaviors, the choice of a = 1 is not optimal with respect to R(g a ).

In particular, one can show that

R(g 2-δ ) = 1 -δ 2 -δ < R(g 1 ) = 2 -2δ 3 -δ .
This is not surprising. By construction, the largest values of H are more likely to be produced by (2 -δ)X than X, especially if δ is small. In this context, the following lemma explains that the risk function R(g) both depend on the upper tail dependence between g(X) and H and their marginal behaviors.

Lemma 2. Assume that the two limits

c(g) := lim u→∞ P(g(X) > u) P(H > u) ∈ (0, ∞) and χ * (g) := lim u→∞ P(g(X) > u | H > u) ∈ [0, 1]
exist. Then, χ * (g) ≤ c(g) and the limiting risk based on (3) has the expression

R(g) = 1 - χ * (g) 1 + c(g) -χ * (g) . (9)
In particular,

R(g) = 0 if and only if c(g) = χ * (g) = 1.
This lemma indicates that c(g) = 1, i.e. ḡ(X) and H are "asymptotically calibrated", is a necessary condition to have R(ḡ) = 0. In the above example with δ = 0, we have c(g 2 ) = 1 and R(g 2 ) = 1/2 = 0. This means that χ * (g 2 ) < 1. Note that c(g) = 1 implies that the constant χ * (g) simply corresponds to the aforementioned tail dependence coefficient χ, denoted χ(g), and so, the case c(g) = 1 simplifies the expression of the risk

R(g) = 1 - χ(g) 2 -χ(g) .
This equality tells us that any asymptotically calibrated classifier with χ(g) = 0 always produces a risk function R(g) = 1. Consequently, any asymptotically independent classifier is as uninformative as the two naive classifiers. A reasonable strategy will be to dismiss all asymptotically independent classifiers and find/construct new asymptotically dependent classifiers with positive χ(g). But, finding asymptotically dependent classifiers can be complex in practice, and in addition, in some not so exotic setups, this is not always possible. To see this, we consider the simple non-linear regression model in the following lemma.

Lemma 3. Assume that the variable H in (1) is generated by the non-linear regression model

H d = f (X) + N,
where d = represents the equality in distribution, N corresponds to a random noise and X corresponds to the explanatory variables, independent of N . If P(f (X) > u) = o (P(N > u)) , then for any classifier of the type defined by (7), we always have

R(g) = 1.
Hence, no classifier can outperform naive classifiers for this regression model.

Note that even if the forecaster knows exactly the function f (.) and has drawn from the explanatory X, the "ideal" classifier g(x) =f (x) will perform badly, i.e. R(f ) = 1. In addition, the classical trick of using ranks to avoid the problem of marginals discrepancy cannot be applied here. For example, suppose that H is unit Fréchet distributed, then transforming the marginals of X into unit Fréchet random variables, say into X, does not remove the issue as the unobserved noise N has still have heavier tails than f ( X) = f (X) for some function f (.). So, a finer risk measure is needed that is able to distinguish different classifiers in case of asymptotic independence.

Conditional risk and hidden regular variation

The choice of the conditioning set in Lemma 1 brings new possibilities to construct finer risk measures for extremes than R (u) . To do so, we opted for the sets: A 1 = {g(X) > u} and B 1 = {H > u}, (when ε = 1) and A 0 = {g(X) > 0} and B 0 = {H > 0} (when ε = 0). But, changing the size of the sets A ε and B ε to make them closer to A 1 and B 1 , will increase the conditional probabilities

P(B 1 | A 1 ∩ B ε ) and P(A 1 | A ε ∩ B 1 ). A simple choice when ε = 0 or ε = 1 is to set A ε = {g(X) > εu} and B ε = {H > εu} with ε > 0.
This modeling strategy is at the core of hidden regular variation and asymptotic independent models. More precisely, we first need to fix marginal features. We assume that both g(X) and H possess regularly varying tails with indices α g > 0 and α H > 0, respectively. This means that for any ε ∈ (0, 1),

lim u→∞ P(g(X) > u | g(X) > εu) = ε αg and lim u→∞ P(H > u | H > εu) = ε αH .
These limits have to be understood with respect to Equation (6), i.e. the terms P(A 1 | A ε ) and P(B 1 | B ε ). To apply (6), the mixing condition (5) needs to be satisfied. To do so, we opt for an extended version of the framework of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF], i.e.

P[g(X) > u, H > v] = L(u, v)(u -αg v -αH ) 1/2η , ( 10
)
where η ∈ (0, 1] indicates the rate of decay of the joint survival function and L(•, •) is bivariate slowly varying function, i.e. there exists a limit function :

(0, ∞) × (0, ∞) → (0, ∞) defined as (s, t) = lim u→∞ L(us, ut) L(u, u) , s, t > 0
and satisfying (cs, ct) = (s, t) for all c, s, t > 0. The parameter η measures the dependence strength. The case η = 1 corresponds to the asymptotic dependence while η < 1 to the asymptotic independence case. In particular, if η = .5, then independence appears in the extremes. If .5 < η < 1 (0 < η < .5) the extremal dependence is said to be positively (negatively) associated. Now, noticing that (10) corresponds to the mixing condition (5), we can apply Lemma 1, see Appendix A for a proof. Proposition 4. Under the Ramos and Ledford model defined by (10), the following risk function

R ε (g) := lim u→∞ P(g(X; u) = Y (u) | Y (εu) = g(X; εu) = 1) P(Y (u) = 1 or g(X; u) = 1 | Y (εu) = g(X; εu) = 1) , (11)
which will henceforth also be called conditional risk, can be expressed as

R ε (g) = 1 - 1 (ε, 1)ε -αg/2η + (1, ε)ε -αH /2η -1 , for any ε ∈ [0, 1).
Note that η ∈ (0, 1] takes a similar role as χ in the case of the unconditional risk R.

For fixed ε ∈ [0, 1), the risk function R ε (g) decreases with increasing η. So, given all parameters are fixed but η, the forecaster should aim at maximizing η. In practice, two forecasters, say g 1 and g 2 , may produce different (., .) and α g . Consequently, the minimization of R ε (g) can also depend, besides η, on other parameters.

Risk function inference

Concerning the estimation of R ε (g) defined by ( 11), the empirical estimator can be easily computed from the sample (X i , H i ) i=1,...,n . The following proposition describes the asymptotic property of such an estimator.

Proposition 5. Assume that the risk function R ε (g) defined by (11) exists for a sequence of

u n → ∞ such that np g,ε (u n ) → ∞ with p g,ε (u n ) := P(max{g(X; u n ), Y (un) } = 1, H > εu n , g(X; εu n ) = 1). If lim n→∞ np g,ε (u n ) P(g(X; u n ) = Y (un) , H > εu n , g(X; εu n ) = 1) p g,ε (u n ) -R ε (g) = 0,
then the empirical estimator based on a sample (g(X i ; u), {H i > u}) i=1,...,n and defined by

R n,ε (g) = n i=1 1{g(X i ; u n ) = Y (un) i , H i > εu n , g(X i ; εu n ) = 1} n i=1 1{max{g(X i ; u n ), Y (un) i } = 1, H i > εu n , g(X i ; εu n ) = 1} . ( 12
)
converges in distribution in the following way

np g,ε (u n ) R n,ε (g) -R ε (g) n→∞ -→ N (0, R ε (g)(1 -R ε (g)))
.

Simulations

A simple linear setup

Our main simulated example is based on a simple linear regression model but with the feature that the explanatory variables X do not have the same tail behavior and the noise is regularly varying, see Lemma 3. More precisely, the multivariate vector X is defined as follows

         X 1 ∼ Pareto(3) X 2 ∼ Pareto(2) X 3 ∼ Exp(1) X 4 ∼ Exp(2) , (13)
where all X i are independent with X 1 and X 2 Pareto distributed with respective tail index 2 and 3, and X 3 and X 4 exponentially distributed with respective scale parameters 1 and 2. The variable of interest H is simply a linear transform of X 1 tainted by an additive noise

H d = X 1 + N, (14)
where N ∼ Pareto(2) represents an independent noise with heavier tail than X 1 . So, given a sample ({X j,i } 1≤j≤4 , H i ) i=1,...,n with n = 10000, our goal is to compare different classifiers in terms of predicting extreme occurrences, here defined as {H > u} with u equal to the 97th percentile of H. In this simulation setup, it is clear from ( 14) that

Method Main features

Linear classifier Simple binary classifier, parameters estimation based on the minimization of the risk function over the set of contributing variables, theoretical value of R(g θ ) inferred from spectral decomposition.

Logistic regression

Parametric linear model with a lasso penalty, coefficients of less contributing variables are set to zero.

Decision trees

Easy to interpret, gives relative importance of each variables, learns simple decision rules inferred from the input.

Random forests

Builds multiple decision trees combined by majority vote, better predictive power than decision trees.

Support vector machines Finds the best hyperplane to separate two overlapping classes, generally performs better than the other classifiers.

all variables but X 1 are useless to explain H. In addition, Lemma 3 tells us that the relevant information contained in the variable X 1 is hidden by the heavier noise N , i.e. we are in the case of asymptotic independence.

An example of such simulation is given in Figure 1. The left panel displays a scatter plot between H (left axis) and g(X) = X 1 (right axis). As expected, no sign of asymptotic dependence can be found in the upper corner. In the right panel, we remove the mass along the axis (gray points) by conditioning on the joint event A ε ∩ B ε = {g(X) > εu} ∩ {H > εu} with ε = .7, see all dark points. The right panel zooms on these black points and highlights a clear dependence between H and X 1 that was hidden by the heavier noise N in H = X 1 + N .

In practice, we do not know the optimal choice for g(.) and we need to introduce different classifiers and compare them.

Classifiers descriptions

Table 1 below provides the list of classifiers that we compare with our metric (12). This list contains some of the most standard classifiers found in the literature (see, e.g. [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]: logistic regression (Logistic), decision tree (Tree), random forest (RF ) and support vector machine (SVM ).

Except for the linear classifier, we apply them with their built-in cost function that is not necessarily fine-tuned to forecast extremes. This is not an issue because our main goal is to compare existing forecasters, and not to create new ones (see, e.g. [START_REF] Jalalzai | On Binary Classification in Extreme Regions[END_REF], for such developments). Still, to fix a baseline in terms of performance, the linear classifier defined as

g θ (X; u) = +1, θ X > u, -1, θ X ≤ u, θ ∈ [0, ∞) d .
should be optimal for the linear model ( 14), especially if the regression parameters are estimated by minimizing our cost function ( 11). In such a context, we expect the linear classifier to be the best. In Appendix B, Proposition 6 provides the condition of the consistency of the estimator θn,un based on minimizing Rn (g θ (•; u n )) under a regularly varying framework. The binary outputs from the decision tree classifier are explained in Figure 2. The light blue and light green regions represent the set of points that are well predicted by the classifier. On the contrary, wrongly classified points belong to the light yellow and light red regions: either an extreme is predicted when it is not (light red region), or an extreme event is missed (light yellow region). The difference between the left and right panels corresponds to the training set based on either {H > u} or on {H > εu}, i.e. mass removed in the latter case, see also (12).

Implementation and results

We split our simulated data set in two: 70% for a training part, over which we train our different classifiers to get good predictive power; 30% for a testing part, which we use to estimate the risks R 0 (g) and R ε (g). Note that each algorithm has the same inputs, in particular the same binary sequence describing the events {H > u} with u set to be equal to the 97th percentile of H. This cross-validation procedure has been repeated 50 times. The sample used to compute our risk function is based on the bivariate binary vector (g({X j,i } 1≤j≤4 ; u), {H i > u}) i=1,...,n where the output of the classifier g is binary.

In addition, the binary outputs of the classifier g are obtained under the threshold u and the threshold εu, so the training part has to be performed twice (once for each threshold). Then, the empirical risk estimator defined by ( 12) can be computed. Figure 3 shows the sensitivity of the classifier ranking with respect to the value of ε.

As expected from Lemma 3, the top-left panel, that corresponds to the case ε = 0, clearly indicates that our five classifiers cannot outperform naive classifiers as all classifiers have a risk near to one, the worst possible value. To start discriminating classifiers, we need to remove the masses along the axes by setting a value for ε. As we increase ε, the size of sets needed to compute (12) becomes smaller, see the values of n ε in the legend of each panel. Hence, the blue box plots becomes wider as ε increases: a classical bias-variance trade-off. As we know from ( 14) that the true generative process is linear, ε = 0.4 appears as a reasonable value to balance the bias-variance trade-off. More importantly, the overall ranking is not sensitive to the values of ε > 0. In all cases, our linear classifier tailored to handle linear asymptotic independence cases outperforms all the other classifiers. Among the other classifiers, decision tree appears to be the best, but it is still far from the optimal linear solution. Other simulations concerning the regular variation case are available upon request.

Danube river discharges

We now apply our assessment approach to summer daily river discharges (measured in m 3 /s) at 31 stations spread over the upper Danube basin, see Figure 4, and recorded over the time period 1960-2010 in June, July and August. These observations have been studied by the EVT community (see, e.g. [START_REF] Asadi | Extremes on river networks[END_REF][START_REF] Mhalla | Causal mechanism of extreme river discharges in the upper danube basin network[END_REF][START_REF] Gnecco | Causal discovery in heavy-tailed models[END_REF]. This dataset was made available by the Bavarian Environmental Agency (http://www.gkd.bayern.de). To remove temporal clustering in extreme river discharges, [START_REF] Mhalla | Causal mechanism of extreme river discharges in the upper danube basin network[END_REF] in their Section 5 implemented a declustering step. Each station then contains n = 428 observations that we will consider temporally independent. In order to reduce the large discrepancy in terms of discharges magnitude among stations, we force the starting value of all 31 time series to equal zero by subtracting to each station its minimum. Then, we re-normalize each time series by its range (i.e. the difference between the maximum and the minimum of each time series).

These post processing treatments are useful to display and interpret the data at hand and do not impact the classifiers performance.

Although all 31 station recordings are available, we can artificially remove one station and try to predict its values from a given subset of other weather stations. In this section, we remove station 1 (downstream) and try to predict its values from only stations 23 and 24, which are indirect tributaries to the main river flow. So, this setup is complex § for two reasons. First, station 1, as a downstream point that accumulates all discharges, has a much heavier tail than the two tributaries. Second, it is difficult to determine if we are in the asymptotic dependent or independent case, see Figure 5 that displays the scatter plot between the hidden station (y-axis with station 1) and the two tributaries (x-axis, stations 23 and 24). In this graph, the threshold u is taken to be equal to the 85th quantile of X 1 and we choose ε = 0.6. Figure 6 summarizes our findings. Removing the mass on the axes when thresholding by εu implies that only 190 points remain from §Section C treats a simpler case where station 1 is predicted from the whole set of remaining stations. In this case, strong dependencies among station 1 and other stations can be observed. So, the main issue is to select these stations, a problem discussed in Section B. the original length of 428 data points per station. This can explain why, looking at Figure 6, the uncertainty in the risk estimate increases when considering R ε (.) (blue boxes) instead of R 0 (.) (red boxes).

Unlike the simulation example in Section 3.1, it is not clear to assess whether our river discharges analysis of our three selected weather stations belongs to the framework of asymptotic independence or not. Still, it is reassuring that the ranking of the classifiers in Figure 6 appears to be insensitive to the values of R 0 or R ε . i.e. whether the data are asymptotically dependent or not. This hints that, among all the classifiers, the logistic regression with lasso penalty seems to perform better than the four other classification methods. This ranking of classifiers is specific to this particular example. No general conclusions about lasso techniques for extremes should be drawn.

Besides this river example, we advocate practitioners to compute risk functions that can both handle the asymptotic dependent and independence cases. This also complements the recent tools used to discriminate between the two cases (see, e.g. [START_REF] Ahmed | Recognizing a spatial extreme dependence structure: A deep learning approach[END_REF]. In addition, the linear classifier could provide a simple benchmark with well understood properties with respect to R 0 , see Proposition 6.

Supplementary Materials

A R package is available on GitHub that implements the empirical estimation of the risk function developed in this paper (https://github.com/jlegrand35/ExtremesBinary Classifier) and can be used either to reproduce the results of the conducted classifier comparisons or to perform new comparisons using other binary classifiers. The data used in the application are available in the R package graphicalExtremes [START_REF] Engelke | graphicalExtremes: Statistical Methodology for Graphical Extreme Value Models[END_REF]. Fig. 5. Summer daily measurements of river discharges from station 1 (y-axis) against station 23 (left x-axis) and station 24 (right x-axis) of Figure 4. The blue dotted lines are the threshold value u (the 85th percentile of station 1). The black dots on the graphs in the top row are the values such that min(X j , X 1 ) > εu where εu = 0.6u. The graphs in the bottom row corresponds to the same data plots but on unit Fr échet scale highlighting potential asymptotic independence in the data. Fig. 6. Estimation of R ε (g) for five different classifiers (cross-validation with 50 repetitions, 70% train and 30% test), threshold is the 0.85 quantile of H. In red (left) are the estimates when ε = 0 and in blue (right) for ε = 0.6. The length of the testing set is equal to 129, this leads to around 60 points such that H > εu and nearly 20 points such that H > u.

A. Proofs

A.1. Proof of Lemma 1:

We can write that

P(A 1 ∪ B 1 | A ε ∩ B ε ) = P((A 1 ∪ B 1 ) ∩ (A ε ∩ B ε )) P(A ε ∩ B ε ) , = P(A 1 ∩ B ε ) + P(A ε ∩ B 1 ) -P(A 1 ∩ B 1 ) P(A ε ∩ B ε )
.

In the same way, we have

P(A 1 B 1 |A ε ∩ B ε ) = P(A 1 ∩ B ε ) + P(A ε ∩ B 1 ) -2 P(A 1 ∩ B 1 ) P(A ε ∩ B ε )
.

Hence, we deduce that

R(A ε , B ε ) = P(A 1 B 1 |A ε ∩ B ε ) P(A 1 ∪ B 1 |A ε ∩ B ε ) = 1 - P(A 1 ∩ B ε ) P(A 1 ∩ B 1 ) + P(A ε ∩ B 1 ) P(A 1 ∩ B 1 ) -1 -1
.

The expression given by (4) follows.

Item (a) of the lemma is based on the following inequality

P(U | V ) ≥ P(U | W ), if the sets U , V and W satisfy U ⊂ V ⊂ W .
For item (b), note that

R(A ε , B ε ) ≤ R(A ε , B ε ) ⇐⇒ (1 -R(A ε , B ε )) -1 ≤ 1 -R(A ε , B ε ) -1 ⇐⇒ P(A 1 ∩ B ε ) P(A 1 ∩ B 1 ) + P(A ε ∩ B 1 ) P(A 1 ∩ B 1 ) ≤ P(A 1 ∩ B ε ) P(A 1 ∩ B 1 ) + P(A ε ∩ B 1 ) P(A 1 ∩ B 1 ) , ⇐⇒ P(A 1 ∩ B ε ) + P(A ε ∩ B 1 ) ≤ P(A 1 ∩ B 1 ) P(A 1 ∩ B 1 ) [P(A 1 ∩ B ε ) + P(A ε ∩ B 1 )]. As P(A 1 ∩ B 1 ) = P(A 1 ∩ B 1 ) and P(A 1 ∩ B ε ) = P(A 1 ∩ B ε ) and P(A ε ∩ B 1 ) ≤ P(A ε ∩ B 1 ), then R(A ε , B ε ) ≤ R(A ε , B ε ).
This provides the second statement (b) since we assume that

A 1 = A 1 and A ε ⊆ A ε .
Item (c) is a direct consequence of (4).

A.2. Proof of Lemma 2

First, we note that, by definition,

χ * (g) = c(g) • lim u→∞ P(H > u | g(X) > u) ≤ c(g).
By taking the limits separately for every single term in (8), we directly obtain (9). From that, we can immediately see that

R(g) = 0 if and only if 1 + c(g) = 2χ * (g).
Now, assume that χ * (g) < c(g). Then, R(g) = 0 would imply that 1 + χ * (g) < 2χ * (g) which is equivalent to χ * (g) > 1 in contradiction to the definition of χ * (g). Thus, R(g) = 0 implies χ * (g) = c(g). Plugging this into (9), we obtain c(g) = χ * (g) = 1. Conversely, c(g) = χ * (g) = 1 obviously leads to R(g) = 0.

A.3. Proof of Lemma 3: From Eq. ( 8), we know that we can only get R(g) < 1, only if H and g(X) are tail equivalent. Thus, this will be assumed in the following. Let w such that δ := lim u→∞ w(u)/u ∈ (0, 1), then for any positive u we can write that P(f (X) + N > u, g(X) > u) P(g(X) > u) ≤ P(f (X) > w(u), g(X) > u) P(g(X) > u)

+ P(N > u -w(u), g(X) > u) P(g(X) > u)

.

Since g(X) and N are independent, the second term reduces to P(N > u -w(u)) which converges to 0 since u -w(u) ∼ (1 -δ)u as u gets large.

For the first term, we rewrite the ratio as follows

P(f (X) > w(u), g(X) > u) P(g(X) > u) ≤ P(f (X) > w(u)) P(g(X) > u) ≤ P(f (X) > w(u)) P(N > w(u)) P(N > w(u)) P(N > u) P(N > u) P(g(X) > u)
.

Since w(u) → ∞ and P(f (X) > u) = o (P(N > u)), the ratio P(f (X)>w(u)) P(N >w(u)) goes to 0 as u gets large. From the assumption w(u) ∼ δu, P(N >w(u)) P(N >u) behaves as a constant when u → ∞.

The only remaining term is P(N >u) P(g(X)>u) which converges to a constant due to tail equivalence. So, lim P(g(X) > u | H > u) = 0.

A.4. Proof of Proposition 4:

In Lemma 1 we fix A ε = {g(X) > εu} and B ε = {H > εu} and A 1 = {g(X) > u} and B 1 = {H > u}.

P(H > u or g(X) > u | min{H, g(X)} > εu) = P(A 1 ∪ B 1 | A ε ∩ B ε ), and P(g(X; u) = Y (u) | min{H, g(X)} > εu) = P(A 1 B 1 |A ε ∩ B ε ).
The Ramos and Ledford model corresponds to the special of case of ( 5)

P[g(X) > εu, H > ε u] = L(εu, ε u)((εu) -αg (ε u) -αH ) 1/2η = P(A ε ∩B ε ) = c ε,ε P a (A ε ) P b (B ε ) with P(A ε ) = P(g(X) > εu) = L g (εu)ε -αg u -αg , P(B ε ) = P(H > εu) = L H (εu)ε -αH u -αH and c ε,ε = L(εu, ε u) (L g (εu)L H (ε u)) 1/(2η) , a = b = 1 2η .
Then, from ( 6)

R(A ε , B ε ) = 1 - c ε,1 c 1,1 (P(A 1 |A ε )) -a + c 1,ε c 1,1 (P(B 1 |B ε )) -b -1 -1
.

Letting u gets large provides the required result. .

A.5. Proof of Proposition 5:

As we assume that R ε (g) exists (for some ε ≥ 0), for u n → ∞ such that np g,ε (u n ) → ∞, we obtain that

np g,ε (u n )   n i=1 1{g(Xi,un) =Y (un) i
,Hi>εun,g(Xi;εun)=1} npg,ε(un) n i=1 1{max{g(Xi,un),Y (un ) i }=1,Hi>εun,g(Xi;εun)=1} npg,ε(un)   -P(g(X,un) =Y (un ) ,H>εun,g(X;εun)=1) pg,ε(un)

1 n→∞ -→ N 0 0 , R ε (g) R ε (g) R ε (g) 1
Provided that the bias is negligible, i.e.

lim n→∞ np g,ε (u n ) P(g(X, u n ) = Y (un) , H > εu n , g(X; εu n ) = 1) p g,ε (u n ) -R ε (g) = 0,
the Delta method yields

np g,ε (u n ) R n,ε (g) -R ε (g) n→∞ -→ N (0, R ε (g)(1 -R ε (g)))
.

B. Linear Classifiers

B.1. Definition, Basic Properties and Inference

In this section, we consider a specific type of classifiers which in this paper is referred to as linear classifiers, i.e. classifiers of the form

g θ (X; u) = +1, θ X > u, -1, θ X ≤ u, θ ∈ [0, ∞) d .
To obtain an optimal linear classifier of g θ (X; u), i.e. some weight vector θ * such that the classification risk R(g θ * ) gets minimal, we need to impose some joint extremal dependence structure on X and H from (2). Even though some of the results can also be obtained in a similar manner in a more general framework for the conditional risk R ε , henceforth, we will focus on the asymptotically dependent case where we might find some optimal classifier with unconditional risk R(g θ * ) < 1. As discussed before, in this case, at least one component of X needs to have a similar tail behavior as H. A natural assumption is therefore that (X, H) is jointly regularly varying on [0, ∞) with index α > 0, i.e. there exist an α-Pareto random variable P and, independently of P , a random vector (Γ, Ω) ∈ [0, ∞) d × [0, ∞), the so-called spectral tail vector, on the unit sphere {x ∈ [0, ∞) d+1 :

x ∞ = 1} such that L (X, H) ∞ u , (X, H) (X, H) ∞ (X, H) ∞ > u → w L(P, (Γ, Ω)), as u → ∞.
Here, we additionally assume that P( Γ ∞ > 0) > 0 and P(Ω > 0) > 0. In this setup, X ∞ and H are tail equivalent in the sense that

lim u→∞ P( X ∞ > u) P(H > u) = lim u→∞ P( X ∞ > u | (X, H) ∞ > u) P(H > u | (X, H) ∞ > u) = P(P • Γ ∞ > 1) P(P • Ω > 1) = E( Γ α ∞ ) E(Ω α ) ∈ (0, ∞),
This is the minimal requirement on the link between the covariates X and the unobserved extremes of H essentially saying that at least one component of X is tail-equivalent to H. It is important to highlight that we do not exclude the case that Γ i = 0 a.s. for some i ∈ {1, . . . , d} which means that X i possesses a lighter tail than H. This property can be read off from the quantity

c i = lim u→∞ P(X i > u) P(H > u) = E(Γ α i ) E(Ω α )
. Thus, Γ i = 0 a.s. if and only if c i = 0. By including this case, we therefore admit that most of the components of X may not contribute to the extremes of the vector H.

Under these conditions, we obtain that, for all θ ∈ [0, ∞) d , the classifier g θ is an extremal classifier as

R(g θ ) = lim u→∞ P(H > u or g(X) > u) -1 P[max{θ X, H} > u] -P[min{θ X, H} > u] = E max{θ Γ, Ω} α -E min{θ Γ, Ω} α E (max{θ Γ, Ω} α ) = 1 - E min{θ Γ, Ω} α E (max{θ Γ, Ω} α ) ∈ [0, 1]. (15)
Equation ( 15) implies that the function θ → R(g θ ) is well-defined and continuous on [0, ∞) d . Its value does not depend on those components θ i for which Γ i = 0 a.s., which is equivalent to c i = 0 as discussed above. Thus, in the following, we will consider this function only on the parameter set

C = {θ ∈ [0, ∞) d : θ i = 0 for all i s.t. c i = 0},
containing all the relevant information -here, note that, in practice, identifying the components i ∈ {1, . . . , d} such that c i = 0 a.s., from a given data set is a necessary step for the correct specification of the set C.

From the consideration in the introduction, it can be easily seen that R(g 0 ) = 1 -the case θ = 0 corresponds to the trivial always optimistic classifier. Furthermore, denoting the set of indices j with c j > 0 by J, we can see that

R(g θ ) ≥ 1 - E(Ω α ) θ α ∞ • min j∈J E(Γ α j ) → 1 (16)
as θ ∞ → ∞. By the continuity of θ → R(g θ ), we obtain that the function attains a global minimum on the domain C.

Proposition 6. Additionally to the assumptions above on the joint distribution of (X, H) with α > 1, assume that there exists a function a(u) with a(u) → 0 as u → ∞ such that P

(u -1 X ∈ A | (X, H) > u) ≤ (1 + a(u)) P(P Γ ∈ A) (17)
for all A ⊂ [0, ∞). Furthermore, let u n → ∞ and n P(H > u n ) → ∞ such that, for every compact subset K ⊂ C,

sup θ∈K n P(H > u n ) P(g θ (X; u n ) = Y un ) P(H > u n ) - E max{θ Γ, Ω} α -E min{θ Γ, Ω} α E(Ω α ) = 0 (18)
and

sup θ∈K n P(H > u n ) P(max{g θ (X; u n ), Y un } = 1) P(H > u n ) - E max{θ Γ, Ω} α E(Ω α ) = 0. ( 19
)
If the function θ → R(g θ ) has a unique minimizer θ * in C, then the estimator

θn,un = argmin θ∈C Rn (g θ (•; u n )).
is consistent, i.e. θn,un → p θ * .

Given the set C, this result provides a strategy to find the optimal θ, i.e., the best linear classifier. Determining the set C requires the identification of the relevant features, i.e. the index set J such that c j > 0 if and only if j ∈ J. This is discussed in more detail in the following subsection.

B.2. Feature Selection

The notion of sparsity quickly comes into play when doing classification. This is all the more true when one is only interested in the extremes. Among the whole data set, only a small proportion will truly contribute to the extremal behavior of the variable of interest. Here, we develop a method to identify the informative signals in terms of extremes among a large data set, assuming that (X, H) is jointly regularly varying. For a comprehensive review of existing methods on sparsity and multivariate extremes we highly recommend the work of [START_REF] Engelke | Sparse structures for multivariate extremes[END_REF].

As we have seen above, for linear classifiers, all the relevant features X i necessarily satisfy c i > 0. Thus, feature selection can be based on estimation of the c i which can be done according to the following proposition. exists, then, we have

n P(H > u n ) n j=1 1{X j,i > u n } n j=1 1{H j > u n } -c i → N (0, c i • [1 -2χ * i + c i ]) .

C. River network

This section deals with a simpler case than the application in Section 4. Here an application could be the following: we want to know which stations should continue to be maintained to prevent extreme floods and maybe some stations are not necessary.

Here we still want to predict the extreme events at station 1. We define an extreme event as an event exceeding the 85th quantile of X 1 and we assume that the whole set of remaining stations is available. In this case, strong dependencies among station 1 and other stations can be observed. So, the main issue is to select these stations following the procedure presented in Appendix B.

As in the simulated study, we identify the stations that may not contribute to the extremes of X 1 by estimating the c i . The estimation of the set C on all the data is presented in Table 2. We found that among the 30 stations, only three stations are relevant: stations 2, 13 and 30. Looking at Figure 4, these stations correspond to the stations closest to X 1 . Figure 7 shows the scatter plots between these stations and station 1.

Once the contributing variables have been identified, we compare the performance of several classifiers, on the one hand keeping all the data and on the other hand keeping If the function θ → R(g θ ) has a unique minimizer θ * , then, necessarily, R(g θ * ) < 1. Now, similarly to the notation above, let J denote the set of indices j with c j > 0, and let us consider θ ∈ C such that θ ∞ > k 0 for some constant k 0 > 0. Then,

R n (g θ ) = 1 - n i=1 1{min(θ X i , H i } > u n )} n i=1 1{max(θ X i , H i } > u n )} ≥ 1 - n i=1 1{H i > u n } min j∈J n i=1 1{k 0 X ij > u n } n→∞ -→ p 1 -max j∈J E(Ω α ) k α 0 E(Γ α j )
where the right-hand side goes to 1 as k 0 → ∞. Thus, as R n (g θ * ) → p R(g θ * ) < 1, we obtain that, for sufficiently large k 0 θ * , with probability going to one,

R n (g θ * ) ≤ min θ∈C\[0,k0] d R n (g θ )
and, consequently, argmin θ∈C R n (g θ ) = argmin θ∈C∩[0,k0] d R n (g θ ).

Now, we note that, by Lemma 8, the bias conditions ( 18) and ( 19) and the functional delta method, R n (g θ ) converges in probability to R(g θ ) uniformly on every compact subset of C. In particular, Taillardat, M., Fougères, A.-L., Naveau, P., and de Fondeville, R. (2019). Extreme events evaluation using CRPS distributions. arXiv:1905.04022. Van der Vaart, A. W. and Wellner, J. A. (1996). Weak Convergence and Empirical Processes. Springer Series in Statistics. Springer, New York.

Fig. 1 .

 1 Fig. 1. Simulated example for Lemma 3, for more details see equations (13) and (14).

Fig. 2 .

 2 Fig. 2. Example of predicted binary output for the decision tree classifier versus the true values of H (on a logarithmic scale). The classifier has been trained twice as explained below. On the left-hand side the classifier is trained with the subset {H > u} and the horizontal blue dashed line represents the threshold value u. On the right-hand side the classifier is trained with the subset {H > εu} and the horizontal red dashed line represents the threshold value εu = 0.4 × u. The results shown are the predicted values made upon the testing set whose sample size is equal to 3000 (30% of the data).

Fig. 3 .

 3 Fig. 3. Estimation of R ε (g) for different classifiers (cross-validation with 50 repetitions). In red (top-left) are the estimates when ε = 0 and in blue for different values of ε > 0 (ε ∈ {0.4, 0.6, 0.8}). At the top of each plot, the value of ε and the number of points such that H > εu from the testing set are given.

Fig. 4 .

 4 Fig. 4. River map of the upper Danube basin, showing sites of the 31 gauging stations along the Danube and its tributaries. Water flows toward gauging station 1. The stations represented by a green triangle shaped dot are the three stations of interest as described in Section 4.

Proposition 7 .

 7 Assume thatc i = lim u→∞ P(X i > u) P(H > u) exists. If u n → ∞ and n P(H > u n ) → ∞, then n j=1 1{X j,i > u n } n j=1 1{H j > u n } → p c i .

Fig. 7 .

 7 Fig. 7.Correlation plot between X 1 and the stations considered to contribute to the extremes of X 1 according to the table 2 (i.e. for which c i = 0). The blue dotted lines represent the threshold u defined by the 85th percentile of X 1 .

Fig. 8 .

 8 Fig. 8. Estimation of R 0 (g) for different classifiers g (cross-validation -70% train, 30% test -with 50 repetitions). The red distributions come from the estimation with all the stations, and the green distributions represent the estimations with only the variables having ĉi = 0.

  sup θ∈C∩[0,K] d R n (g θ ) -R(g θ ) → p 0. Thus, argmin θ∈C∩[0,K] d R n (g θ ) → p argmin θ∈C∩[0,K] d R(g θ ) = θ * .

Table 1 .

 1 Summary and key features of the different classifiers studied. See for example[START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF] for a comprehensive review of the last four classification methods.

Acknowledgment

Part of this work was supported by the DAMOCLES-COST-ACTION on compound events, the French national program (FRAISE-LEFE/INSU and 80 PRIME CNRS-INSU), and the European H2020 XAIDA (Grant agreement ID: 101003469) and funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2075 -390740016. The authors also acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-20-CE40-0025-01 (T-REX project), the ANR-Melody and the Stuttgart Center for Simulation Science (SimTech).

0.00 X10 0.00 X11 0.00 X12 0.00 X13 0.30 X14 0.00 X15 0.00 X16 0.00 c i X17 0.00 X18 0.00 X19 0.00 X20 0.00 X21 0.00 X22 0.00 X23 0.00 X24 0.00 X25 0.00 X26 0.00 X27 0.00 X28 0.00 X29 0.00 X30 0.02 X31 0.00 only the informative stations. Since there is a strong dependence between the data, we assume that it is sufficient here to look at the risk R 0 . Comparison results are shown in Figure 8.

By definition of the linear classifier, the estimation is already done by keeping only the informative variables, which is why the estimates are identical for this specific classifier. As for the other classifiers, we see some improvements when keeping only the informative variables: the risk estimates are slightly smaller. This means that even if we remove a lot of information by going from 30 explanatory variables to 3, these 3 remaining stations contain all the information in terms of extremes of station 1.

D. Proofs of the appendix

D.1. Proof of Proposition 6

The proof is based on the following lemma which is proven in Subsection D.2. Lemma 8. Under the assumptions from Proposition 6, for every compact subset K ⊂ C, the sequences of processes {A n (θ), θ ∈ K} and {B n (θ), θ ∈ K} defined by

converge to centered Gaussian processes {A(θ), θ ∈ K} and {B(θ), θ ∈ K}, respectively, weakly in ∞ (K).

D.2. Proof of Lemma 8

We will proof the lemma by applying the Central Limit Theorem 2.11.9 in Van der Vaart and Wellner (1996). To this end, we define the function spaces A = {a θ , θ ∈ K} and B = {b θ , θ ∈ K} where

Then, with

for l = 1, . . . , n, we have that

for all l = 1, . . . , n and n ∈ N. Consequently, we check the Lindeberg condition: For k ∈ N, we have

as n P(H > u n ) → ∞ by definition. For k = 2, we obtain a Lindeberg type condition that ensures convergence of A n and B n to A and B, respectively, in terms of finite-dimensional distributions. For k = 1, we obtain the Lindeberg type condition of Theorem 2.11.9 in Van der Vaart and Wellner (1996).

It remains to check the equi-continuity condition. In the following, to simplify notation, we assume that

, 1} and the probability that any of those two expressions is equal to one is bounded by the probability

where we use that b X > u n implies that X > u n /( K d) with K = sup x∈K x ∞ . Making use of the fact that P ( (X, H) > u n /( K d)) ≤ C 0 ( K d) α P(H > u n ) for some constant C 0 > 0 and the bound given by Equation ( 17), we obtain that

which tends to 0 as δ → 0. From this inequality, it can also be seen that any partition of K into hypercubes with length ε 2 /(2C 0 ) leads to a valid ε-bracketing, i.e. the number N ε ∝ 1/ε 2d grows with a power rate and, so, log(N ε ) is integrable. Thus, by Theorem 2.11.9, the processes A n and B n converge to Gaussian processes A and B, weakly in ∞ (K).


Jerome Dubail
email: jerome.dubail@cea.fr

Jesper Lykke Jacobsen
email: jesper.jacobsen@ens.fr

Hubert Saleur
email: hubert.saleur@cea.fr

Jérôme Dubail

Bulk and boundary critical behaviour of thin and thick domain walls in the two-dimensional Potts model

Keywords: conformal field theory, correlation functions, loop models and polymers, solvable lattice models

come

Bulk and boundary critical behaviour of thin and thick domain walls in the two-dimensional Potts model

Introduction

Many of the key developments in the study of two-dimensional critical phenomena originate from the study of just a very few lattice models [START_REF] Wu | Exactly Solved Models: A Journey in Statistical Mechanics[END_REF][START_REF] Baxter | Exactly Solved Models in Statistical Mechanics[END_REF]. These include dimer coverings, the Ising and six-vertex models, the O(n) model, and the Q-state Potts model. The latter two models are particularly important, as they can be formulated in terms of geometrical degrees of freedom that in turn describe extended fluctuating objects, such as domain walls in magnets, percolation clusters, and polymers adsorbed on √ walls and interfaces [START_REF] Eisenriegler | Polymers Near Surfaces[END_REF]. These objects are in fact conformally invariant whenever n or Q lies in the range 2 < n 2, and their fluctuations are characterized by critical exponents.

The purpose of this paper is to revisit the geometrical formulation of the twodimensional Potts model. In particular, we shall define a set of geometrical observables which are different from those considered in most of the existing literature. These observables are simply the domain walls in the formulation of the Potts model in terms of Q-component spins. Not only are they in many ways more natural and physically more relevant than the Fortuin-Kasteleyn (FK) clusters [START_REF] Fortuin | [END_REF] usually considered, but they also

Σ - → - ≤
turn out to have a richer critical behaviour and a more complete set of critical exponents. These features stem in part from the fact that there are actually two different kinds of spin domain walls (thin and thick), as we shall explain shortly.

The Q-state Potts model is defined by the partition function

Z = σ (ij)∈E exp Kδσ i,σj , ( 1
)
where K is the coupling between spins σi = 1, 2 , . . . , Q along the edges E of some lattice L. For simplicity we shall take L to be the square lattice in the computations below, whereas we use the triangular lattice in the figures. The Kronecker delta function δσ i,σj equals 1 if σi = σj, and 0 otherwise. The usual route [START_REF] Fortuin | [END_REF] is to write the obvious identity exp Kδσ i,σj = 1 + vδσ i,σj , (

with v = e K 1, and to expand Z in powers of v. The result is an expression of Z as a sum over FK clusters with weight v per unit length and fugacity Q per connected component.

Alternatively one can think of the FK clusters in terms of their √ surrounding hulls, which are self-avoiding and mutually avoiding loops with fugacity n = Q [START_REF] Baxter | Exactly Solved Models in Statistical Mechanics[END_REF]. Most features of these FK clusters and loops, in the critical regime 2 < n 2, are by now under complete control, thanks to the combined powers of conformal field theory (CFT) and Schramm-Loewner evolution (SLE) [5,6]. In particular, the loops behave like the SLE trace in the continuum limit [6], and viewing them as contour lines of a (deformed) Gaussian free field leads to the Coulomb gas (CG) approach to CFT [5,[START_REF] Nienhuis | Loop Models[END_REF]. Our understanding of the critical properties of FK clusters and loops can be considered almost complete, although some of their more intricate observables-such as the backbone [8] and shortest-path [9] dimensions-are still unknown.

The expansion in powers of v has some pleasing features, notably that the FK clusters coincide with percolation clusters in the formal limit Q 1. But apart from that it is somewhat artificial in view of the original formulation [START_REF] Wu | Exactly Solved Models: A Journey in Statistical Mechanics[END_REF] in terms of Q-component spins. In particular, although two lattice sites belonging to the same FK cluster will necessarily have the same spin value, the converse is not true. It would therefore seem more natural to perform the expansion in powers of e K and consider as basic geometrical observables the connected domains with a constant value of the Potts spin, henceforth referred to as spin clusters.

The properties of these spin clusters have as a rule remained ill understood. This is particularly frustrating, since those are the very clusters that one would observe in an actual experiment on a magnetic alloy in the Potts universality class. The case of the Ising model Q = 2 is an exception to this rule [10,11], but this is due to its 'coincidental' equivalence to the O(n) vector model with n = 1. Indeed, defining the Ising spins on the triangular lattice, the corresponding O(n) model is described by self-avoiding and mutually avoiding loops on the hexagonal lattice, and such loops are readily treated by CFT and SLE techniques.

For general Q, the salient feature of Potts spin clusters is that the domain walls separating different clusters undergo branchings and crossings (see figure 1). These phenomena are however absent for Q = 2 with the above choice of lattice. It is precisely these branchings and crossings that make the application of exact techniques-such as CG mappings or Bethe ansatz diagonalization-very difficult, if not impossible. The belief that spin clusters are indeed conformally invariant for other values of Q = 2 in the critical regime has even been challenged at times, but seems however well established now [12,13].

Some progress has been accomplished in the Q = 3 case [11,14] by speculating that the spin clusters in the critical Potts model would be equivalent to FK clusters in the tricritical Potts model [15,14] having the same central charge. This equivalence was stated in [16] for general (non-integer) Q [0, 4], in the context of a model in which FK and spin clusters were constructed in a unified way by connecting neighbouring vertices with a variable bond-percolation probability p. With the exception of the Q = 2 case [10], the evidence for this equivalence is however mainly numerical and for Q = 3 only [16], and moreover restricted so far to the simplest geometrical questions [17], such as that of the fractal dimension of a bulk cluster. The equivalence can also be understood as a relationship with the dilute O(n) model [18].

The Potts model has been used recently to build a new class of 2D quantum lattice models that exhibit topological order [19]. Both FK clusters and domain walls between spin clusters are important in these models. We should also mention that spin domain walls have been studied numerically-and their difference from the FK clusters pointed out-in the context of the ZN parafermionic models [20].

Apart from issues of branching and crossing, another major hurdle in the study of Potts spin clusters has come from the lack of a formulation that can be conveniently extended to Q a real variable. In the case of loops surrounding the FK clusters, this formulation led naturally to the introduction of powerful algebraic tools via the Temperley-Lieb (TL) algebra, and to the equivalence with the six-vertex model-the eventual key to the exact solution of the problem [5]. Factors of Q then appear naturally through a parameter in the TL algebra, or-via a geometrical construction-as complex vertex weights in the six-vertex formulation. Also for spin clusters, Q can be promoted to an arbitrary variable: the weight of a set of spin clusters is simply the chromatic polynomial of the graph dual to the domain walls. From the point of view of the TL algebra, the domain walls are composite (spin-1) objects, and hence more complicated than the FK clusters and loops. Recent work on the related Birman-Wenzl-Murakami (BWM) algebra [22,21] suggests that this formulation might be amenable to the standard algebraic and Bethe ansatz techniques, although such a lofty goal has not been achieved so far.

In this paper we report major progress towards the understanding of Potts spin clusters, in both the bulk and the boundary case. A brief account on the bulk case has recently appeared elsewhere [23]. Our results are of two kinds. On the one hand, we develop a transfer matrix technique which allows the formulation and numerical study of the spin clusters for all real Q. On the other hand, we identify the geometrical events that give rise to conformal correlations, and provide exact (albeit numerically determined) expressions for infinite families of critical exponents, similar to the familiar 'L-legs' exponents [5,[START_REF] Nienhuis | Loop Models[END_REF] for TL loops. Surprisingly, we find that geometrical properties of spin clusters encompass all integer indices (r, s) in the Kac table hr,s. An analytical derivation of our results appears for now beyond reach, in part because the algebraic properties of our transfer matrix are still ill understood. We do however provide some exact results based on an approach which does not involve a CG mapping, but rather the use of a massless scattering description.

Apart from the bulk critical exponents [23], we report here the boundary critical exponents for free, fixed, and mixed boundary conditions. This means that the Potts spins on a segment of the boundary are restricted to taking Q1 values, where Q1 can assume any real value and is in general different from the number of states Q taken by bulk spins. These results can be viewed as a further step in the programme [24] of classifying non-unitary boundary conditions in 2D geometrical models.

Domain wall expansion

The domain wall (DW) expansion of (1) involves all possible configurations of domain walls that can be drawn on the dual of (see figure 2). A DW configuration is given by a graph G (not necessarily connected). The faces of G are the spin clusters. Since we do

{ } - - ∝ ≤ ≤ × ∝
not specify the colour of each of these clusters, a DW configuration has to be weighted by the chromatic polynomial χ G ˆ(Q) of the dual graph G ˆ. Initially χ G ˆ(Q) is defined as the number of colourings of the vertices of the graph G ˆ, using colours 1, 2, . . . , Q , with the constraint that neighbouring vertices have different colours. This is indeed a polynomial in Q for any G, and so can be evaluated for any real Q (but χ G ˆ(Q) is integer only when Q is integer). For example, the chromatic polynomial of the graph G ˆ on figure 2 is 7 . The partition function (1) can thus be written as a sum over all possible DW configurations

Q(Q 1) 7 (Q 2)
Z = e N K Σ e -K length(G) χ (Q) ( 3
)
G
where N is the number of spins, and length (G) denotes the total length of the domain walls.

The fundamental geometric object that we consider is a connected part of a domain wall that separates two clusters. One can ask how the probability that a certain number of such DW connect a small (in units of the lattice spacing) neighbourhood A to another small neighbourhood B decays when the distance x between A and B increases. Each DW separates two spin clusters which connect A and B. There are in fact two kinds of such DW, depending on the relative colouring of the two clusters that are separated. If the two clusters have different colours, they can touch, so the DW is thin (see figure 3(a)). If the two clusters have the same colour, then they cannot touch (otherwise they would not be distinct), so the DW has to be thick (see figure 3(b)).

Several different geometries are of interest. In the bulk case, the neighbourhoods A and B are at arbitrary, but widely separated, locations in the infinite plane. This is conformally equivalent to an infinitely long cylinder-a strip with periodic boundary conditions-with A and B situated at the two extremities. In the boundary case, the geometry is that of the upper half-plane, with A being at the origin and B far away from the real axis. The boundary conditions are taken to be free, fixed or mixed on the positive real axis, and free on the negative real axis. In other words, the Potts spins can take Q1 (resp. Q) values on the positive (resp. negative) real axis. This geometry is illustrated in figure 4; it is conformally equivalent to an infinitely long strip, with A and B situated at respectively the upper and lower extremities, and spins on the left (resp. right) rim taking Q1 (resp. Q) values. We shall use the cylinder and strip geometries below to conduct our numerical calculations.

We can now state the central claims of this paper. Consider the 2D Potts model for any real Q in the critical regime 0 Q 4. Then the probability P that the two regions A and B, with separation x 1, are connected by l1 thin DW and l2 thick DW decays algebraically. In the bulk case (geometry of the plane) the corresponding critical exponent is denoted as h(Q, l1, l2), and we have P

x -4h(Q,l1 ,l 2 ) . Equivalently, on a long cylinder of size L l with l L, and A and B identified with the opposite ends of the cylinder, the decay is exponential: P e -4π(l/L) h(Q,l1,l2) . In the boundary case (geometry of the halfplane) the critical exponent for free boundary conditions is denoted as h ˜(Q, l1, l2).

Below we check these assertions numerically, and we observe that the numerical values of the exponents match the formulae where we have used the Kac parametrization of CFT (r

J. Stat. Mech. (2010) P12026 G ˆ h(Q, l1, l2) = hl 1-l2 ,2l1 , h ˜(Q, l 1 , l 2 ) = h 1+2(l1-l2),1+4l1 , (4)

≤ ≤ ∈

-sκ/4) 2 -(1 -κ/4) 2 h r,s = (5) κ and 2 κ 4 parametrizes Q = 4(cos(κπ/4)) 2 [0, 4].
It remains to give our results for the mixed boundary conditions with parameter Q1 (for which fixed boundary conditions corresponds to the special case Q1 = 1). These follow by conformal fusion of the free-to-mixed boundary condition changing operator Φr 0 ,s0 that has been worked out in [24] with the operator Φ1+2(l 1 -l2),1+4l1 that inserts the required number of DW with free boundary conditions. Details about this will be given below in section 5.

We close this section with two remarks about the precise interpretation of the main result [START_REF] Fortuin | [END_REF].

First, in the expression (4) for the bulk exponent h(Q, l1, l2), the case (l1, l2) = (0, 1) of a single DW (which must then necessarily be thick, due to the periodic boundary conditions) is special. In that case (4) remains valid only if the DW is forbidden to wrap -≡ - around the neighbourhoods A and B (in the geometry of the plane), or around the periodic direction (in the equivalent cylinder geometry). Without that restriction being imposed, we obtain another result:

h(Q, 0, 1) = h0,1/2 (wrapping allowed). (6)
Second, in the boundary case the numbers l1 and l2 appearing in (4) must be defined more carefully. Indeed, when there are l l1 + l2 propagating spin clusters, only the nature ('thin' or 'thick') or the l 1 DW separating them is clear from the above definition, whereas it is not yet clear how to characterize the 'half-domain walls' that separate the two outermost spin clusters from the boundaries. The correct and unambiguous definition of l1 and l2 is as follows: let the leftmost spin cluster contribute one unit to l1, and let each of the l 1 subsequent spin clusters contribute one unit to l1 (resp. l2) if it has a different colour to (resp. the same colour as) the cluster on its left.

Transfer matrix formulation

The DW expansion (3) may appear unwieldy and difficult to study numerically for noninteger Q. There nevertheless exist several Monte Carlo methods for studying the Potts model when Q is not an integer [25]- [27], all of which are roughly speaking based on the FK cluster representation. It is possible to reconstruct the spin clusters from the Chayes-Machta algorithm [26] by performing the bond-adding step at zero temperature. Following our brief account [23], this method has very recently been used [28] to verify a special case of (4) corresponding to h(Q, 0, 1) = h-1,0 (with wrapping forbidden).

As in [23] we shall however take a quite different route and resort to a transfer matrix construction. This has the advantage of linking up more easily with CFT [29], and giving very precise numerical results for all the boundary conditions outlined above. Moreover, the ensuing transfer matrix formalism is no more complicated than the one [30] routinely used in the study of the FK clusters.

State space

Consider a strip of the square lattice of width L spins (boundary conditions will be detailed later). The basis states on which the row-to-row transfer matrix T acts contain one colour label ci per spin. By definition, one has ci = cj if and only if σi = σj (i.e., the two spins

Σ B L - ≈ ∼ { } ∈ { } -
on sites i and j have the same colour). The colour labels ci contain less information than the spin colours σi themselves. For instance, any configuration in which the first and third spins have the same colour, no matter which one, and no other spins have identical colours, is represented by

. ( 7
)
With these conventions we are always able to recognize whether two spins have the same colour or not, even if we do not know the precise value of this colour. But this is all that is needed to determine the Boltzmann weights in [START_REF] Wu | Exactly Solved Models: A Journey in Statistical Mechanics[END_REF].

For a row of L = 3 vertices, and for any (non-integer) Q, there are precisely five basis states

. ( 8
)
Note that the last state would carry zero weight for Q = 2, but apart from that the number of basis states for any given L will be finite and independent of Q.

In general, the number of states is equal to the Bell number BL of partitions of L objects, with exponential generating function Let us take the time evolution direction to be upwards, so that horizontal (resp. vertical) edges are space-like (resp. time-like). We can write T as a product of elementary transfer matrices, each represented symbolically as a rhombus surrounding a single lattice edge, and corresponding to the addition of that edge to the lattice. This edge links spins (shown as solid circles) on diametrically opposite sites of the rhombus. On an L = 4 square lattice with periodic boundary conditions this reads T = .

A rhombus corresponds to a vertical edge, and acts on a basis state s as follows. If exactly Qs distinct colour labels {ck} are used in s, then the new colour label c J i of the spin σi can be unchanged c J i = ci (with weight e K ), or any one of the other labels already in use c J i = ck (each with weight 1), or a new one c J i / ck (with weight Q Qs). Note that this latter weight is in general non-integer, and is responsible for the correct computation of the chromatic polynomial χ G ˆ(Q). For an example with Qs = 3 this reads explicitly

= + + + 1 -δc 1 ,c J 3 -δc 2 ,c J 3 -δc 3 ,c 3 J (Q -Qs) .
The last state could be written as , but this is equivalent to in the basis (8).

v = ⎝ ≤ ≥ - • • • • ⎛ ⎜ ⎞ ⎟ 0 Q 2 Q 2 0 e K + Q 3 ⎠ ⎟ ⎠ K ⎞ ⎝ ⎠
A rhombus adding a horizontal edge between vertices i and i+1 corresponds simply to a diagonal matrix, with a weight e K if ci = ci+1, and 1 otherwise.

With these rules at hand, one can write the periodic (cylinder geometry) L = 3 transfer matrix for arbitrary Q in the basis (8) as an instructive example:

T = h1 • h2 • h3 • v1 • v2 • v3 with e K 0 0 1 0 0 e K 1 0 1 1 0 1 e K 0 1 Q -1 0 0 e K + Q -2 0 -
-and h1 = diag(e K , 1, e K , 1, 1). The remaining matrices can be obtained from those given by cyclic permutations of the colour labels:

⎛ e K 1 0 0 0 ⎞ ⎜ Q -1 e K + Q -2 0 0 0 ⎟ 2 ⎜ ⎝ ⎛ 0 0 e K 1 1 0 0 1 e K 1 0 0 Q -2 Q -2 e K + Q -3 e 0 1 0 0 0 e K 0 1 1 (10) v3 = ⎜ Q -1 0 e K + Q -2 0 0 ⎟ 0 1 0 e K 1 0
Q -2 0 Q -2 e K + Q -3 and h2 = diag(e K , 1, 1, e K , 1), h3 = diag(e K , e K , 1, 1, 1).

To get a non-periodic strip geometry instead of a periodic cylinder, one needs simply to omit the last horizontal edge which is responsible for the periodic boundary conditions. For instance, with L = 3 one would have T = h1 h2 v1 v2 v3. This corresponds to free boundary conditions. Mixed boundary conditions are obtained by constraining the spins on the left boundary of the strip to belong to a subset of Q1 different colours, with 0 < Q1 Q. Spins on the right boundary remain free. For Q1 integer this can be coded in the transfer matrix by considering that the number of colours in use in a given state is always at least Q1 (and so in particular Qs Q1). These first Q1 colours are considered fixed, whereas other colours Q1 + 1, Q1 + 2 ,... are defined only relative to the fixed ones, as described above. In particular, the choice Q1 = 1 corresponds to fixed boundary conditions.

The leading eigenvalue Λ0 of the transfer matrix T with periodic boundary conditions gives the ground state free energy f0 = (1/L) log Λ0. This f0 coincides precisely with that of the usual FK transfer matrix [30], even when Q is non-integer 5 . Its finite-size corrections possess a universal L -2 term whose coefficient determines the central charge of the corresponding CFT [29].

Correlation functions

However, to obtain the desired two-point correlation functions of DW, we need to construct a variant transfer matrix Tl 1 ,l2 which imposes the propagation of l1 thin and l2 thick DW 5 The reader may wish to check this fact on the L = 3 example given above.

J. Stat. Mech. (2010) P12026

v = - { }
along the time direction of the cylinder (or strip). From its leading eigenvalue Λl 1 ,l2 one can determine the energy gap Δfl 1,l2 = (1/L) log(Λl 1,l2 /Λ0) whose finite-size scaling in turn determines the critical exponents h(Q, l1, l2) and h ˜(Q, l1, l2) [29,5].

To this end, the basis states (8) need to be endowed with some additional information about the connectivity of the spin clusters, ensuring the propagation of the desired number and type of DW. The crucial point is that we need to know whether two spins having the same colour label ci also belong to the same cluster. Thus the states that we use in the final transfer matrix have the form (11) meaning that two spins belong to the same spin cluster if and only if they are linked up in this pictorial representation. Of course, only spins with a common colour label can be linked up. Thus, in the left state of ( 11), the spins on vertices 3 and 4 are in the same cluster, but not in the same cluster as the spin 1. In the right state the spins 1, 3 and 4 are all in the same cluster. These two states are different. In the transfer matrix evolution, each time two neighbour vertices correspond to the same colour, the corresponding clusters are linked up.

Note that the possible ways of linking up the spins must respect planarity. For instance, in the L = 4 state with colour labels 1, 2, 1, 2 , one cannot simultaneously link the first and third spins, and the second and fourth, since this possibility is disallowed by planarity.

To construct Tl 1 ,l2 we modify the linked basis states (11) so that precisely l1 + l2 spin clusters are marked. The colour labels of the marked clusters must respect the chosen values of l1 and l2. In order to conserve the marked clusters in the transfer matrix evolution, none of the marked clusters must be 'left behind', and two distinct marked clusters must not be allowed to link up. When a marked and an unmarked cluster link up, the result is a marked cluster.

To summarize, the final transfer matrix thus keeps enough information, both about the mutual colouring of the sites and about the connectivity of the clusters, to give the correct Boltzmann weights to the different configurations, even for non-integer Q, and to follow the evolution of the boundary of a particular set of clusters. These boundaries are precisely the domain walls (see figure 3).

Numerical results

We have numerically diagonalized the transfer matrix in the DW representation for cylinders and strips of width up to L = 11 spins. We verified that the leading eigenvalue Λ0 in the ground state sector coincides with that of the FK transfer matrix, including for non-integer Q. As to the excitations Λl 1,l2 , we explored systematically all possible colouring combinations for up to four marked spin clusters, for a variety of values of the parameter κ, and for several different boundary conditions (periodic, free, fixed, and mixed).

Finite-size approximations of the critical exponents h(L) and h ˜(L) were extracted from the leading eigenvalue in each sector, using standard CFT results [29,5], and fitting both for the universal corrections in L -2 and the non-universal L -4 term. These

- | | - | | - - →∞ ≤ Table 1.
Bulk critical exponents corresponding to five different DW configurations (l1, l2), as functions of the parameter p = κ/(4 κ), along with the conjectured exact expression (4). The last column labelled (0, 1) * is for a single DW which is allowed to wrap around the periodic direction. The table entries give the value of ρ , when ( 5) is rewritten as hr,s = (ρ 2 1)/(4p(p + 1)), with error bars shown in parentheses. approximations were further extrapolated to the L limit by fitting them to first-and second-order polynomials in L -1 , gradually excluding data points corresponding to the smallest L. Error bars were obtained by carefully comparing the consistency of the various extrapolations.

p = κ/(4 -κ) (2, 0) (0, 2) (3, 0) (2, 1) (0, 3) (0,
Representative final results for the bulk case (periodic boundary conditions) are shown in table 1. Recall that (4) is only valid for (l1, l2) = (0, 1) if the spin cluster is forbidden to wrap around the periodic direction. The sector where such wrapping is allowed is denoted as (0, 1) * in table 1, and the corresponding critical exponent turns out to be h0,1/2.

Final results for the boundary case with free boundary conditions are shown in table 2. The convention for the labels (l1, l2) is the one stated below equation (6).

Mixed boundary conditions

We recall that mixed boundary conditions have been defined by imposing that the spins on the left boundary of the strip belong to a subset of Q1 different colours, with 0 < Q1 Q, while the spins on the right boundary remain free. In particular, Q1 = 1 (resp. Q1 = Q) describes fixed (resp. free) boundary conditions. As above, the convention for the labels (l1, l2) is the one stated below equation (6).

An important point is that the fixed colour of the leftmost marked spin cluster may or may not be one of those allowed on the left boundary. For each choice of colours of the marked clusters, we define a sign ϵ = 1 (resp. ϵ = 1) if the colour of the leftmost cluster is k1 Q1 (resp. k1 > Q1), i.e., if the leftmost cluster is allowed to (resp. not allowed to) touch the left boundary. This sign closely parallels a construction in [24], where the case ϵ = 1 (resp. ϵ = 1) was referred to as the blobbed (resp. unblobbed) sector. For free boundary conditions, only the choice ϵ = 1 is meaningful (at least for 0 < Q1 Q, although it might be possible to give a meaningful definition of the model outside this range, by a suitable analytical continuation). The issue of the sign ϵ = 1 is illustrated in figure 5.

One would expect the mixed boundary conditions to be given by the conformal fusion of the operator Φ1+2(l 1 -l2),1+4l1 that inserts the required number of DW, and the free-tomixed boundary condition changing operator Φr 0 ,s0 . The latter operator is responsible for the shift from Q to Q1 allowed colours on the boundary. It has been worked out in [24] in the context of a different set of geometrical observables (TL loops), but since it should be representation independent we can take it over here. The dominant exponents for mixed boundary conditions therefore follow from the standard CFT fusion rules [5] as

h r0+2є(L1-L2),s0+4єL1 , (12)
where hr 0,s0 is the dimension of the free-to-mixed boundary condition changing operator Φ r0,s0 . The values of (r0, s0) follow from [24], e.g., ( 1, 2) for Q1 = 1 and ( 1/2, 1) for Q1 = 2. We have verified these and other cases by explicit numerical calculations. In particular, the possible ambiguity of ϵ versus ϵ in ( 12) is ruled out by the numerical results. Moreover (12) agrees with (4) for free boundary conditions (ϵ = 1) as it should.

A massless scattering description

While it is far from obvious how to derive (4) by CG methods, minor progress can be achieved using a rather different set of ideas. We restrict the discussion in this section to the bulk properties.

We can learn about the dynamics of DW in the critical theory by using known information about the low-temperature (K > Kc) phase of the Potts model. Although nonintegrable on the lattice, the corresponding deformation by the operator Φ21 is integrable in the continuum [31]. It can be described using a basic set of kinks Kab These kinks scatter with a known S-matrix related to the BWM algebra [22]. Importantly, the dynamics conserves the number of kinks: the process KabKbc Kac is forbidden (as in any elastic relativistic scattering theory), although kinks do appear as bound states in kink-kink processes. Many properties of these kinks can be calculated using integrability techniques. When the mass m 0 (i.e., K Kc), the S-matrix provides a 'massless scattering' description [START_REF] Fendley | Proc. Trieste Summer School in High Energy Physics and Cosmology[END_REF] of some of the degrees of freedom of the critical theory itself. It is not entirely clear what a kink, which is well defined for K Kc, becomes at Kc, but it is natural to expect thick DW to be described by the propagation of two (or more) kinks such as KabKba. As for thin DW, the potential existence of regions where they are reduced to a single edge-that is a single kink-suggests that they have to do with more processes involving two kinks merging into one. We will not discuss them further here.

It is an easy exercise to obtain the scaling dimension of thick DW using the massless scattering description. Indeed, the fact that the S-matrix satisfies relations from the BWM algebra allows us to re-express it in terms of the a (2) or Bullough-Dodd S-matrix [START_REF] Efthimiou | [END_REF],

for which the thermodynamic Bethe ansatz was studied in [34]. This describes also the dynamics of the field theory

1 S = 8π d x[(∂x Φ) 2 + (∂ y Φ) 2 + g(2e -(i/ √ 2)βΦ √ +e )] (13)
with β 2 = κ/4. Giving each kink the fugacity Q -1 produces the correct central charge

3 (κ -4) 2 c = 1 -2 , ( 14
)
κ
where each kink has a U (1) charge equal to 0, ±1. The scaling of the sector with charge j produces a gap Δj = j 2 /4κ, so the leading dimension is Δj -(1c)/24 = hj/2,0. This agrees with h-l 2,0 with l2 = j/2, so there are two kinks per thick DW.

This argument validates (4) for l1 = 0.

Discussion

The results ( 4) and ( 6) can be used to predict the fractal dimension of various geometrical objects related to spin clusters. It should be possible to observe these dimensions in Monte Carlo simulations, and possible in real experiments.

Fractal dimensions

The dimension of a spin cluster follows from (6) as

⎧ ⎨ 2 for κ ∈ [2, 8 ] d = min(2, 2 -2h0,1/2) = (8 + κ)(8 + 3κ) (15)
in agreement with [10,11].

⎩ 32κ for κ ∈ [ 3 , 4] b b b b ≥
The boundary of a spin cluster has dimension

κ db = 2 -2h-1,0 = 1 + 8 . ( 16
)
This is an example of a duality relation: db for a spin cluster with SLE parameter κ equals the dimension d FK of the boundary of an FK cluster at the dual parameter κ * = 16/κ. It is moreover known that d FK at parameter κ coincides with the dimension d EP of the FK cluster's external perimeter (with fjords filled in) at parameter κ * [17]. Combining these, (16) means that db = d EP , i.e., the dimension of the boundary of a spin cluster equals that of the external perimeter of an FK cluster, at the same κ. This latter fact has very recently been verified numerically by Monte Carlo simulations, including for non-integer values of Q [28].

In the half-plane geometry, the intersection of the spin cluster containing the origin with the real axis has the dimension

ds = min(1, 1 -h3,5 ) = min 1, 8 8 -κ - 3κ , (17) 2
as follows from (4).

Thin and thick domain walls

We have initially attached the epithets 'thin' and 'thick' to the two kinds of domain walls on the basis of considerations in the microscopic model. Indeed, thin DW can narrow down to a single lattice edge on the dual lattice, separating spins of different colour. Since (4) implies that the two kinds of DW scale differently in the continuum limit, one would expect to find that the continuum geometrical objects can also somehow be characterized as 'thin' and 'thick'. We now show that this is indeed the case.

According to (4), the set of points where a thin DW has minimal width-i.e., one lattice spacing in the microscopic formulation, or the two clusters separated by the DW 'come close' in the continuum limit in the sense of the small neighbourhoods defined previously-has dimension

3 d1 = 2 -2h2,4 = 8κ (4 -κ)(5κ -4). ( 18
)
This is analogous to the dimension of so-called 'red' or 'pivotal' bonds in the theory of percolation and of FK clusters. The corresponding result for a thick DW reads

d2 = max(0, 2 -2h-2,0) = 0. ( 19
)
The fact that d1 0 means that thin DW are indeed thin, in the sense that they have a macroscopic number of loci of zero thickness in the continuum limit. Likewise, the fact that d2 = 0 means that the loci where thick DW have zero thickness is a set of measure zero.

To summarize, there is therefore a consistency between the distinction of the two kinds of DW microscopically and in the continuum limit, and this ensures that they can (and do) scale differently.

One can also obtain dimensions analogous to d1 and d2 at the boundary. Indeed, define d ˜1 (resp. d ˜2) to be the fractal dimension of the set which is the intersection of the real axis with the set of points where a thin (resp. thick) DW has minimal width. In other

κ → → → → → → → → ⎪ ⎨ 5 1 (3 -κ)(5κ -8) for κ ∈ [ 12 , 3]
words, these are the dimensions of pivotal bonds on the boundary. By setting l1 = 2, l2 = 0 and l1 = 0, and l2 = 2 in (4) we have then, respectively,

⎧ ⎪ 1 for κ ∈ [2, 12 ]
and

⎪ ⎩ 0 for κ ∈ [3, 4] d ˜2 = max(0, 1 -h-3,1) = 0. ( 21
)
7.3. The limit Q = 4

In the limit Q 4 (or κ 4), we have d1 0 and d2 0. This means that the distinction between thin and thick DW disappears in that limit. This is of course consistent with the fact that κ = 4 is a fixed point of the duality transformation κ * = 16/κ, and so spin clusters and FK clusters have identical properties.

A further manifestation of the indistinguishability of thin and thick DW is that the bulk exponent ( 4

) becomes (l1 + l2) 2 h l1-l2,2l1 = . ( 22) 4
Indeed, this formula is invariant under the permutation of l1 and l2. Moreover, the exponent for l FK clusters reads l 2 /4 for Q = 4 [5,[START_REF] Nienhuis | Loop Models[END_REF], and so a spin DW (whether thin or thick) behaves as a FK cluster, as expected.

The case Q = 2

For the Ising model Q = 2 (or κ = 3) the absence of branchings means that one thick DW equals two thin DW. Indeed the bulk exponent for (4) becomes 4l 2 -1 h l1-l2,2l1 = (23) 48 with l = l1 + 2l2 in that case, and this agrees with the exponent hl/2,0 for l loop strands [5,[START_REF] Nienhuis | Loop Models[END_REF] in the dilute O(1) model.

The limit Q → 1

The limit Q 1 (or κ 8/3) is often studied in the context of FK clusters, since these become then bond-percolation clusters. Setting Q = 1 for spin clusters has a more trivial meaning: all the spins are simply in the same state, σi = 1. Therefore we should have d = 2 and ds for Q 1. This is indeed satisfied by (15) and (17).

Apart from that, the Q 1 limit for spin clusters is more subtle. One can tackle it by considering configurations of bond percolation on the square lattice where one picks out one single percolation cluster. The spins on the vertices sitting on this cluster are given a certain colour. All the other spins in the system are given another single colour. Then the boundary of the coloured spin cluster is a single loop with fractal dimension df = 4/3, as expected from the duality κ * = 16/κ. It is, however, remarkable that we )) = (20)

J. Stat. Mech. (2010) P12026

d ˜1 = min(1, max(0, 1 -h5,9

≤ ≤ can find such a simple construction of the spin cluster for Q = 1 on the lattice. For more interfaces one has to pick out more percolation clusters, and give them different colours. Finally, we note that for Q = 1 we have l 2 -1 h l1-l2,2l1 = , (24) 24 which agrees with h0,l/2 if one sets l = l1 +3l2. This seems to indicate that for percolation, a thick DW equals three thin DW (or TL loop strands)-a curious result for which we have no convincing explanation at present.

Conclusion

We have defined a set of geometrical observables based on the branching domain walls of the Potts model for any real 0 Q 4. These observables are defined starting from the DW expansion (3) of the partition function. We have studied these objects numerically for noninteger Q using a transfer matrix formulation. Our results are compatible with conformal invariance of these observables and we have given a set of conjectures for the corresponding bulk and boundary critical exponents.

It appears straightforward to extend our transfer matrix construction to deal also with the tricritical Potts model, by including vacancies with the appropriate fugacity. We strongly expect our expressions (4) for the critical exponents to carry over to that case, via the usual analytic continuation of the Coulomb gas coupling constant: simply apply (5) with 4 ≤ κ ≤ 6 instead of 2 ≤ κ ≤ 4.

These new families of exponents for the tricritical model shou √ ld then also apply (by universality) to the dilute phase of the O(n) model with n = Q, and (by the SLE duality [17,16]) to FK clusters in the usual critical Potts model. It is however not clear at present how to identify the corresponding microscopic observables in the well-known lattice realizations of the dilute O(n) model and the critical Potts model.

Figure 1 .

 1 Figure 1. A configuration of the Q = 3 Potts model, and the corresponding set of branching domain walls.

L

Figure 2 .

 2 Figure 2. A domain wall configuration corresponding to a graph G, and its dual graph G ˆ.

Figure 3 .

 3 Figure 3. The two different kinds of DW, here shown in the bulk case (geometry of the infinite plane). A thin DW corresponds to the interface between two clusters of different colours (a), while for a thick DW the two clusters have the same colour. An illustration for the Q = 3 Potts model is given (left) as well as a schematic picture for non-integer Q (right).

Figure 4 .

 4 Figure 4. A thin DW in the half-plane geometry for Q = 3. The left (resp. right) panel illustrates the case where the positive real axis supports free (resp. fixed) boundary conditions.

  BL differs from the Catalan number CL of planar partitions, familiar from the FK cluster representation. In particular, 1, 2, 1, 2 is a valid state for L = 4. For L 1 one has asymptotically CL 4 L , whereas BL grows super-exponentially. When Q is integer, the number of states truncates to Q L /Q!.

Figure 5 .

 5 Figure 5. Mixed boundary conditions in the half-plane geometry for Q = 3 and Q1 = 2. The left (resp. right) panel illustrates the case ϵ = 1 (resp. ϵ = 1) where the propagating cluster, shown in blue colour, is allowed to (resp. not allowed to) touch the positive real axis.

  , i.e., ordered regions where the dominant value of the spin is a or b.

Table 2 .

 2 Boundary critical exponents with free boundary conditions,

		2, 0)	(1, 2)	(2, 1)	(3, 0)
	3	7.00(3) 19.0(2) 15.16(4) 11.13(4)
	4	10.89(7) 24.8(2) 20.9 (1) 16.74(5)
	5	14.5 (1) 31.1(2) 26.7 (2) 22.2 (3)
	Exact	4p -5	6p + 1	6p -3	6p -7

corresponding to four different DW configurations (l1, l2), as functions of the parameter p = κ/(4 κ), along with the conjectured exact expression (4). The table entries give the value of ρ , when (

5

) is rewritten as hr,s = (ρ 2 1)/(4p(p + 1)), with error bars shown in parentheses. p = κ/(4κ) (

doi:10.1088/1742-5468/2010/12/P12026

Acknowledgments 17

Acknowledgments

We thank G Delfino and P Fendley for discussions. This work was supported by the Agence Nationale de la Recherche (grant ANR-06-BLAN-0124-03).