Quantile-constrained Wasserstein projections for robust interpretability of numerical and machine learning models

Robustness studies of black-box models is recognized as a necessary task for numerical models based on structural equations and predictive models learned from data. These studies must assess the model's robustness to possible misspecification of regarding its inputs (e.g., covariate shift). The study of black-box models, through the prism of uncertainty quantification (UQ), is often based on sensitivity analysis involving a probabilistic structure imposed on the inputs, while ML models are solely constructed from observed data. Our work aim at unifying the UQ and ML interpretability approaches, by providing relevant and easy-to-use tools for both paradigms. To provide a generic and understandable framework for robustness studies, we define perturbations of input information relying on quantile constraints and projections with respect to the Wasserstein distance between probability measures, while preserving their dependence structure. We show that this perturbation problem can be analytically solved. Ensuring regularity constraints by means of isotonic polynomial approximations leads to smoother perturbations, which can be more suitable in practice. Numerical experiments on real case studies, from the UQ and ML fields, highlight the computational feasibility of such studies and provide local and global insights on the robustness of black-box models to input perturbations.


Introduction
Multiple engineering fields require models for prediction and phenomenological understanding. Machine learning (ML) and uncertainty quantification (UQ) of numerical models are two essential approaches to developing and manipulating such models. Because they require, for their enlightened use, an adequate understanding of their characteristics, they share fundamental similarities. These two frameworks feed off each other through the duality of sensitivity analyses (SA), a fundamental methodological corpus in UQ, and ML interpretability methods, as explained by [90,61]. In particular, recent advances in explainable ML leverage tools from SA to produce meaningful interpretations of black-box models [36,13], and novel SA estimation schemes are heavily based on the construction of suitable ML models [18,12]. Both SA and ML interpretability especially rely on the definition, estimation, and manipulation of diagnostics related to the characteristics of a model and how its behavior depends on its inputs [28,74,95]. Formally speaking, let a model f be defined as a mapping between inputs X ∈ X and outputs Y ∈ Y where (X , Y) are two metric spaces: In a ML context, f is defined as a predictive model (e.g., penalized linear regression, neural network) linking a feature (or covariate) instance X to a prediction Y [54]. In the UQ framework, a so-called computer model f represents the numerical implementation of a hypothetical-deductive link (e.g., by systems of ordinary differential equations, by finite element methods) between X and Y [100].
In both fields, the input X ∈ X is generally assumed to be random, inducing a general framework for handling uncertainties about the latter. Let P be the distribution of X. In the ML context, P is defined implicitly by an empirical measure: given a set of observations x (1) , . . . , x (n) ∈ X , where δ denotes the Dirac measure. On the other hand, in the UQ setting, P is often explicitly chosen based on observations of X, expert assessment (domain knowledge), or stochastic inversion from observations of Y [105]. The diagnostics mentioned above correspond to estimations of key interpretable features of Y , or quantities of interest (QoI). In the SA literature, such a QoI is often referred as the score [94], while ML researchers rather speak about predictive performance [81]. Recall that SA aims to rank the dimensions of X according to their influence on this QoI [20]. For instance, in local SA, it is usually computed by way of a differential operator with respect to (w.r.t.) the dimensions of X [28]. In global SA, it can be typically chosen as the output's variance or its quantiles [36,72]. More general objects characterizing the distribution of Y (e.g., a kernel embedding [9]) can also be of interest, possibly at the cost of a less immediate interpretability. In ML interpretability, local methods focus on a particular prediction instance, letting the QoI be the identity function [104] or a local linear decomposition of f [113]. Global ML interpretability often relies on QoI defined as performance metrics (e.g., accuracy, loss value) of f computed over a training dataset [49,26,60]. Accordingly, the use of these diagnostics to allow different levels of interpretability is subject to the same robustness problem: they must remain relevant when P suffers from misspecification. It would improve the confidence in both the usage and the insights that an ML model offers [10]. In this framework it is fundamentally connected to problems of domain adaptation and transfer learning [19,70] (e.g., when the data used for the design of f suffer from selection bias w.r.t. operational data), including in particular the robustness to covariate shift [51,109] (see Figures 1 and 2 for an illustration). More generally, whether in UQ or ML, this misspecification is due to the epistemic uncertainty affecting the knowledge of the properties of P (e.g., support, geometry, topology) due to the finiteness of the available information (e.g., data, expertise, boundary conditions) [58,105].
This epistemic uncertainty about P thus impacts the generalization capability of f . It is essential : Two-dimensional illustration of covariate shift in the UQ framework. Graph (a) displays the isodensity curves of a joint distribution P on X = (X 1 , X 2 ). Graph (b) displays the result of a simultaneous quantile shrink for X 1 and quantile dilatation for X 2 while preserving the dependence structure.
to define a practical application domain [93], or usage domain of f , allowing to extrapolate beyond a validity domain based on selected situations (e.g., a domain including unobserved data points in ML, or unsure structural mechanisms in numerical models). Generalization bounds in ML usually depend on notions of capacities (e.g., the Vapnik-Chervonenkis dimension [111]). They rely on a sufficient number of observations to obtain generalization, i.e., linked to epistemic uncertainty reduction. Hence, generalization power and epistemic uncertainties are intimately linked [110]. Therefore, robustness studies of the application domain's uncertainty are particularly important for both fields. It motivated the search for adversarial situations in ML [4] and the use of imprecise probabilities in numerical model exploration [11]. The need to conduct robustness studies raises the question of how to reasonably vary P in a perturbation class C ⊂ P(X ), where P(X ) is the set of probability measures supported on X , by modifying the mass given to different parts of X in an interpretable way. Such variations allow checking the consistency of the information delivered by these diagnostics. Here, reasonableness means two things: a perturbed distribution Q must remain close to the original P in some formalized sense, and the computations required for the robustness analyses must have a moderate cost.
Several authors in both fields have conducted this type of study. In the ML setting, finding adversarial situations or detecting bias affecting a QoI can stand on moving the original data (e.g., by optimal transport [89,47]), resulting in an empirical shift of P . Researchers in UQ distinguish between distributional sensitivity analysis (DSA) [2,76] which aims at providing a robustness diagnosis on a QoI used in SA (or connected similar studies, e.g., in financial risk analysis [25]), and other robustness approaches whose goal is to study the sensitivity of a given QoI to a gradual modification of P in C. For instance [69], [43] and [63] proposed such indices for exceedance probabilities, quantiles or superquantiles of the distribution of Y under various perturbation schemes.
Many choices have been proposed for perturbation classes. They are often indexed by parameters θ and denoted C(θ). For instance, in the SA setting, they can be defined as contamination classes for Bayesian SA [50,92], as classes of generalized moments [69,6] or as classes defined by Fisher-Rao derivatives [67,43], among others. A unified view of these approaches is to define the perturbed distribution Q as the solution to the projection problem where D is a discrepancy between probability measures and C(θ) denotes a fixed perturbation class. The Kullback-Leibler (KL) divergence is often chosen as a suitable discrepancy (e.g. in [69,6]), leading to entropic projections [27]. For instance in those developments, whenever X ⊂ R d (d ≥ 1), C(θ) can be defined as the distributions belonging to P(X ) with a deviated mean η, different from E P [X]. Reusing the term proposed by [6], θ = η − E[X] can be called a perturbation intensity, where θ belongs to an ordered set Θ. Hence, for a fixed perturbation intensity θ, the solution Q of (2) is the optimally perturbed distribution of P w.r.t. C(θ).
Although rational, moment-based perturbations of P presents significant issues that limit the conclusions of interpretability analyses, and the subsequent unification of UQ and ML methodologies. First, the choices of θ, C(θ), and Θ remain subjective and submitted to strong assumptions. For instance, in [69], some generalized moment of Q are required to exist, and the moment deviation intensity must subjectively chosen by the practitioner. When dim(X) is large, such assessments seem challenging to make. In [43], Θ comprises the values taken by the derivative of the Fisher metric, requiring regularity properties on P , which is often assumed to be in a restrictive parametric family, for computational reasons. In general, determining the boundaries of Θ remains challenging.
Second, choosing the KL divergence or the derivative of the Fréchet metric as in [53] as a discrepancy imposes strong restrictions on the space of reachable probability distributions. It result in a search on a restricted part of P(X ) in (2). For instance, using the KL divergence implies that Q should be absolutely continuous w.r.t. P , which does not allow to consider continuous perturbed distributions Q given an empirical initial distribution P . Finally, the behavior of these discrepancies (purely related to information geometry) remains uneasy to explain to non-specialists.
This article, therefore, aims to address both of these issues, improving the connections between the interpretability analyses conducted in UQ and ML settings, and their robustness. More precisely, our contributions are twofold: (a) We propose a meaningful and generic approach to perturb distributions through marginal quantile constraints, selecting the 2-Wasserstein distance for D. These choices allow to solve the problems mentioned above. A key point of the proposed methodology is that no strong regularity assumptions must be assumed on f and it does not rely on f having a particular structure (e.g., belonging to a particular family of predictive models or based on specific physical equations). This effectively unifies UQ and ML approaches.
(b) We demonstrate the computational tractability of this methodology by implementing it on different types of QoI, on both numerical and predictive models, and studying some robustness indicators to perturbations in the particular case where f is considered to be a black box.
This article is organized as follows. In Section 2 we first introduce useful notations and definitions. Section 3 develops and motivates the choice of quantiles as an interpretable and generic basis for defining meaningful perturbations and defines perturbation schemes relevant for different robust interpretability studies. Section 4 presents the general framework of probability measure projection using the 2-Wasserstein distance and proposes analytical solutions and numerical optimization schemes for solving the input perturbation problem through the help of isotonic polynomials. Section 5 showcases our method on two use-cases from both the ML and UQ fields, from which local and global robustness insights are highlighted. A discussion section ends this article, opening avenues for improvement. All proofs of technical results are postponed to a dedicated appendix.

Notations and main definitions
Let us introduce some useful notations and definitions. Let d and p be two positive integers. Let X be a subset of R d , and P p (X ) be the subset of P(X ) of measures with finite p-th moment. The initial probability measure P is either defined through an explicit distribution, or empirically, as in (1). We will denote by Q ∈ P(X ) the optimally perturbed distribution of P .
Furthermore, let us denote Ω X ⊂ X the application domain. It is the subset of X where f is intended to be used for predictions [93]. Both in ML and UQ, given a set x n = {x 1 , . . . , x n } of training, validating or testing examples, the convex hull of x n or a broader span of x n are common candidates to define Ω X [116]. More generally, Ω X is an extrapolation domain where f is assumed to generalize well (e.g., a paving of a compact subspace of X selected by tree-based classification [56], confidence measures or cross-validation schemes [57,78,114,66]). In ML specifically, including out-of-distribution data in Ω X remains an open problem [55,114,99].
To echo the classical assumptions in ML and UQ, we assume that Ω X is the union of compact subsets of X . These subsets can be defined under some uncertainties, typically on their bounds. In a robustness analysis perspective, assuming that the dependence structure is maintained, the uncertainties on Ω X can be interpreted as variations on the values (or thresholds) of the extreme quantiles of marginal distributions. Figure 3 illustrates a typical situation for a univariate marginal of X.
The upcoming developments deal with perturbations on the marginal distributions of P . The following definitions recall classical results on cumulative distribution functions (cdf) and generalized quantile functions (gqf) of univariate probability measures. Denote F the space of univariate (a.) Figure 3: Application domain Ω X ⊂ X in the UQ (a.) and ML (b.) settings. In the UQ setting, X is the support of the explicitly chosen density (in green). In the ML setting, X is the interval between the minimum and maximum observed values. In both cases, Ω X is included in X , although it is not mandatory. distribution functions: For any P ∈ P(R), F contains the cdf of P , defined as: The gqf of any probability measure P can be formally defined as follows [91,33,64].
Definition 1 (Generalized quantile function). Let P ∈ P(R) with cdf F P .
(i) The generalized quantile function (gqf ) of P is the unique left-continuous, non-decreasing generalized inverse of F P , defined, for every a ∈ (0, 1), as: (ii) The unique right-continuous non-decreasing generalized inverse F → P of F P , almost-everywhere equal to F ← P , is defined, for every a ∈ (0, 1), as: It is important to note that when F P admits an inverse F −1 P in the traditional sense (e.g., it is continuous, strictly increasing), then the following equality holds: which echoes and generalizes the traditional definition of a quantile function as the inverse of a cdf.
In the remainder of this work, any probability measure is assumed to have a finite 2-nd order moment; hence, their gqf is assumed to be square integrable. Accordingly, the perturbation θ considered in the projection problem (2) corresponds to quantile constraints placed on one-dimensional marginal variables. Therefore it is assumed that Θ ⊂ R.

Quantile perturbations
In this section, we define a particular set of perturbation classes. They are characterized as a set of constraints on the (generalized) quantiles of the marginal components of the inputs (or features) variables X of f . Before presenting a formal definition of such classes, we argue their relevance for practical studies.

Motivations
First, generalized quantiles always exist, unlike generalized moments considered by [69]. It is also what motivated [6], who proposed to define the perturbation range Θ in (2) as a compact set bounded by empirical quantiles. More precisely, recalling Definition 1, a gqf is the left-generalized inverse of a cdf. Moreover, they allow to uniquely characterize probability measures in P(R). Remark 1. The unique link between left-generalized inverses of functions in F and probability measures can be formalized as follows. Let F ← be a function in the set: Then there exists a unique probability measure in P(R) with cdf F , and such that F ← is its gqf. This classical result seems to be widely known in the literature. For completeness, a proof sketch is provided in the appendices.
Since generalized inverses of functions in F always exist, perturbing marginal quantiles in (2) do not require additional assumptions either on the initial probability measure P or on the shape of the target probability measure Q. Hence, it allows for generic, well-defined perturbations, which is key for merging both SA and ML interpretability. Second, constraints placed on (generalized) quantiles can be interpretable. Considering b i ∈ R as a representative magnitude of a real component X i of X, then P(X i < b i ) = α i ∈ (0, 1) if and only if b i minimizes the expected pinball cost function [23]: This means that α i can be interpreted as the ratio c 1 /(c 1 + c 2 ) where c 1 |x i − b|1 {Xi<b} and c 2 |x i − b|1 {Xi≥b} are relative first-order costs associated to estimating x i by b. When X i is an influential input, it seems relevant to produce rules to modify these costs (then the value of α i ), shifting the distribution of X i . In particular, this information-theoretic interpretation of quantiles promotes their variation in sensitivity studies conducted in economic or financial contexts (see, e.g., [85]). Moreover, Bayesian statisticians have long recognized this representation as one of the most appropriate formal approaches to incorporating expert knowledge into statistical model inference [44,73]. Third, in many applied problems, (generalized) quantile specifications are often key to studying the influence of input variables on a decision-making output. In both the UQ and ML framework, inputs X can themselves be partially derived from calculations from upstream learning or numerical models (e.g., for multi-physics problems). P is then often calibrated by quantile matching, which may introduce uncertainties on some of its marginal features [103]. Numerous applications dealing with economic stress tests or risk mitigation against natural hazards use quantiles as influential inputs of decision-helping models. For instance, in the drought risk studies in [35], the association between soil wetness, climatic, seismic, and socioeconomic variables (e.g., city-level description) is often carried out using marginal quantiles that play the role of features for cost predictive models. Input variations of daily value-at-risk percentiles, computed from legacy data, were recently required by the European Banking Authority for generating macroeconomic scenarios used for EU-wide stress tests [5]. Reverse SA studies for financial risk management, such as those conducted in [86], are primarily based on moving values-at-risk, which are quantiles.
Let us end this subsection with two motivating examples. They offer two additional concrete illustrations of using quantiles for influence analysis. They also illustrate two different quantile perturbation schemes: quantile shifting and application domain dilatation. These schemes are formally introduced in Section 3.3.
Example 1 (Economic stress test). Inspired by [16], assume that an economic shock happens in an abstract country. Prospective analyses announce a $200 drop in the population median wage. Before the shock, the population wage distribution P is known (or observed), thanks to some annual census data. This distribution has a median wage of $2000. The new population wage distribution is unknown due to the lack of recent data. The economists would like to know if they can be confident in their predictive macro-economic model f w.r.t. this sudden change. A way to answer this problem would be assessing the behavior of the model f on a distribution Q, such that: Example 2 (River water level). This example is inspired from [62] and more deeply studied in Section 5.2. The safety of an industrial site located near a river is studied through the prediction of the water level Y = f (X) where f is a numerical hydrodynamic model, and X gathers the physical features of the river. A key dimension of X is the Strickler roughness coefficient for the upstream water level [42], which is modeled as a truncated Gaussian distribution on Ω X = [20,50]. However, this application domain is tainted with epistemic uncertainties on the actual nature of the riverbed (e.g., more or less subject to shrubby vegetation). The practical use of f would require assessing its predictive power under a wider interval Ω X = [5,65]. A way to express this prospective study is to assess the model's behavior on a distribution Q, such that:

A formal definition of quantile perturbation classes
Focusing on a univariate input X ∼ P ∈ P(R), let us first recall the formal definition of a quantile. For P ∈ P(R) and X ∼ P , for α ∈ [0, 1], an α-quantile of P is a number p α ∈ R such that: In certain cases, an α-quantile is not unique. For instance, assuming that P is purely atomic (e.g., an empirical measure) and that its cdf F P takes the constant value α on an open interval (t 0 , t 1 ) (i.e., it is the case if t 0 and t 1 are both atoms of an empirical probability measure), then any t ∈ (t 0 , t 1 ) is an α-quantile. By convention, the gqf of P defined by (4) is the smallest of the α-quantiles of P (in this case, F ← P (α) = t 0 ). As a result, given a chosen b ∈ X , defining a perturbed version Q of P through the equality constraints F ← Q (α) = b seems somewhat arbitrary. It would implicitly imply constraining the smallest α-quantile value of Q. Arguably, the value of b should be constrained to be a part of the set of all possible α-quantiles: In the case where F Q is invertible, it becomes a traditional equality constraint: any α-quantile is uniquely defined (i.e., F ← Q (α) = F → Q (α)). Constraints of the form (7), which are referred to as quantile constraints. They are the basis to define quantile perturbation classes. Definition 2 (Quantile perturbation class). Let K ∈ N * be the cardinality of a collection of quantile constraints defined by α = (α 1 , . . . , α K ) ∈ [0, 1] K and b = (b 1 , . . . , b K ) ∈ R K . The quantile perturbation class Q ⊂ P(R) is the set of probability measures defined as: Under simple conditions on the values α and b, quantile perturbation classes are non-empty.
Lemma 1. Let Q be a perturbation class defined over α ∈ [0, 1] K and b ∈ R K , which are assumed to be ordered, without loss of generality. If then Q is non-empty.
Since F ← Q can oftentimes be discontinuous, smoothing restrictions can be envisioned. It entails restricting Q to probability measures whose quantile functions are smooth (e.g., continuous, derivable).
Smooth quantile perturbation classes can be introduced as follows.
Definition 3 (Smooth quantile perturbation class). Let K ∈ N * and let α = (α 1 , . . . , Additionally, let V ⊆ F ← be a given set of smooth non-decreasing functions. The smooth quantile perturbation class Q V ⊂ P(R) is the set of probability measures defined as: These classes are further discussed and illustrated in Section 4.2.2. Note that smooth perturbation classes generalize perturbation classes since Q = Q F ← . Therefore, without loss of generality, perturbation classes are denoted Q V in the following. In the next few paragraphs, echoing the illustrative examples in Section 3.1, we formalize two types of quantile constraints: quantile shifts and application domain dilatations.

Two key quantile perturbations 3.3.1. Quantile shift
The quantile shift perturbation aims at increasing as well as decreasing some of the quantiles of the initial distribution. Formally, given a quantile level α ∈ [0, 1], and an initial α-quantile p α = F ← P (α), the quantile shift entails defining values for b in (7) ranging over a compact interval [η 0 , η 1 ] ⊆ Ω X such that η 0 < p α < η 1 . The next lemma formalizes an expression for b as a function of a perturbation intensity θ standardized on Θ = [−1, 1].
We refer to Figure 4 (a.) for an illustration of this perturbation scheme, and to Section 5.1.1 for a real-world application. The quantile shift perturbation class can, for a given initial quantile level α ∈ [0, 1], and valid interval bounds η = (η 0 , η 1 ), η 0 < p α < η 1 , be formally defined as the collection of perturbation classes indexed by the intensity θ ∈ [−1, 1]. Each choice of θ induces a perturbation class for which the projection problem in (2) must be solved.

Application domain dilatation
Application domain dilatation consists in perturbing the bounds of the application domain of an input. For a univariate X ∼ P with Ω X = [ω 0 , ω 1 ], the dilatation process amounts to widening or narrowing the width (or diameter diam(Ω X )) of Ω X . It is similar to perturbing extreme quantile levels (α ∈ {0, 1}) of P while preserving the midpoint of Ω X . The dilatation is characterized by a parameter η > 1 controlling the rescaling magnitude of Ω X while preserving its midpoint. In other words, one aims at finding a distribution Q with support is equal to the midpoint of Ω X , but such that diam(Q) := diam(Supp(Q)) is rescaled compared to diam(Ω X ). Similarly to quantile shift, the next lemma formalizes expressions for these two bounds as a function of a perturbation intensity θ standardized on Θ = [−1, 1].
We refer to Figure 4 (b.) for an illustration of this perturbation scheme. The initial application domain is displayed in magenta and is subject to a dilatation of parameter η = 2. The red constraints halve its width, and the blue constraints double it. One can additionally check that in both cases, the midpoint of the original validity domain is preserved. Section 5.2.1 showcases the usage of application domain dilatation perturbation in practice.
Given a perturbation class Q V defined over some modeling constraints, a collection C V (θ) of perturbation classes driven by θ can be easily defined: in the case of application domain dilatation perturbations, or T (η, θ) as defined in (9) in the case of a quantile shift.
Many perturbation settings can be defined by combining quantile shifts and domain dilatations. However, for the sake of simplicity, quantile shifts and domain dilatations are studied independently in Section 5.

Perturbing multiple inputs
In the previous sections, perturbation classes have been defined marginally, i.e., on uni-dimensional probability measures. Whenever multiple inputs (or features) are involved, an independence assumption is often assumed, allowing to perturb marginal distributions independently, as proposed by [69]. While this hypothesis facilitate the interpretation and use of models, it is questionable in practice. Therefore, one of the main challenges in ML interpretability and SA is to account for the potential dependence structure between the inputs (or features) [88]. Dependencies can often provide useful information which must be preserved through a simultaneous perturbation (e.g., to maintain the plausibility of perturbed information and avoid creating meaningless patterns [14]). In ML, perturbation problems dealing with data privacy are thus particularly concerned by the preservation of dependencies [71,84]. In UQ and ML frameworks, marginal quantile perturbations on P must obey this requirement: the dependence structure between the initial and perturbed distributions must remain the same.
Dependencies between random variables can be modeled using copula-based representations [77]. Let us recall the definition of a copula (or dependence function) C P : [0, 1] d → [0, 1]. Let the random d-dimensional vector of inputs X = (X 1 , . . . , X d ) ∼ P , where P ∈ P(X ) and X ⊆ R d . Assume that the marginal cdfs F Pi , i = 1, . . . , d are continuous. Denote U 1 , . . . , U d the random variables defined as: The copula of X, denoted C p is defined as: We refer to the proof of Remark 2 for a proper definition of copula for empirical measures. In the following remark, we show that by means of particular monotone transportation maps inspired by optimal transportation theory, the marginal inputs can be optimally perturbed while preserving their copula. Moreover, this result is suitable for both ML and UQ applications since it also holds whenever P is an empirical measure related to an observed set of data points.
Remark 2. Let X ⊆ R d , for d a positive integer, and P ∈ P(X ). Let Q i be the solution of the optimal projection problem (2) with C = C V , for every marginal distribution P i of P , i = 1, . . . , d, and where V ⊆ ⊗ d j=1 F ← j . Let the random vectors (i) If P is an empirical measure (i.e., X represents a dataset), then X and the perturbed dataset X have the same empirical copula. Moreover, the empirical measure of every perturbed marginal sample X i converges towards Q i , i = 1, . . . , d.
(ii) If P is atomless, and assuming additionally that V is such that every F ← Qi , i = 1, . . . , d is strictly increasing, then the random vectors X and X have the same copula. Moreover, each perturbed marginal X i ∼ Q i .
In other words, applying the perturbation map (2) to the inputs allows for preserving their copula and hence their dependence structure. If only an initial dataset is observed, applying T to every observation results in a perturbed dataset with, for instance, the same Spearman correlation matrix. Moreover, these transportation maps achieve optimality for various univariate transportation costs (see the proof). Moreover, in the UQ framework, sampling from the perturbed inputs is as simple as applying these maps to simulated samples of P , which can naturally benefit from the large literature available on copula. However, it requires that the marginal gqfs {F ← Qj } 1≤j≤d are accessible, which obviously depend on the projection problem (2). As shown in the next sections, and among many other practical and theoretical motivations, the particular choice of the 2-Wasserstein distance as a projection metric allows for characterizing the optimally perturbed marginal probability measures through their quantile functions and thus simplifies their accessibility.

Wasserstein projections
This section motivates the choice of the 2-Wasserstein distance as a projection metric for the perturbation problem (2). This distance is deeply rooted in optimal transportation theory [112] and has been used successfully in many ML and deep learning applications [40,3]. It has also been extensively studied as a tool for guaranteeing distributional robustness to adversarial attacks in ML [32]. In SA, it has been used to produce novel sensitivity indices [38,17].
Using the Wasserstein distance imposes to work on the subset P p (R d ) ⊂ P(R d ) of probability measures with finite p-th moment. The Wasserstein distance between multi-dimensional probability measures can be computationally expensive to evaluate [87]. However, as stated in Section 3.4, the fact that one wishes to preserve the copula between P and its optimally perturbed counterpart Q greatly simplifies the projection problem. Let P, Q ∈ P p (R d ) be two multi-dimensional probability measures, with marginals P 1 , . . . , P d and Q 1 , . . . , Q d respectively. Leveraging the work in [1], if P and Q share the same copula, one can rewrite their 2-Wasserstein distance as: Hence, finding Q that minimizes W p p (P, Q) under constraints on the marginals Q 1 , . . . , Q d is equivalent to solving d independent univariate projection problems with the p-Wasserstein distance as a projection metric. Furthermore, the transportation map defined in (10) is indeed optimal [1] (but not unique), in the sense that it minimizes (11).
In other words, finding Q that minimizes (11) such that C P = C Q is equivalent to projecting each P i under its relevant constraints, and applying the transportation map (10). As the problem reduces to perturbations of univariate marginal distributions, the p−Wasserstein distance may be easily computed, as recalled in the following definition.
Definition 4 (Wasserstein distance on the real line). Let p ∈ N * and P, Q ∈ P p (R) be two probability measures on R admitting F P and F Q as probability distribution functions, respectively. Then, the p-Wasserstein distance between P and Q is: has been defined in Definition 1. The following subsections argue on the specific choice of the 2-Wasserstein distance. First, we highlight its attractive properties for conducting robust interpretability analyses. Then we investigate the solution of the perturbation problem (2), with and without regularity constraints, the latter enforced using isotonic, piece-wise polynomial approximations.

The 2-Wasserstein distance as a suitable perturbation discrepancy
The special choice of the 2-Wasserstein distance to instantiate the perturbation problem (2) is based on the following rationale. First, W 2 metricizes weak convergence on P 2 (R) ( [112], Section 7.2). It means that W 2 is a measure of proximity on a broad set of probability measures. In other words, for any P ∈ P 2 (R), and a sequence of probability measures (Q n ) n∈N * ∈ P 2 (R): where d → denotes the convergence in distribution (or weak convergence). That is, W 2 allows for assessing the point-wise proximity between two probability measures, as long as both admit finite second-order moments (a current assumption in both SA and ML fields). Contrary to the KL divergence, no additional conditions on the absolute continuity of Q n w.r.t. P are needed. When it comes to the perturbation problem (2), two practical advantages in favor of the 2-Wasserstein distance can be drawn, compared to entropic projections (i.e., using the KL divergence): if P is an empirical measure (i.e., purely atomic), then Q is not restricted to be purely atomic; conversely, if P admits a density, then it does not restrict Q to admit a density.
These benefits are key in unifying the frameworks of SA and ML interpretability: the flexibility of W 2 allows for greater explicit control (e.g., through smoothing restriction) on the nature of Q, independently of that of P . Results of entropic projections entail a re-weighting of the atoms when P is empirical [6], or a perturbed density of Q proportional to the one of P when it is absolutely continuous w.r.t. the Lebesgue measure [69]. Using W 2 allows, for instance, to add additional atoms, or to allow Q to admit a density, independently of the regularity of P .
Second, W 2 facilitates the projection of P onto a quantile perturbation class C V (θ). Indeed, the next proposition shows that the optimization problem is equivalent to a projection in L 2 of its gqf onto V under interpolation constraints Proposition 1. Let P ∈ P 2 (R), and C V (θ) be a non-empty perturbation class defined by a subset V ⊆ F ← . Consider for V the constraints system defined in § 3.2-3.3 associated with couples (α i , b i ) 1≤i≤K . In this frame, the solution Q of the perturbation problem (2), i.e., is characterized as the unique Lebesgue-Stieljes measure induced by the cdf F Q with gqf F ← Q ∈ F ← : The equivalent projection problem in (13) echoes with the result in Proposition 1. Projecting a measure w.r.t. the 2-Wasserstein distance is equivalent to projecting its gqf in L 2 . One can then leverage the vast literature on function approximations in L 2 , especially on monotonic approximations. Moreover, as eluded at the end of Section 3, the proposed perturbation scheme depends heavily on the knowledge of the gqf of Q. Solving (13) grants direct access to the gqf of Q, and thus allows applying perturbations fairly easily.

Quantile constrained Wasserstein projections
This subsection presents and discusses the main results of this paper. If no smoothing constraints on F ← Q ∈ F ← is enforced, the perturbation problem (13) leads to a unique analytical solution. However, many studies require F ← Q to be smooth (e.g., continuous). We propose enforce continuity by using isotonic interpolating piece-wise polynomials, which lead to a well-defined convex constrained quadratic program, easily solvable in practice.

Analytical solution without smoothing restrictions
The following proposition provides a convenient way to solve the perturbation problem (13) in the case when no smoothing constraint on F ← Q is enforced. Proposition 2. Let P be a probability measure in P 2 (R). Let C be a non-empty perturbation class characterized by a set of K quantile constraints. Assume, without loss of generality, for i = 1, . . . , K, that α 1 < · · · < α K along with b 1 < · · · < b K . Let β i = F P (b i ) for i = 1, . . . , K. Define the intervals A i = (c i , d i ] for i = 1, . . . , K, such that: Then the problem (13) has a unique solution which can be written as, for any y ∈ [0, 1]: In order to interpret this result, illustrated in Figure 5, let us recall that when a quantile function is constant on an interval, it implies that its related probability measure admits an atom at the constant value taken by the gqf. Moreover, the mass allocated to this atom is equal to the length of the interval. Additionally, each jump of the quantile function induces an interval with no mass. The solution displayed in (14) shows that on A, both initial and perturbed quantile functions are equal. However, they differ on every interval A i in the following fashion: • Q have atoms at each constraint point b i , i = 1, . . . , K; • Each of these atoms have mass • Each open interval I i ⊂ R defined as with, by convention, b 0 = −∞ and b K+1 = ∞, has no mass. To put it briefly, Q(I i ) = 0 for every i = 1, . . . , K.
In other words, whenever an α-quantile p α is shifted up to a value b, the perturbation entails sending every possible values in the range (p α , b) to b. Hence, every value in (p α , b) cannot be sampled according to Q. Moreover, the singleton {b} now admits a probability of being observed equal to the initial probability of this interval, i.e., Q({b}) = P (p α , b) . When an α-quantile is shifted down to b, the interval becomes (b, p α ), and the same reasoning can be done.
The statement of Proposition 2 is rather intuitive. Indeed, the Wasserstein distance quantifies the amount of work needed to transform a probability measure into another one [96]. When using W 2 , the amount of work is quantified using the Euclidian distance, i.e., transporting a point x 0 to x 1 requires (x 0 − x 1 ) 2 units of work. This intrinsic point-wise way of quantifying similarities can be sensed in the previous result: perturbing an α-quantile entails giving the initial mass of an interval adjacent to b to the singleton {b} in order to satisfy the constraint.

Projection solution under smoothness restrictions
The analytical solution provided in Proposition 2 presents a significant drawback: part of the application domain Ω X of the perturbed input receives no mass. From a practical standpoint, it cannot be explored in a robustness study. It is because F ← contains discontinuous functions. Ensuring a smoother solution relies on specifying a smooth perturbation class Q V where V is a set of continuous, non-decreasing functions. To solve the perturbation problem, we can take advantage of its L 2 ([0, 1]) equivalent formulation 13) to project F ← P onto V. However, the main challenge is to ensure that V only contains monotonic functions.
We propose to projection F ← P onto a space of piece-wise continuous polynomials. It implies that the support of Q must be bounded. These bounds are made explicit using extremal quantile constraints (i.e., F ← Q (0) and F ← Q (1) are constrained to take finite values). Formally, we aim to find a piece-wise polynomial of the form under the continuity constraints at each knot of the grid α 1 < · · · < α K Here, each G j ∈ R[x] ≤p , for j = 0, . . . , K where we recall that R[x] ≤p denotes the set of all real polynomials of degree at most equal to p. Let S p denote the space of functions defined by (16). Restricting the solution of the perturbation problem (13) leads to the following optimization problem Hence, the smooth perturbation class is defined by V = F ← ∩ S p . Since the design enforced the polynomials in S p to be defined on the grid α 0 < α 1 < · · · < α K < α K+1 = 1, solving (17) reduces to solve several sub-problems on each sub-interval [α i , α i+1 ], i = 0, . . . , K of [0, 1]. The optimization problem is indeed separable into K + 1 independent optimization sub-problems. Each of them defines an optimal component G i of the piece-wise polynomial G as defined in (16).
Any of these problems can be formulated generically as follows. Let [t 0 , t 1 ] ⊂ [0, 1], and z 0 , z 1 ∈ R be interpolation values at t 0 and t 1 respectively. We aim to find the solution to the optimization sub-problem This optimization sub-problem is nothing more than the L 2 isotonic (i.e., monotonic, in this case non-decreasing) polynomial approximation on a compact interval [75], with interpolation constraints at the boundaries. The interpolating polynomials have been extensively studied in the literature [39], as well as isotonic polynomial regression and approximation [97,115]. However, to our knowledge, this specific optimization problem does not seem to have been particularly studied. We propose to solve (18) using the sum-of-squares (SOS) polynomials [68] representation of non negative polynomials. Furthermore, we leverage in particular the representation of SOS polynomials using semi-definite positive (SDP) matrices [82,83,101]. A similar characterization of isotonic polynomials has been proposed in [101]. The following result, the second main outcome of this article, shows that the problem to solve falls into the category of strictly convex programs: the solution in (21) is unique [15].
and denote r ∈ R d+1 the moment vector of F → P (x), i.e., for i = 0, . . . , d Then, the vector s * = (s 0 , . . . , s d ) ∈ R d+1 of coefficients characterizing the polynomial S in (18) is the solution of the following convex constrained quadratic program where K is an identifiable closed convex subset of R p+1 (for the sake of conciseness, K is characterized within the proof ).
As the computation of s * is a convex constrained quadratic program, it can be addressed efficiently using devoted solvers. The problem (17) can be addressed by solving K + 1 optimization problems of the form (21). Furthermore, computations can be done in parallel, leading to fast computational times. Notice that (21) can be formulated and solved using CVXR. This is an R package for disciplined convex programming [41]. The pretty generic low-level logic behind the optimization scheme can be found in Algorithm 1.
While computing the Lebesgue moment matrix M on each sub-interval of [0, 1] is straightforward, computing strategies for r, the moment vector of F ← P , can vary depending on the nature of P . Additional computational details are given in Appendix Appendix B. The set-up of the CVXR constraints is detailed in the accompanying GitLab repository 1 .
To provide a frame of reference for the practical usage of our method, the empirical computational time of solving one element of G, w.r.t. the polynomial degree is studied as follows. Values t 0 , t 1 ∈ [0, 1], and z 0 , z 1 ∈ Ω X are randomly selected, and an isotonic interpolating piece-wise continuous polynomial is fitted (i.e., solving (21)). Polynomials of degrees ranging from 2 to 50 are fitted for each experiment, repeated 150 times. The execution time has been recorded and is displayed in Figure 6. One can notice that the mean computational time seems to be linear w.r.t. the polynomial degree. However, the higher the degree, the wider the 90% time coverage seems to be, which may be caused by the complexity of the underlying optimization problem. In our limited testing, further numerical experiments showed that small polynomial degrees (≤ 7) often appear sufficient to obtain good approximations and that the approximation error tends to stabilize, w.r.t. the polynomial degree, rather rapidly.
Remark 3. The numerical solver used is SCS V3.2.1 [79]. To improve numerical stability, the quantile functions have been mapped to take values between [−1, 1]. All the figures and all obtained optimal perturbations have been computed by performing this pre-processing step first.

Robustness diagnostics to distributional perturbations
We illustrate the previous method on two use cases. First, in an ML context, we assess the robustness to feature perturbations of a classification model (i.e., a one-layer neural network) trained on an acoustic fire extinguisher dataset. Then, in the UQ framework, we extend Example 2, studying the impact of input perturbation on the output of a numerical hydrological model predicting a river water level.
Remark 4. In the following applications, isotonic polynomial smoothing is applied with an arbitrarily high degree. We chose the degree based on an empirical inspection of the solutions and by ensuring that the approximation error remains relatively the same w.r.t. higher degrees.
Our method is in line with the general SIPA (Sampling, Intervention, Prediction, Aggregation) framework for model-agnostic interpretation proposed in [98]. More precisely, the four step can be broken down as follows: 1. Sampling: In Section 5.1, we have access to an i.i.d. sample of the inputs. In Section 5.2, a sampling strategy is used to simulate data. 2. Intervention: In both use-cases, we define meaningful quantile perturbations, and solve the subsequent perturbation problem. Then we apply the transformation in (10), resulting in perturbed samples.

Prediction:
We predict using the available black-box model (a neural network in 5.1, and a numerical model in Section 5.2), resulting in perturbed outputs. 4. Aggregation: The resulting perturbed outputs are aggregated w.r.t. the perturbation intensity, leading to global robustness metrics, or are simply plotted against the initial and perturbed of the perturbed input values, allowing for local robustness assessments.

ML application: Acoustic fire extinguisher dataset
The acoustic fire extinguisher dataset is composed of 15390 experiments of fire extinguishing tests of three different liquid fire fuels. Amplified subwoofers are placed in a collimator with an opening.  When activated at different frequencies, the acoustic waves produce an escape of air through the opening, which is used to extinguish fires. Three features are set using a design of experiment (DoE), and two are measured using appropriate equipment. For more details on the experiments settings, one can refer to the in-depth descriptions in [65,107].  For each experiment, a binary output variable Y is measured, representing the result of the experiment, i.e., whether the fire has been put out (Y = 1) or not (Y = 0). The two output classes are relatively balanced (i.e., 48.97% of the observations describe effectively put out fires). The distribution, correlation structure, and relationship of the continuous features with the output are represented in Figure 7. Some variables seem fairly correlated (in Spearman's sense, i.e., the linear correlation of the rank-transformed data), such as Frequency and Decibel, as well as Distance and Airflow.
The classification black-box model is a one-layer neural network (composed of 100 neurons), trained on 500 epochs, with a learning rate of 10 −4 , similar to the study conducted in [106]. 5% of the data has been randomly selected to serve as validation data. The model resulted in a good prediction accuracy: 95.15% of the training data and 94.26% of the validation data are correctly classified. Figure 8 depicts the ROC curve and confusion matrix of the trained black-box model. The model's predictive performance can be validated globally with an AUC of 0.992 and less than 3% of type 1 and 2 prediction errors.
However, global predictive performance only focuses on effectively observed data points. It is mandatory to study the model's behavior on predictions outside of these points to improve confidence in its usage. Hence, one can be interested in the robustness of the model w.r.t. perturbations on its inputs. Note that ground truths cannot be observed for perturbed data. However, the impact, either globally or locally, of these perturbations on the predictive behavior of the model can still be assessed using predictions on the perturbed data. In the following, the feature perturbation scheme is detailed and motivated, and then the model's behavior is studied under these perturbations.

Perturbation strategy
We propose a straightforward perturbation strategy. Only the Airflow feature is perturbed. The perturbation is composed of the K = 14 constraints: • The application domain of the feature is preserved by setting both the 0 and 1-quantiles to the minimum and the maximum observed value of the dataset.
• The left tail of the distribution is preserved by constraining every quantile of level 0.1 to 0.6 with a step of 0.05 to interpolate the empirical quantile function of the feature.
• A quantile shift perturbation is put on the 0.8-quantile of the feature, with an initial value of F ← P (0.8) = 12, being shifted between 9.5 (θ = −1) and 14.5 (θ = 1). Additionally to these perturbations, smoothness is enforced using piece-wise continuous isotonic polynomials, as described in Section 4.2.2. The degree of each monotone polynomial has been arbitrarily chosen to be up to 9. The constraints and the resulting quantile-constrained Wasserstein projections are illustrated in Figure 9 for intensity values −1, 0, and 1.
The perturbed quantile level has been chosen in relation to the model's decision boundary: no observation in the initial dataset with an Airflow value exceeding 12.3m/s is classified by the model as not extinguishing the fire, regardless of the values taken by the other features. Perturbing the 0.8quantile of the Airflow variable allows for exploring the model's behavior in regions close to this decision boundary. More importantly, it allows assessing the predictive robustness of the neural network in this region under perturbations of varying magnitude. Generally, this quantile shift regime can be understood as a perturbation on the right tail of the initial distribution, i.e., on values higher than the 0.6-quantile.

Model robustness assessment
First, we are interested in assessing the robustness of the neural network model in a global fashion. The left plot of Figure 10 presents the proportion of perturbed observations with predictions of 1    w.r.t. to the intensity of the perturbation. Notice that the proportion is increasing, along with θ. Hence, decreasing the value of the initial 0.8-quantile tend to result in a lower number of predicted put-out fires, and increasing its value results in an increasing number of predicted put-out fires. This interpretation is rather intuitive: all other things being equal, a higher Airflow value entails a higher chance of predicting Y = 1. The right plot of Figure 10 presents the proportion of prediction shift with respect to θ. Notice that the higher the magnitude of the perturbation (either positively or negatively), the more predictions tend to change, and the closest θ is to 0, the fewer predictions shift. This observation informs on the predictive stability in the vicinity of the decision boundary of the model: small perturbations tend to result in less prediction shift than bigger perturbations. Second, we study the robustness of global SA results. Figure 11 presents the target Shapley effects [59], a global SA input importance measure for binary black-box model outputs with dependent inputs, w.r.t. the perturbation intensity parameter θ. These indices have been computed using the nearestneighbor (KNN) approach proposed in [18] (with an arbitrarily chosen number of neighbors equal to 6). Recall that our perturbation method allows preserving the empirical copula between the features, justifying the use of the KNN approach. Studying the behavior of importance measures allows to verify if the feature importance order shifts due to the perturbations, i.e., if the importance hierarchy between the inputs changes due to perturbations around the model's decision boundary. The left barplot presents the initial target Shapley effects, computed on the model's prediction on the observed data, and the right plot presents their behavior under the airflow perturbation. One can notice that the importance indices remain stable w.r.t. θ. This result indicates that the global SA of the neural network is robust to the distributional perturbations driven by θ. Hence, we can be confident in the importance measures under uncertainties in the region near the model's decision boundary.
Finally, the robustness of the neural network can also be assessed locally. Figure 12 allow visualizing whether a prediction has shifted w.r.t. to the effective magnitude of the perturbation. The black line indicates no perturbation change: the airflow value of an observation has been mapped to itself. For a fixed initial airflow datapoint, its vertical distance to the black line indicates the (signed) magnitude of the applied perturbation. Red points indicate that the prediction has shifted w.r.t. the initial dataset, and blue points indicate no predictive change. One can note the presence of red dots close to the black line around the prediction boundary of the model. Small perturbations for observations with airflow values around 12, all other features being equal, can lead to a prediction change. Hence, the confidence in predictions on observations in this region can be questioned. However, notice the lack of red dots near the black line for airflow values on the interval [13,17] and on the interval [7,10]. Hence, we can be confident in the model's predictions for Airflow values on these intervals, which seem to be robust w.r.t. the quantile shift. One may notice the presence of small perturbation resulting in prediction changes for small airflow values. However, since the perturbation scheme focuses on exploring the model's behavior around its Airflow decision boundary, their interpretation is voluntarily omitted: a different perturbation scheme involving perturbing the left tail of the airflow distribution would be advised.
In summary, besides the model's good prediction accuracy, it also seems globally robust to distributional perturbation focused around its decision boundary. Moreover, we can be confident in the feature importance indices since they remain relatively similar under perturbation. Locally, the model prediction seems stable w.r.t. small perturbations, except on a small interval around its decision boundary (a behavior generally expected in ML applications). In conclusion, this robust interpretability analysis further assesses the model's behavior beyond classical accuracy metrics and provides additional arguments for its validation.

SA application: Simplified hydrological model
This use-case focuses on a simplified model of the water level of a river. This model has been extensively used in the safety and reliability of industrial sites, where the occurrence of a flood can lead to dramatic human and ecological consequences. It consists of a substantial simplification of the one-dimensional Saint-Venant equation, with a uniform and constant flow rate, inspired from [62,42]. The maximal annual water level from sea level is modeled as: where the description of each input variable and their explicit marginal probabilistic structure is detailed in Table 2.
Additionally, similarly to [21], a dependence structure is modeled using a Gaussian copula, with the covariance matrix Echoing Example 2, we are interested in uncertainties on the application domain of the K s input, i.e., the Strickler riverbed roughness coefficient (which is the inverse of the Manning coefficient). Its value can range from around 3 (proliferating algae) to around 90 (smooth concrete). We refer the interested reader to the in-depth study in [42] for more details on the determination and inference of the Strickler coefficient for realistic rivers. In this use-case, initially, the application domain Ω X of the Strickler coefficient is set between the values of 20 and 50, corresponding to situations from   very cluttered riverbeds to earthen channels. However, to illustrate our robustness method, epistemic uncertainties are assumed to affect this application domain.

Perturbation strategy
In this use case, the three following inputs are perturbed. The river maximum annual water flow rate Q, the river length L, and the upstream river level Z m are subject to the following quantile constraints: • Quantile perturbations on Q: The initial input distributions, their application domain, and the optimally perturbed results are illustrated in Figure 13. These constraints are mainly enforced to illustrate that multiple inputs can be perturbed simultaneously while preserving their dependence structure.
In addition to these constraints, the Strickler coefficient K s is subject to an application domain dilatation perturbation, with a scaling parameter η = 2. Each perturbation intensity represents a degree of uncertainty on the type of riverbed roughness. When θ = −1, the width of the initial application domain is halved, i.e., from [20,50] to [27.5, 42.5], which can be interpreted in a situation where the epistemic uncertainty on the riverbed roughness is narrower, between a slow winding natural river, up to a plain river without shrub vegetation. When θ = 1, the epistemic uncertainty on the riverbed is much wider, with an application domain equal to [5,65], which depicts a range of riverbed roughness from proliferating algae up to smooth concrete. Figure 14 illustrates the initial K s distribution, along with the optimally perturbed quantile functions for θ being equal to −1 and 1.
Additionally, the perturbations' smoothness is enforced using piece-wise continuous isotonic polynomials of degree up to 12, chosen arbitrarily.

Robustness of the sensitivity analysis
From a global standpoint, one can be interested in the impact of the distributional perturbations on key statistics of the random output of the river water level model. Figure 15 presents estimated values for the mean, standard deviation, 0.025 and 0.975-quantiles (shown by the 95% coverage), and minimum and maximum values of the random output, computed on 10 5 Monte Carlo samples, w.r.t. the dilatation intensity θ. These values are compared to the reference ones according to the initial distribution of the inputs, estimated on a 2 × 10 5 Monte Carlo sample.
Notice that the expectation, standard deviation, 95% coverage quantiles, and minimum value of the model output remain stable under the distributional perturbations on the application domain of the Strickler coefficient. However, the estimated upper bound of the output support increases exponentially for positive values of θ. Widening the uncertainty on the type of riverbed allows for relatively rare events of high river water levels since the 0.975-quantile does not seem to be dramatically affected by the distributional perturbations.  Figure 16 presents the Shapley effects [80], which are global SA importance measure for real-valued model outputs with dependent inputs. These indices have been computed using a double Monte Carlo scheme as depicted in [102], with fixed simulated sample sizes, for each perturbed distribution Q driven by a value of θ, N v = 10 4 for estimating Var Q (Y ), as well as N o = 10 3 and N i = 100 to estimate E Q [Var Q (Y | X A )] for every subset X A , A ⊆ {1, . . . , d} of variables. Additionally, the reference Shapley effects have been computed under the initial distribution with sample sizes N v = 10 5 , N o = 3 × 10 3 and N i = 300.
Note that the distributional perturbations have an impact on the importance measures. More precisely, increasing the range of the uncertainty of the riverbed roughness increases its importance for positive values of θ. Conversely, the importance of both Q and Z v decreases accordingly. However, the variable importance hierarchy induced by the Shapley effects is preserved. It is also essential to notice that both Q and Z v tend to be considered equally important as θ gets large. Hence, this SA does not seem robust to distributional perturbations and, more precisely, to a widening of the support of the Strickler coefficient in combination with the quantile perturbations put on Q, L, and Z m .

Conclusions
These two use cases illustrate the usefulness of our method in both UQ and ML studies. On the ML side, for classification tasks, it allows assessing the global behavior of black-box models under input perturbations. This assessment is quantified either through studying the prediction shifts due to the perturbation, or through the behavior of feature importance metrics. Locally, it allows the detection of low-stability regions of interest (regions where small perturbations induce a classification change).
Overall, in addition to classical accuracy metrics, our method can be used to assess confidence in a predictive model. On the UQ side, it allows for studying the impact of distributional perturbations (whose  intensity can be tuned to represent epistemic uncertainties) on the model output, even in situations where inputs are correlated. Furthermore, in a SA context, the behavior of classical sensitivity indices under those perturbations can also be studied, and their robustness (for instance, the preservation of the input importance hierarchy) w.r.t. the probabilistic modeling on the inputs can be assessed.

Discussion and perspective
Obtaining a robustness diagnosis on the influence of input variables and the behavior of a model considered a black box is essential for its acceptance and use. Such models can either implement solutions of mechanistic equations or be learned from data. In both cases, the sensitivity of key indicators to a misspecification of the probabilistic input model must be evaluated. It is essential to have specific intelligible rules and well-defined computational tools to achieve this goal. This paper provides an answer to this question by proposing to modify the distributions of the input variables, seen as features. These perturbations modify the quantile of marginal distributions while preserving the dependence structure. We project the initial distribution under a Wasserstein cost to design optimal perturbations. The proposed approach allows the inclusion of regularity conditions through piece-wise isotonic polynomials. The robustness analyses conducted on real case studies illustrate its potential flexibility and speed, which are essential for scaling up (increasing dimension and size of data). These methods and examples have been implemented within a dedicated openly accessible GitLab repository.
We will continue this work by examining the following points, which offer specific avenues for further research. The first technical problem is determining rules for the degree of the nonnegative polynomials in smooth approximations. It can be understood as the injection of prior information on the order of differentiability of the sought-after perturbed gqf. If F ← P is supposedly regular (on all or some parts of [0, 1]) and its differentiability can be estimated, a rule-of-thumb heuristic could be to impose the same differentiability on the resulting perturbed gqf F ← Q . In an ML framework, nonparametric approaches to isotonic regression of the marginal gqfs of P can provide answers through statistical testing [34,29] or criteria enforcing a trade-off between approximation error and sparsity (e.g., inspired from AIC or BIC).
The second problem to be addressed is providing rules for modifying the dependence structure between features. In an UQ context, classes of parametric high-dimensional (vine) copulas have been recently explored by [108,14,88,8] for several tasks, including SA. In the light of the dedicated literature, association and concordance measures appear as the most interpretable tools for these multivariate distributions (and therefore frequently used to incorporate expert opinion) [24,117,14]. Relevant classes of parametric vine copulas could be defined objectively by the Fréchet-Hoeffding bounds on the conditional bivariate dependence structures of such copulas or by computing extreme dependence structures relevant for the phenomenon of interest, as proposed by [14]. The latter authors use a sparsity hypothesis to avoid prohibitive computational costs related to the curse of dimensionality. It echoes the strong need for sparsity assumptions about feature dependence in an ML framework, to which preliminary dimensionality reduction steps usually respond. In this context, many ML methods pay particular attention to the search for decorrelated, if not independent, features or small groups of features. It ensures the good behavior of the learning models and the relevance of the indicators allowing their interpretability [46]. For instance, it is well known that permutation-based indicators for random forests provide relevant information on the importance of input features solely when they are independent [48]. New solutions have recently been developed to overcome this limitation [12]. In addition to computational limitations, it seems difficult, if not impossible, to produce understandable indicators when dealing with a large number of features.
An alternative approach to account for feature dependence could be based on multivariate quantile extensions. Among the many approaches to defining such a notion, the most theoretically accomplished today is the one resulting from the concept of center-outward distribution function, based on optimal transportation ideas. It was first proposed and studied by [22] and [52]. The resulting center-outward quantile function offers desirable continuity and invertibility properties, which ensures the existence of closed and nested quantile contours [37]. Suggesting that characteristics of this multivariate quantile function can be used to define variation classes seems natural to generalize our approach to multidimensional perturbations instead of focusing on marginal perturbations. Therefore we suggest that studying the interpretability of these features and the computational work required to handle quantile contours can be useful to improve our approach.
On the subject of isotonic polynomials, one can notice that the solution of the proposed approximation scheme can not be differentiable at the points of the grid. To circumvent this drawback and enforce the derivability of the smoothed gqf on [0, 1], one could enforce additional constraints on the derivate of the polynomials. It would result in smoother perturbed gqf. Doing so would bear resemblance with splines [54], and in particular isotonic splines, which have been extensively studied in the literature [97,39,115].
As stated in the proof of Theorem 1, we chose to use the SOS representation of nonnegative polynomials and their subsequent characterization using semi-definite positive matrices. Works in [31,75] could allow for less convoluted representations, leading to faster solving times. However, our developments allow for more apparent extensions to the definition of multivariate nonnegative polynomials. This decision has been made to account for using multivariate polynomials for smooth perturbation classes of quantile contours, which generalizes our approach to multidimensional perturbations.
Other spaces of functions can also be used for smoothing purposes. Following the work of [7], abstract reproducing kernel Hilbert space of nonnegative functions can be reached through particular kernels. Hence, it would allow accessing different sets of nonnegative functions, whose regularities can be assessed through a thorough study of these kernels.
Finally, one of the primary motivations for using the 2-Wasserstein distance as a projection metric is that it metricizes weak convergence on a broad set of probability measures. Other distances between probability measures are endowed with similar properties, such as the Prokhorov-Levy distance. An exciting improvement of our method would be to leverage the different relationships between such distances (see [45]) to assess their use's relevancy for robustness studies.
Proof of Lemma 1. Notice that if (8) is respected, then the constraints are non-decreasing. Then, there exists at least a function F ← in F ← such that the constraints are respected (e.g., the linear interpolant of the constraints). So that, using Proposition 1, there exists a probability measure with F ← as a generalized quantile function.
Equivalently, one can express b in terms of θ, which directly provides the expression of b α (η, θ).
Proof of Lemma 3. Preserving the midpoint of Ω X while perturbing its width requires that, for any couple Using the transformation allows to define the formulas for b 0 and b 1 provided in the lemma statement. The perturbation intensity θ can be interpreted as follows: • If −1 ≤ θ < 0, the application domain is narrowed, its width being divided by up to η; • If θ = 0, the application domain boundaries are not perturbed; • If 0 < θ ≤ 1, the application domain is widened, its width being multiplied by up to η. where R j,k denotes the rank of x j,k in {x j,i } 1≤i≤n . Perturbing only the marginals of P when solving the optimal projection problem (2) amounts apply to each marginal sample {x j,i } 1≤i≤n the transportation maps and consequently the marginal empirical probability measure of the perturbed samples are defined as wherex j,i = T j (x j,i ) are the perturbed datapoints. Since the marginal cdf and gqf in (A.2) are nondecreasing (ie., monotonic) functions, their composition is non-decreasing as well. Since monotonic functions preserve orders, the ranks of x j,i andx j,i are invariant, for (i, j) ∈ {1, . . . , n} × {1, . . . , d}. From (A.1), it implies that the empirical copula of Q isĈ Q =Ĉ P . This traduces, for instance, the invariance of Spearman correlation matrices between the initial and perturbed datasets. The transportation map T j defined by (10) is not the unique non-decreasing transportation map between P and Q ( [96], Chap. 2). However, T j is indeed an optimal perturbation plan (in Monge's sense) for a wide variety of transportation costs between the marginals P j and Q j ( [87], Remark 2.30). Moreover notice that for any j = 1, . . . , d (ii) Now, for the sake of conciseness, let P denote any univariate probability measure and Q its optimally perturbed counterpart. If F ← Q is strictly increasing then from [33], for all u ∈ [0, 1] Now denote the transportation map Let X ∼ P and define U P = F P (X). Then, one has, if P is atomless, that T (X) ∼ Q and Now, for a multivariate probability measure P with marginals P 1 , . . . , P d , and respective optimally perturbed probability measures Q 1 , . . . , Q d , this entails that, for i = 1, . . . , d: and hence, the random vectors U Q and U P are equal almost surely. Subsequently, for any u ∈ [0, 1] d , C P (u) = C Q (u). Hence the X and T (X) have the same copula.
Proof of Proposition 1. First, one can notice that, for any probability measure G ∈ P(R): where µ denotes the Lebesgue measure on [0, 1] (see, [30]). This entails that: by a continuous composition of µ-a.e. equal functions, and by identification with the definition of the W 2 distance for probability measures supported on R.
In the following of this proof, one refers to the projection problem (12) as the "Wasserstein projection", and the projection problem (13) as the "L 2 projection".
The proposition can be proven in two steps. First, if Q is the solution of the Wasserstein projection, then its quantile function F ← Q is necessarily the solution of the L 2 projection. This is due to the fact that F ← Q is uniquely characterized by Q, as well as the particular case of the 2-Wasserstein distance for measures supported on R (see Definition 4).
Second, let Q be any probability measure in P 2 (R), and denote F ← Q its characterizing quantile function. Assume that F ← Q is the solution of the L 2 projection. Since F ← Q characterizes uniquely Q, thanks to Proposition 1, one has that Q is necessarily the solution of the Wasserstein projection.
The integral can be decomposed as follows: Since the quantile constraints are of the form: we can always write L(y) = b i + h(y) for y ∈ A i , and where h is an non-decreasing, left-continuous function. Moreover, note that: • h(y) is non-negative, and F → P (y) − b i ≤ 0 if A i falls in cases 2. and 4. • h(y) is non-positive, and F → P (y) − b i ≥ 0 if A i falls in cases 1. and 3. Then we have: since h(x) and F → P (x) − b i have different sign. Due to the constraint and the left-continuous nondecreasing nature of L, this bound is tight and is attained if and only if h(y) = 0 for all y ∈ A i . Globally, this entails that and this tight bound is uniquely attained by the left-continuous non-decreasing function defined as Proof of Theorem 1 (ingredients). This proof relies on the following results that can be found in [82,83,68], and further recalled in [101]. They involve sum-of-squares (SOS) polynomials, which can be defined as follows.
Definition 5 (SOS polynomials). A polynomial S of even degree p is said to be a SOS polynomial if, for m ∈ N * , there exists s 1 , . . . , s m polynomials of degree at most equal to d 2 , and such that, ∀x ∈ R: Theorem 3. Let t 0 , t 1 ∈ R such that t 0 < t 1 , and let p ∈ N * .
(i) A univariate polynomial S of even degree d = 2p is non-negative on [t 0 , t 1 ] if and only if it can be written as, ∀x ∈ [t 0 , t 1 ] where Z is a SOS polynomial of degree at most equal to d, and W is an SOS polynomial of degree at most equal to d − 2.
(i) An univariate polynomial S of odd degree d = 2p + 1 is non-negative on [t 0 , t 1 ] if and only if it can be written as, ∀x ∈ [t 0 , t 1 ] where Z, W are SOS polynomials of degree at most equal to d.
It is important to note that Theorem 3 is quite general, in the sense that it allows for extensions to multivariate polynomials (i.e., polynomials taking values from R d ). As pointed out in [31] (Thm. 1.4.2), nonnegative polynomials on compact intervals can also be defined as a linear combination of squared polynomials. It may facilitate the identification of the nonnegative polynomials' coefficients, as done in [75] in the context of statistical learning. However, for the sake of potential future genericity, we chose to leverage the direct powerful link between SOS polynomials and semi definite positive matrices, as expressed in the following theorem.
Theorem 5. Let S n the subspace of real-valued symmetric matrices, in the vector space of square matrices. The set of symmetric SDP matrices Σ N is a proper cone in S n , and thus is a closed convex set.
A few results on the preservation of convexity of sets under transformations are also required. These lemmas can be found in [15].
Lemma 4 (Linear maps preserve convexity). Let V, W be two vector spaces over the same field F . Let T : V → W be a linear map, and let C ⊂ V be a convex set. Then the image of C under T , i.e., : is also a convex set.
Lemma 5 (Cartesian product of convex sets is a convex set). Let C 1 be a subset of R m and C 2 be a convex subset of R n . Then, the Cartesian product C 1 × C 2 is a convex subset of R m × R n .
Two additional results, proven beneath, are required before proceeding to the proof of Theorem 1.  Hence T is a linear map between S p and R 2p .
Lemma 7. Let S be a univariate polynomial of degree d and s = (s 0 , . . . , s d ) ∈ R d+1 its coefficients. Let S be its derivative, i.e., a polynomial of degree d − 1, with coefficientss = (s 1 , . . . , s d ) ∈ R d . Let Z and W be SOS polynomials, with coefficients z and w, and assume that S is non-negative on [t 0 , t 1 ] as a combination of Z and W as in Theorem 3. Moreover, let D = diag(1, 2, . . . , d) be the (d × d) diagonal matrix with (1, . . . , d) as a diagonal elements and denote the bloc-matrices where 0 i,d denotes the (i × d) matrix of zeros, and I d be the (d × d) identity matrix. If d is odd, then z ∈ R d and w ∈ R d−2 and furthermores = Az + Bw where A and B are (d × d) and (d × d − 2) matrices, respectively. If the degree d of S is even, one has that z, w ∈ R d−1 and furthermore:s = Cz + Dw.
where C and D are (d × d − 1) matrices. More specifically, Proof of Lemma 7. First, assume that S is a polynomial of odd degree d = 2p + 1, meaning that its derivative, S , is a polynomial of even degree 2p. From Theorem 3, one has that S (x) is positive on an interval [t 0 , t 1 ] if and only if it can be expressed as : where Z is an SOS polynomial of degree at most equal to d − 1 and W is an SOS polynomial of degree at most equal to d − 3. Denotes = (s 1 , . . . , s d ) ∈ R d the coefficients of S and z = (z 1 , . . . , z d ) ∈ R d and w = (w 1 , . . . , w d−2) ∈ R d−2 the coefficients of Z and W respectively. One has that : (j + 1)s i+1 x i and if S is assumed to be non-negative on [t 0 , t 1 ] w j+1 x j leading to the following identification : or, written in a matrix form: If S is assumed to be a polynomial of even degree d = 2p, S is necessarily odd degree. From Theorem 3, one has that S (x) is positive on an interval [t 0 , t 1 ] if and only if it can be expressed as : where Z and W are SOS polynomials of degree at most equal to d − 2 with z = (z 1 , . . . , z d−1 ) ∈ R d−1 and w = (w 1 , . . . , w d−1) ∈ R d−1 as coefficients, respectively. It leads to the following identification: We can now proceed to prove Theorem 1.
Proof of Theorem 1 (rationale). This rationale can be broken down in two steps: (a) proving that the objective function (18) can indeed be written in a quadratic form, and:(b) proving that the problem constraints form a feasible set in R d+1 which is closed and convex.
(a) Notice first that the initial objective function where L ∈ R[x] ≤d with coefficients s ∈ R d+1 , can be rewritten as: Note that and further notice that M is thus positive definite since, for any u ∈ R d+1 , u M u = x i F → P (x)dx = s r where r ∈ R d+1 is the moment vector of G with respect to the Lebesgue measure on [t 0 , t 1 ], defined for i = 0, . . . , d as: Since a polynomial is completely characterized by its coefficients, searching for:  can be written as where, for a ∈ R, one denote a d the vector of powers of a up to d, i.e., a d = (1, a, . . . , a d−1 , a d ) ∈ R d+1 . Moreover, by letting: where T is a (2 × d + 1) bloc-matrix, the constraint can be written as: Furthermore, notice that  Moreover, denote the following sets: and notice the polynomial Z (resp. W ) is SOS if and only its coefficients z (resp. w) are in Z (resp. W) thanks to Theorem 5. In addition again, notice that, thanks to Lemma 4, and due to the fact that Σ p is a closed convex set in S p as per Theorem 5, both Z and W are convex subsets of R 2p and R 2q respectively. Besides, thanks to Lemma 5, the set Z × W is a convex subset of R 2p × R 2q as well.
Moreover, let and note that it is a convex subset of R d+1 due to the fact that T 0 is a linear map. Finally, since both C 0 and C 1 are convex sets, their intersection: is as well, and note that any element s ∈ K are the coefficients of a polynomial respecting both equality and monotonicity constraints. In other words, K is the feasible set of coefficients of the initial optimization problem.