Jared Miller
email: miller.jare@northeastern.edu

Matteo Tacchi
email: matteo.tacchi@gipsa-lab.fr

Mario Sznaier
email: msznaier@coe.neu.edu

Ashkan Jasour
email: jasour@jpl.caltech.edu

Peak Value-at-Risk Estimation for Stochastic Differential Equations using Occupation Measures

This paper proposes an algorithm to upper-bound maximal quantile statistics of a state function over the course of a Stochastic Differential Equation (SDE) system execution. This chance-peak problem is posed as a nonconvex program aiming to maximize the Value-at-Risk (VaR) of a state function along SDE state distributions. The VaR problem is upper-bounded by an infinite-dimensional Second-Order Cone Program in occupation measures through the use of one-sided Cantelli or Vysochanskii-Petunin inequalities. These upper bounds on the true quantile statistics may be approximated from above by a sequence of Semidefinite Programs in increasing size using the moment-Sum-of-Squares hierarchy when all data is polynomial. Effectiveness of this approach is demonstrated on example stochastic polynomial dynamical systems.

Introduction

This paper analyzes maximal (1 -ϵ)-quantile statistics of a state function p(x) for Stochastic Differential Equation (SDE) trajectories evolving in a compact set X. An example of this type of quantile statistic for trajectory analysis is in establishing that there exists at least one time with a 1% chance of the aircraft exceeding a height of 100 meters. This task of quantile estimation is related to peak and Value-at-Risk (VaR) estimation, and will also be referred to as the 'chance-peak' problem.

The ϵ-VaR is the value at which there is an ϵ-probability of exceedance [START_REF] Jorion | Value at Risk[END_REF]. Control and portfolio design typically aims to minimise the VaR. One specific VaR-upper-bounding coherent risk measure [START_REF] Artzner | Coherent measures of risk[END_REF] that results in convex programs is the conditional VaR risk measure [START_REF] Pflug | Some Remarks on the Value-at-Risk and the Conditional Value-at-Risk[END_REF][START_REF] Rockafellar | Conditional Value-at-Risk for General Loss Distributions[END_REF]. The conditional VaR has been utilized for stochastic optimal control in [START_REF] Miller | Optimal control of conditional value-at-risk in continuous time[END_REF], and for approximation of discrete-time riskbounded sets using exponential and logarithmic inequalities with Markov Decision Processes in [START_REF] Chapman | Risksensitive safety analysis using conditional value-at-risk[END_REF]. In contrast, the chance-peak approach upper-bounds maximum VaR of the continuous-time SDE state distribution of x(t) across all times. We will solve this problem by maximizing the Cantelli [START_REF] Cantelli | Sui confini della probabilita[END_REF] and Vysochanskij-Petunin (VP) [START_REF] Vysochanskij | Justification of the 3σ rule for unimodal distributions[END_REF] upper bounds for the VaR.

Chance constraints are an adjacent topic to VaR optimization, in which a probability inequality must hold as a hard constraint. Chance-constrained programs have a wide variety of application in control theory [START_REF] Charnes | Cost Horizons and Certainty Equivalents: An Approach to Stochastic Programming of Heating Oil[END_REF][START_REF] Lagoa | Probabilistically Constrained Linear Programs and Risk-Adjusted Controller Design[END_REF][START_REF] Ben-Tal | Robust Optimization[END_REF], and are generally intractable to solve explicitly. Approximation methods for chance constraints include the Cantelli [START_REF] Cantelli | Sui confini della probabilita[END_REF] and VP [START_REF] Vysochanskij | Justification of the 3σ rule for unimodal distributions[END_REF] inequalities, and application of these tailbounds in control include [START_REF] Wang | Non-gaussian chance-constrained trajectory planning for autonomous vehicles under agent uncertainty[END_REF][START_REF] Han | Non-Gaussian Risk Bounded Trajectory Optimization for Stochastic Nonlinear Systems in Uncertain Environments[END_REF]. The scenario approach for randomized constraint generation will converge in probability to the chance-constrained optimum, but carries a risk of failure and may require a large number of samples [START_REF] Calafiore | Uncertain convex programs: randomized solutions and confidence levels[END_REF]. The moment-Sum of Squares (SOS) hierarchy of Semidefinite Programs (SDPs) will converge to the chance-constrained optimal solution under appropriate boundedness conditions [START_REF] Jasour | Semidefinite Programming For Chance Constrained Optimization Over Semialgebraic Sets[END_REF].

The chance-peak problem is also related to a family of optimal stopping problems which can be solved using occupation measures. The work in [START_REF] Lewis | Relaxation of Optimal Control Problems to Equivalent Convex Programs[END_REF] expressed optimal control problems of Ordinary Differential Equations (ODEs) as an infinite-dimensional Linear Program (LP) in an initial, terminal, and occupation measure. The peak estimation problem to maximize a state function p(x) is an instance of optimal control with free terminal time and zero running cost. The work in [START_REF] Cho | Linear Programming Formulation for Optimal Stopping Problems[END_REF] generalizes this LP to the stochastic case to find the maximum expectation of p(x) when dynamics are phrased in terms of their infinitesimal generator (Feller process). Such LPs will converge to the true solution of the stopping problem under mild convergence, regularity, and well-posedness assumptions. The moment-SOS hierarchy of finite-dimensional SDPs will converge to the infinite-dimensional LP optimum if all problem data (e.g., dynamics, constraint sets) are polynomial-representable [START_REF] Lasserre | Moments, Positive Polynomials And Their Applications[END_REF]. This convergent SDP approach has been used for optimal control [START_REF] Henrion | Moment-SOS hierarchy and exit time of stochastic processes[END_REF], peak estimation [START_REF] Fantuzzi | Bounding Extreme Events in Nonlinear Dynamics Using Convex Optimization[END_REF] including compact-valued uncertainty [START_REF] Miller | Peak Estimation for Uncertain and Switched Systems[END_REF], expectation-maximization of Lévy processes [START_REF] Kashima | An optimization approach to weak approximation of Lévy-driven stochastic differential equations[END_REF], and option pricing [START_REF] Lasserre | Pricing a class of exotic options via moments and SDP relaxations[END_REF]. Other instances of the moment-SOS hierarchy used to solve stochastic safety problems include Barrier certificates [START_REF] Prajna | A Framework for Worst-Case and Stochastic Safety Verification Using Barrier Certificates[END_REF], infinite-time averages [START_REF] Fantuzzi | Bounds for Deterministic and Stochastic Dynamical Systems using Sum-of-Squares Optimization[END_REF], and Reach-Avoid sets [START_REF] Xue | Reach-Avoid Analysis for Stochastic Differential Equations[END_REF].

The contributions of this paper are:

• An infinite-dimensional Second-Order Cone (SOCP) that upper-bounds the chance-peak program using VaR inequalities [START_REF] Cantelli | Sui confini della probabilita[END_REF][START_REF] Vysochanskij | Justification of the 3σ rule for unimodal distributions[END_REF] • A convergent set of SDPs using the Moment-SOS hierarchy to the SOCP upper bound

• A verification of this approach on example polynomial SDE systems This paper has the following structure: Section 2 gives an overview of notation, SDEs, and occupation measures. Section 3 upper-bounds the chance-peak problem using an infinite-dimensional SOCP in occupation measures. Section 4 reviews the moment-SOS hierarchy and presents a hierarchy of SDPs that approximate the infinte-dimensional chance-peak SOCP. Section 5 provides numerical examples of the chance-peak problem on ODE and SDE systems. Section 6 concludes the paper.

Preliminaries

Notation

The n-dimensional real Euclidean space is R n . The set of natural numbers is N, and its subset of natural numbers between 1 and N is 1..N . The set of n-multi-indices is N n . The degree of a multi-

index α ∈ N n is |α| = n i=1 α i . A monomial is a term x α = n i=1 x αi i with degree deg x α = |α|. A polynomial p(x) ∈ R[x]
may be uniquely represented in terms of multi-indices α and coefficients p α as p(x) = α∈J p α x α for some finite set J ∈ N n . The degree of a vector of polynomials (f ∈ (R[x]) N ) is the maximum degree of any coordinate (deg f = max i∈1..N deg f i ). The vector space of polynomials of degree at most d is R[x] ≤d and its dimension is n+d d . The Second-Order Cone (SOC) (or Lorentz cone) is

Q n = {(s, κ) ∈ R n ×R + : κ ≥ ∥s∥ 2 }, where ∥s∥ 2 = (s 2 1 +. . .+s 2 n ) 1/2 denotes the Euclidean norm.
The vector space of continuous functions over a topological space X is C(X), and its nonnegative subcone is C + (X). The topological dual of a Banach space V is V * . The cone of (nonnegative) Borel measures supported over X is M + (X) and the vector space of signed Borel measures supported on X is M(X) = M + (X) -M + (X). When X is compact, C(X) and M(X) are topological dual spaces that have a duality product by Lebesgue integration: for f ∈ C(X), µ ∈ M(X) the duality product is ⟨f, µ⟩ = X f (x)dµ(x). This duality product also induces a duality pairing between C + (X) and M + (X). As a slight abuse of notation, we extend this duality product to all Borel measurable functions f : ⟨f, µ⟩ = X f (x)dµ(x). The set of k-times continuously differentiable functions over X is C k (X).

The indicator function of a set A ⊆ X is I A : X → {0, 1}, and has the values I A (x) = 0 for x ̸ ∈ A and I A (x) = 1 for x ∈ A. The measure of A with respect to µ ∈ M + (X) is defined as µ(A) = ⟨I A , µ⟩. The mass of a measure µ ∈ M + (X) is µ(X) = ⟨1, µ⟩, and µ is a probability measure if this mass is 1. The support of a measure µ is the set of all points x such that each open neighborhood N x of x obeys µ(N x ) > 0. The Dirac delta δ x ′ supported only at the point x ′ is a probability measure such that ⟨f, δ x ′ ⟩ = f (x ′ ) for all f ∈ C(X). Given two measures µ ∈ M + (X) and ν ∈ M + (Y ), the product measure µ ⊗ ν is the unique measure that satisfies

∀A ⊆ X, B ⊆ Y : (µ ⊗ ν)(A × B) = µ(A)ν(B). The pushforward of a function Q : X → Y along a measure µ(x) is Q # µ(y) and satisfies the relation ∀g ∈ C(Y ) : ⟨g(y), Q # µ(y)⟩ = ⟨g(Q(x)), µ(x)⟩.
The operator ∧ will be used to denote the minimum of two quantities (stopping times) as a ∧ b = min(a, b). The adjoint of a linear operator L :

X → Y is L † : Y * → X * .

Probability Tail Bounds and Value-at-Risk

Let ξ be a univariate probability measure ξ(ω) ∈ M + (R) for a coordinate ω ∈ R, with ⟨1, ξ⟩ = 1 and |⟨ω, ξ⟩|, ⟨ω 2 , ξ⟩ < ∞ (finite first and second moments). In this paper, we define the ϵ-VaR of ξ as follows:

V aR ϵ (ξ) = sup {λ ∈ R | ξ([λ, ∞)) ≥ ϵ} . (1)
Let σ 2 = ⟨ω 2 , ξ⟩ -⟨ω, ξ⟩ 2 be the variance of the probability distribution ξ.

The Cantelli bound for VaR is [START_REF] Cantelli | Sui confini della probabilita[END_REF] V aR ϵ (ξ

) ≤ σ 1/(ϵ) -1 + ⟨ω, ξ⟩ = V aR cant ϵ (ξ). ( 2a
)
The VP bound for the VaR is [START_REF] Vysochanskij | Justification of the 3σ rule for unimodal distributions[END_REF] V aR ϵ (ξ

) ≤ σ 4/(9ϵ) -1 + ⟨ω, ξ⟩ = V aR V P ϵ (ξ). ( 2b
)
The Cantelli bound is applicable for any probability distribution ξ(ω) and value ϵ ∈ [0, 1]. The VP bound is sharper than the Cantelli bound, but is only valid when ξ is unimodal and ϵ ≤ 1/6.

Stochastic Differential Equations

Let (Ω, F, P) be a probability space with time-indexed filtration F t , X ⊂ R n be a compact set, and w be n-dimensional Wiener process. An Itô SDE with a drift function f and diffusion function g is [START_REF] Gardiner | Stochastic methods[END_REF] dx = f (t, x)dt + g(t, x)dw.

In this paper, trajectories will start from an initial set X 0 ⊆ X and will remain within X in times t ∈ [0, T ] by virtue of stopping at the boundary ∂X. Define τ X as a stopping time (random variable) corresponding to the time at which the process (3) starting from X 0 first touches the boundary ∂X for the first time. A process of (3) starting from an initial condition x(0

) ∈ X 0 in times t ∈ [0, T ] is x(t) = x(0) + τ X ∧T t=0 f (t, x)dt + τ X ∧T t=0 g(t, x)dw. (4)
Solutions of ( 4) are unique if there exists finite constants C, D > 0 such that for all (t, x, x ′ ) ∈ [0, T ] × X 2 , the following Lipschitz and Growth conditions hold [START_REF]Stochastic Differential Equations: An Introduction with Applications[END_REF]:

D∥x -x ′ ∥ 2 ≥ ∥f (t, x) -f (t, x ′ )∥ 2 + ∥g(t, x) -g(t, x ′ )∥ 2 C(1 + ∥x∥ 2 ) ≥ ∥f (t, x)∥ 2 + ∥g(t, x)∥ 2 .
(

) 5
The Lipschitz and Growth conditions will hold if (f, g) are locally Lipschitz and the set X is compact. Distributions of the densities of (4) may be computed by solving a Fokker-Planck equation with absorbing boundary conditions on ∂X [START_REF] Fürth | Einige Untersuchungen über Brownsche Bewegung an einem Einzelteilchen[END_REF][START_REF] Hwang | The Fokker-Planck equation with absorbing boundary conditions[END_REF].

The generator L associated with the SDE is a linear operator that satisfies ∀v(t, x) ∈ C 2 ([0, T ]× X) [START_REF]Stochastic Differential Equations: An Introduction with Applications[END_REF]:

Lv(t, x) = ∂ t v + f (t, x) • ∇ x v + 1 2 g(t, x) T ∇ 2 xx v g(t, x). (6)
The ∇ 2 xx v term arises from the Itô Lemma. Let τ be a stopping time adapted to the filtration, defined by τ = τ X ∧ T . The occupation measure µ ∈ M + ([0, T ] × X) corresponding to the stopping time τ , initial distribution µ 0 ∈ M + (X 0 ), and dynamics (3

) is ∀A ⊆ [0, T ], B ⊆ X is µ(A × B) = X0 τ t=0 I A×B ((t, x(t | x 0 ))) dt dµ 0 (x 0 ). ( 7
)
The initial measure µ 0 ∈ M + (X 0 ), the occupation measure µ from [START_REF] Cantelli | Sui confini della probabilita[END_REF], and the terminal measure µ τ ∈ M + ([0, T ] × X) defined by following the SDE (3) from initial conditions x 0 ∼ µ 0 until the stopping time τ , are all related by Dynkin's formula [START_REF] Dynkin | Markov Processes[END_REF] ⟨v,

µ τ ⟩ = ⟨v(0, x), µ 0 (x)⟩ + ⟨Lv, µ⟩ ∀v ∈ C 2 . ( 8
)
Dynkin's formula is an SDE generalization of the Liouville equation for ODEs. Equation ( 8) may be equivalently written in weak form as

µ τ = δ 0 ⊗ µ 0 + L † µ. (9)
An expectation-maximizing optimal stopping problem for the SDE in (3) with a reward function of p(x) in the region [0, T ] × X, when starting at the initial condition x(0) ∼ µ 0 ∈ M + (X 0 ), is

P * = sup E µ0 [p(x(τ ))].
The work in [START_REF] Cho | Linear Programming Formulation for Optimal Stopping Problems[END_REF] presents an infinite-dimensional LP in measures to solve this stopping problem

p * = sup ⟨p, µ τ ⟩ (10a) µ τ = δ 0 ⊗ µ 0 + L † µ (10b) ⟨1, µ 0 ⟩ = 1 (10c) µ, µ τ ∈ M + ([0, T ] × X) (10d)
µ 0 ∈ M + (X 0 ). ( 10e
)
Any µ that is part of a feasible solution (µ, µ 0 , µ τ ) for (10b)-(10e) will be referred to as a relaxed occupation measure. Program (10) satisfies p * ≥ P * , and tightness (p * = P * ) is achieved under the assumptions of Lipschitz continuity and Growth (5), compactness of [0, T ] × X, and continuity of p(x).

Peak Value-at-Risk Estimation

This section will present the chance-peak problem statement, and will also derive the infinitedimensional SOCP to upper bound the chance-peak quantile statistic.

Problem Statement

Let ϵ ∈ [0, 1] be a value for the quantile statistic, X be a compact set, X 0 ⊆ X be a set of initial conditions, and (4) be the solution to an SDE evolving from x(0) ∈ X 0 that remains within X until it stops. For a given initial probability distribution µ 0 ∈ M + (X 0 ), and for all t ∈ [0, T ], let x(t) be the stochastic process of (4) at time t, and let µ t ∈ M + (X) be its probability distribution (with x(t) stopping at ∂X).

Assumptions

The following assumptions will be posed throughout this paper, A1 The set [0, T ] × X is compact and X 0 ⊆ X.

A2 The functions (f, g) satisfy (5).

A3 The state function p(x) is continuous on X.

A4 The initial measure µ 0 ∈ M + (X 0 ) is a given probability distribution (⟨1, µ 0 ⟩ = 1).

VaR Problem

Problem 3.1. The chance-peak problem to find the ϵ-VaR of p(x) is

P * = sup t * ∈[0,T ] V aR ϵ (p # µ t * ) (11a) dx = f (t, x)dt + g(t, x)dw (11b) from t = 0 until a stopping time of τ X ∧ t * (11c) x(0) ∼ µ 0 . (11d)
The pushforward p # µ t * from (11a) is the univariate probability distribution of p(x) at the state distribution x ∼ µ t * .

Tail-Bound Upper Bound

Let r be the constant factor multiplying σ in (2) such that

r cant = 1/(ϵ) -1 r V P = 4/(9ϵ) -1. ( 12
)
It is further assumed that the VP-bound will only be used if its conditions are satisfied (ϵ ≤ 1/6, unimodal). The distribution of p(x) with respect to the state distribution µ t * is univariate, for which the relation in ( 1) and the constants in (2) can be used to upper-bound on Problem 3.1. We will use the notation ⟨p 2 , µ t * ⟩ to refer to ⟨p(x) 2 , µ t * (x)⟩.

Problem 3.2. The tail-bound program that upper-bounds the chance-peak (11) with constant r is

P * r = sup t * ∈[0,T ] r ⟨p 2 , µ t * ⟩ -⟨p, µ t * ⟩ 2 + ⟨p, µ t * ⟩ (13a) dx = f (t, x)dt + g(t, x)dw (13b)
from t = 0 until a stopping time of τ X ∧ t * (13c)

x(0) ∼ µ 0 . (13d)

Nonlinear Measure Program

Problem 3.2 can be upper-bounded by an infinite-dimensional nonlinear program in a given initial probability distribution µ 0 , terminal measure µ τ , and relaxed occupation measure µ, using the generator L in (6) as

p * r = sup r ⟨p 2 , µ τ ⟩ -⟨p, µ τ ⟩ 2 + ⟨p, µ τ ⟩ (14a) µ τ = δ 0 ⊗ µ 0 + L † µ (14b) µ τ , µ ∈ M + ([0, T ] × X). (14c)
Theorem 3.3. Program 14 is an upper bound on (13) with p * r ≥ P * r under A1-A4. Proof. Let t * be a stopping time in [0, T ], and let x 0 ∈ X 0 be an initial condition. Measures (µ 0 , µ, µ τ ) that satisfy (14b) may be constructed from this (t * , x 0 ) by µ t * as the the state distribution of (4) at time t * given µ 0 , and µ as the occupation measure in [START_REF] Cantelli | Sui confini della probabilita[END_REF] associated to this SDE trajectory with distribution µ 0 . Because the feasible set to constraint (14b) contains measures induced by all possible provided SDE trajectories starting from µ 0 , it holds that p * r ≥ P * r .

Remark 1. The initial distribution µ 0 ∈ M + (X 0 ) may be optimized to find a supremal p * r over all probability distributions in X 0 by adding µ 0 as a variable and adding the constraint ⟨µ 0 , 1⟩ = 1 to (14).

Measure Second-Order Cone Program

The nonlinear measure program ( 14) may be recast as an infinite-dimensional convex SOCP.

Q 3 = {([s 1 , s 2 , s 3 ], κ) ∈ R 3 × R + | ∥s∥ 2 ≤ κ} is an SOC cone): sup (a,b)∈C a + r b -a 2 (15)
sup

(a,b)∈C, z∈R a + rz : ([1 -b, 2z, 2a], 1 + b) ∈ Q 3 . ( 16
)
Proof. The new variable z is introduced under the constraint √ b -a 2 ≥ z, implying that z 2 +a 2 ≤ b. The SOCP equivalence follows from the power-representation of √ b -a 2 from [START_REF] Alizadeh | Second-order cone programming[END_REF][START_REF] Lofberg | Working with square roots[END_REF], with the steps of

([1 -b, 2z, 2a], 1 + b) ∈ Q 3 (17a) (1 -b) 2 + 4(z 2 + a 2 ) ≤ (1 + b) 2 (17b) (1 + b 2 ) -2b + 4(z 2 + a 2 ) ≤ (1 + b 2 ) + 2b (17c) 4(z 2 + a 2 ) ≤ 4b. ( 17d
)
Theorem 3.5. An infinite-dimensional SOCP with the same optimal value and set of feasible solutions as [START_REF] Calafiore | Uncertain convex programs: randomized solutions and confidence levels[END_REF] given µ 0 is

p * r = sup rz + ⟨p, µ τ ⟩ (18a) µ τ = δ 0 ⊗ µ 0 + L † µ (18b) u = [1 -⟨p 2 , µ τ ⟩, 2z, 2⟨p, µ τ ⟩] (18c) (u, 1 + ⟨p 2 , µ τ ⟩) ∈ Q 3 (18d) µ, µ τ ∈ M + ([0, T ] × X), z ∈ R, u ∈ R 3 . (18e)
Proof. This results from an application of Lemma 3.4 to the objective term (14a). The optimization variables are now (µ τ , µ, z, u).

Corollary 1. Program (18) is convex.
Proof. The objective (18a) is affine in (z, µ τ ). Constraints (18b)-(18e) are convex (affine for (18b) and SOC for (18d)), ensuring convexity of (18).

Remark 2. Problem (18) has an infinite-dimensional affine constraint in (18b) and a finitedimensional SOC constraint in (18d).

This section will upper-bound (18) utilizing a converging hierarchy of SDPs of increasing size.

Review of Moment-SOS Hierarchy

All content from this subsection is referenced from [START_REF] Lasserre | Moments, Positive Polynomials And Their Applications[END_REF]. For any multi-indexed sequence m = {m α } α∈N n ∈ R N n , we define the Riesz functional L m : R[x] → R as follows:

p(x) = α∈N n p α x α -→ L m p = α∈N n p α m α . ( 19a
)
Let µ ∈ M + (X) be a measure. The α-moment of µ for α ∈ N n is m α = ⟨x α , µ⟩. The collection of moments m = {m α } α∈N n is a moment sequence and has the following property:

∀p ∈ R[x], L m p = ⟨p, µ⟩. (19b)
A key result to build the moment-SOS hierarchy is the characterization of sequences m that correspond to moment sequences, i.e. such that (19b) holds for some X ⊂ R n and some µ ∈ M + (X):

For m ∈ R N n and h ∈ R[x], we define the localizing bilinear functional

L hm : R[x] × R[x] -→ R by L hm = (p, q) -→ L m (hpq). ( 20a
)
Equipping R[x] with a linear basis (e i ) i∈N (e.g. e i (x) = x αi with {α i } i∈N = N n , an ordering of monomials such that |α i | < |α j | ⇒ i < j) yields an infinite size matrix representation of L hm , which we call the localizing matrix M[hm] = (L m (h e i e j )) i,j∈N . For instance, if h(x) = β∈N n h β x β , using a basis of monomials with nondecreasing degrees yields, for all i, j ∈ N

M[hm] i,j = L m (hx αi x αj ) = β∈N n h β m αi+αj +β . (20b)
A Basic Semialgebraic (BSA) set is a set defined by a finite number of bounded-degree inequality constraints such as

K = {x | h k (x) ≥ 0 : k = 1..N c }. Assuming "ball constraints" h 1 (x) = 1 and h Nc (x) = R -∥x∥ 2 2 (this can always be enforced if K is compact in R n , up to adding redundant constraints), m ∈ R N n has a representing measure µ ∈ M + (K) such that (19b) holds if, for all k = 1..N c , the bilinear functional L h k m is positive semidefinite, i.e. ∀p ∈ R[x], k = 1..N c , L m (h k p 2 ) ≥ 0. (21a) or, equivalently, ∀d ∈ N, k = 1..N c , M d [h k m] ⪰ 0. ( 21b
)
where

M d [h k m] is the top left block of size n+d d of M[h k m],
which corresponds to the matrix representation of L hm in the finite dimensional space R[x] ≤d .

For notational convenience, we define the block diagonal synthetic matrix

M d [Km] = diag(M d-⌈d k /2⌉ [h k m]) k=1..Nc (22)
where d k = deg(h k ). This synthetic matrix has two important properties, deduced from (20b) and (21b):

• it exactly involves all the terms m α for |α| ≤ 2d

• (21) holds if and only if M d [Km] ⪰ 0 for all d ∈ N.
The process of increasing the degree d → ∞ when posing Positive Semidefinite (PSD) constraints on M d [Km] is called the moment-SOS hierarchy.

Moment Program

The following assumptions are required to utilize the moment-SOS hierarchy in approximating (18): A5 The sets X 0 and X are BSA sets with ball constraints. A6 The functions f (t, x), g(t, x) are polynomial vector fields and p(x) is a polynomial.

Given an initial measure µ 0 ∈ M + (X 0 ), let (m, m τ ) be moment sequences corresponding to the measures (µ, µ τ ) respectively. For each monomial x α t β with α ∈ N n , β ∈ N, define the operator D αβ (m, m τ ) as the moment counterpart of the operator involved in Dynkin's formula ( 9)

D αβ (m, m τ ) = m τ αβ -L m (L(x α t β )). (23)
Define the dynamics degree D as

D = d + ⌈max(deg f -1, 2 deg g -2)/2⌉. (24)
so that for (α, β)

∈ N n+1 , |α| + β ≤ 2d ⇒ deg(L(x α t β )) ≤ 2D.
p * r,d = max rz + L m τ p (25a) z ∈ R, m ∈ R ( 2D+n+1 n+1 ) , m τ ∈ R ( 2d+n+1 n+1 ) (25b) D αβ (m, m τ ) = δ β0 ⟨x α , µ 0 ⟩ ∀(α, β) ∈ N n+1 s.t. |α| + β ≤ 2d (25c) s = [1 -L m τ (p 2 ), 2z, 2L m τ p] (25d)
(s, 1 + L m τ (p 2 )) ∈ Q 3 (25e) M d [([0, T ] × X)m τ ] ⪰ 0 (25f) M D [([0, T ] × X)m] ⪰ 0, (25g)
where δ β0 denotes the Kronecker symbol that is 1 if β = 0 and 0 otherwise. Note that constraint (25c) is a finite-dimensional truncation of the infinite-dimensional (18b).

The following boundedness result is required to ensure convergence: Lemma 4.2. All of (µ, µ τ , z) are bounded in (18) under A1-A3.

Proof. A sufficient condition for a measure to be bounded (in the sense that all moments are bounded) is that it has finite mass and is supported on a compact set. Compactness of [0, T ] × X holds by A1. Assumption A4 imposes that ⟨1, µ 0 ⟩ = 1. By substituting v(t, x) = 1 (18b), it holds that ⟨1, µ τ ⟩ = ⟨1, µ 0 ⟩ = 1. Performing the same step with v(t, x) = t yields T ≥ ⟨t, µ τ ⟩ = ⟨1, µ⟩. It therefore holds that ⟨p, µ τ ⟩ and ⟨p 2 , µ τ ⟩ are bounded. The SOC constraint (18d) ensures that z is finite, demonstrating that all variables are bounded. Proof. This convergence will occur by Corollary 8 of [START_REF] Tacchi | Convergence of Lasserre's hierarchy: the general case[END_REF] (when extending to the finite-dimensional SOC case) through Lemma 4.2.

Remark 3. The relation p * d ≥ p * r ≥ P * r will still hold when [0, T ] × X is noncompact (violating A1 and A5), but it may no longer occur that lim d→∞ p * r,d = p * r (the conditions Lemma 4.3 will no longer apply).

Computational Complexity

In Problem [START_REF] Fantuzzi | Bounds for Deterministic and Stochastic Dynamical Systems using Sum-of-Squares Optimization[END_REF], the computational complexity mostly depends on the number and size of the matrix blocks involved in LMI constraints (25f,25g), which in turn depend on the number and degrees of polynomial inequalities describing X (the higher d k = deg(h k ), the smaller M d-⌈d k /2⌉ [hm]). At order-d, the maximum size of localizing matrices is n+1+D D . Problem [START_REF] Fantuzzi | Bounds for Deterministic and Stochastic Dynamical Systems using Sum-of-Squares Optimization[END_REF] must be converted to SDP-standard form by introducing equality constraints between the entries of the moment matrices in order to utilize symmetric-cone Interior Point Methods (e.g., Mosek [START_REF] Aps | The MOSEK optimization toolbox for MATLAB manual[END_REF]). The per-iteration complexity of an SDP involving a single moment matrix of size n+d d scales as n 6d [START_REF] Lasserre | Pricing a class of exotic options via moments and SDP relaxations[END_REF]. The scaling of an SDP with multiple moment and localizing matrices generally depends on the maximal size of any PSD matrix. In our case, this size is at most n+1+d d with a scaling impact of (n + 1) 6d . The complexity of using this chance-peak routine increases in a jointly polynomial manner with d and n.

Numerical Examples

MATLAB (2022a) code to replicate experiments is available at https://github.com/Jarmill/ chance_peak. Dependencies include Mosek [START_REF] Aps | The MOSEK optimization toolbox for MATLAB manual[END_REF] and YALMIP [START_REF] Lofberg | YALMIP : a toolbox for modeling and optimization in MATLAB[END_REF]. Monte Carlo (MC) sampling to empirically find VaR estimates is conducted over 50,000 SDE paths under antithetic sampling with a time spacing of ∆t = 10 -3 . All experiments contain a table of chance-peak bounds as well as solver-times to compute these bounds.

Two States

Example 1 of [START_REF] Prajna | Stochastic safety verification using barrier certificates[END_REF] is the following two-dimensional cubic polynomial SDE

dx = x 2 -x 1 -x 2 -1 2 x 3 1 dt + 0 0.1 dw. (26)
This example performs chance-peak maximization of p(x) = -x 2 starting at the point (Diracdelta initial measure µ 0 ) X 0 = [1, 1] with X = [-1, 2] × [-1, 1.5] and T = 5. Trajectories of (26) are Figure 1: Trajectories of ( 26) with ϵ = 0.5 (dashed red) and ϵ = 0.15 (solid red) bounds displayed in cyan in Figure 1 starting from the black-circle X 0 , and four of these trajectories are marked in non-cyan colors. The ϵ = 0.5 row of Table 1 displays the bounds on the mean distribution as solved through finite-degree SDP truncations of [START_REF] Lagoa | Probabilistically Constrained Linear Programs and Risk-Adjusted Controller Design[END_REF]. The bounds at ϵ = {0.15, 0.1, 0.05} are obtained through the VP expression in (2b) and solving the SDPs obtained from [START_REF] Fantuzzi | Bounds for Deterministic and Stochastic Dynamical Systems using Sum-of-Squares Optimization[END_REF]. The dotted and solid red lines in Figure 1 are the ϵ = 0.5 and ϵ = 0.15 bounds respectively at order 5.

Three States

An SDE modification of the Twist system from [START_REF] Miller | Bounding the Distance of Closest Approach to Unsafe Sets with Occupation Measures[END_REF] This second example performs chance-peak maximization of p(x) = x 3 starting at the point X 0 = [0.5, 0, 0] with X = [-0.6, 0.6] × [-1, 1] × [-1, 1.5] and T = 5. VP bounds from solving the SDEs from ( 10) and ( 25) are recorded in Table 3 in the same manner as in Table 1. Figure 2 plots trajectories and bounds of ( 27) starting from the black-circle X 0 point, with four of these trajectories visibly distinguished. The solid red plane in Figure 2 is the ϵ = 0.15 bound on x 3 at order 6, and the transluscent red plane is the ϵ = 0.5 bound on x 3 (also at order 6).

is, dx =   -2.5x 1 + x 2 -0.5x 3 + 2x 3 1 + 2x 3 3 -x 1 + 1.5x 2 + 0.5x 3 -2x 3 2 -2x 3 3 1.5x 1 + 2.5x 2 -2x 3 -2x 3 1 -2x 3 2   dt +   0 0 0.1   dw. (27)

Conclusion

This paper considered the chance-peak problem, which involved finding upper bounds on the quantiles of state functions p(x) achieved by SDE systems. The true (1 -ϵ)-quantile statistic P * (11) is upper-bounded by the Cantelli/VP approximation P * r [START_REF] Han | Non-Gaussian Risk Bounded Trajectory Optimization for Stochastic Nonlinear Systems in Uncertain Environments[END_REF], which in turn is upper bounded by an infinite-dimensional SOCP p * r (18) and its moment-SOS finite-dimensional SDPs yielding p * r,d with lim d→∞ p * r,d = p * d . Each of these upper-bounds contribute valuable information towards the analysis of SDEs.

Future work includes finding conditions under which the measure-based upper-bounding does not add conservatism (e.g. cases where p * r = P * r ), and utilizing higher-moment tail-probability inequalities to obtain closer estimates to the VaR [START_REF] He | Bounding Probability of Small Deviation: A Fourth Moment Approach[END_REF]. Other work involves studying the duality structure of [START_REF] Lasserre | Moments, Positive Polynomials And Their Applications[END_REF] with respect to the achievement of strong duality. The techniques introduced in this paper can also be extended to application domains such as distance estimation [START_REF] Miller | Bounding the Distance of Closest Approach to Unsafe Sets with Occupation Measures[END_REF], Lévy processes [START_REF] Kashima | An optimization approach to weak approximation of Lévy-driven stochastic differential equations[END_REF], and exit-time statistics [START_REF] Henrion | Moment-SOS hierarchy and exit time of stochastic processes[END_REF]. Another avenue involves developing stochastic optimal control strategies to minimize quantile statistics.

Lemma 3 . 4 .

 34 Let J r (a, b) = r √ b -a 2 + a be the objective (14a) with a = ⟨p(x), µ τ ⟩ and b = ⟨p(x) 2 , µ τ ⟩. For any convex set C ∈ R × R + with (a, b) ∈ C, the following pair of programs have the same optimal value (in which

Problem 4 . 1 .

 41 For d ≥ deg(p), the order-d moment problem that upper-bounds problem[START_REF] Lasserre | Moments, Positive Polynomials And Their Applications[END_REF], given µ 0 is

Theorem 4 . 3 .

 43 Under A1-A6, the optima in (25) will converge to (18) as lim d→∞ p * r,d = p * r .

Figure 2 :

 2 Figure 2: Trajectories of (27) with ϵ = 0.5 (dashed red) and ϵ = 0.15 (solid red) bounds

Table 1 :

 1 Chance-Peak estimation of the Stochastic Flow System (26) to maximize p(x) = -x 2

	order	2	3	4	5	6	MC
	ϵ = 0.5	0.8817 0.8773 0.8747 0.8745 0.8744 0.8559
	ϵ = 0.15 1.2210 1.1817 1.1689 1.1657 1.1642 0.9142
	ϵ = 0.1	1.5009 1.4520 1.4361 1.4323 1.4303 0.9279
	ϵ = 0.05 2.1306 2.0613 2.0387 2.0332 2.0305 0.9484

Table 2 :

 2 Solver time (seconds) to compute Table1

	order	2	3	4	5	6
	ϵ = 0.5	0.417 0.431 1.963 1.659 3.641
	ϵ = 0.15 0.216 0.325 1.592 4.178 7.464
	ϵ = 0.1	0.213 0.316 1.651 1.339 4.225
	ϵ = 0.05 0.222 0.366 0.936 2.446 5.298

Table 3 :

 3 Chance-Peak estimation of the Stochastic Twist System (27) to maximize p(x) = x 3

	order	2	3	4	5	6	MC
	ϵ = 0.5	0.9100 0.8312 0.8231 0.8211 0.8203 0.7206
	ϵ = 0.15 1.1729 1.0813 1.0171 0.9891 0.9755 0.7685
	ϵ = 0.1	1.4358 1.3180 1.2217 1.1859 1.1694 0.7801
	ϵ = 0.05 2.0288 1.8497 1.6866 1.6299 1.5894 0.7970

Table 4 :

 4 Solver time (seconds) to compute Table3

	order	2	3	4	5	6
	ϵ = 0.5	0.526 1.880 5.213 28.269 81.736
	ϵ = 0.15 0.397 1.117 3.876 23.302 132.539
	ϵ = 0.1	0.210 1.080 6.313 27.829 116.733
	ϵ = 0.05 0.185 1.352 4.999 20.863 137.833

Acknowledgements J. Miller and M. Sznaier were partially supported by NSF grants CNS-1646121, ECCS-1808381 and CNS-2038493, AFOSR grant FA9550-19-1-0005, and ONR grant N00014-21-1-2431. J. Miller was in part supported by the Chateaubriand Fellowship of the Office for Science & Technology of the Embassy of France in the United States and by the International Student Exchange Program from AFOSR. M. Tacchi was supported by the French company Réseau de Transport d' Électricité, as well as the Swiss National Science Foundation under the "NCCR Automation" grant n • 51NF40 180545.


Aleksey Buzmakov
email: avbuzmakov@hse.ru

Tatiana Makhalova
email: t.makhalova@gmail.com

Sergei O Kuznetsov
email: skuznetsov@hse.ru

Amedeo Napoli
email: amedeo.napoli@loria.fr

Δ-Closure Structure for Studying Data Distribution

Keywords: pattern mining, closed sets, equivalence class, generators, data distribution

In this paper, we revisit pattern mining and study the distribution underlying a binary dataset thanks to the closure structure which is based on passkeys, i.e., minimum generators in equivalence classes robust to noise. We introduce ∆-closedness, a generalization of the closure operator, where ∆ measures how a closed set differs from its upper neighbors in the partial order induced by closure. A ∆-class of equivalence includes minimum and maximum elements and allows us to characterize the distribution underlying the data. Moreover, the set of ∆classes of equivalence can be partitioned into the so-called ∆closure structure. In particular, a ∆-class of equivalence with a high ∆ is supported by more observations and thus is more stable. In the experiments, we study the ∆-closure structure of several real-world datasets and show that this structure is very stable for large ∆ and does not substantially depend on the data sampling used for the analysis.

I. INTRODUCTION

In this paper, we are interested in pattern or itemset mining in tabular data. There is a considerable work on this subject, especially regarding algorithms and search for interesting patterns [START_REF]Frequent Pattern Mining[END_REF]. Here we rather focus on the distribution underlying the dataset under study thanks to closed itemsets, their equivalence classes, and the associated generators.

Many pattern mining approaches produce a particular set of patterns, which provides a certain "view" of the intrinsic structure underlying the data. However, this view, usually, is not directly related to the distribution underlying the data. In this paper we rely on the closure structure which was recently introduced [START_REF] Makhalova | Introducing the closure structure and the GDPM algorithm for mining and understanding a tabular dataset[END_REF] and which characterizes the intrinsic structure of a dataset. The closure structure reveals the distribution of the itemsets in the data in terms of frequency and also of stability, i.e., how a closed set depends on its content, and in addition supports an interpretation of the content of a dataset. Here we propose a methodology for computing and understanding the intrinsic structure of a dataset in two main ways: (i) we revisit the closure structure, which reveals the distribution of the itemsets in the data in terms of frequency, (ii) we generalize Sergei O. Kuznetsov would like to thank Basic Research Program of the HSE University for its support.

the closure structure to Δ-closure structure, which is a more robust closure structure related to the distribution underlying the data. The closure structure and the Δ-closure stracture support an interpretation of the content of a dataset. The closure structure is based on "closed itemsets" and minimum elements of their "equivalence classes" -related to levels-, which are computed independently of any interestingness measure or set of constraints. The closure structure and its levels provide a representation of the complexity of the content of a dataset. We propose a formalization of the closure structure in terms of Formal Concept Analysis [START_REF] Ganter | Formal Concept Analysis -Mathematical Foundations[END_REF], which is well adapted to the study of closed itemsets, equivalence classes, and data topology in terms of closed sets.

Then we generalize the notion of closure to ∆-closedness which allows us to work with stability along with frequency. Actually, ∆ measures how much a closed set differs from its upper neighbors in the partial order of closed sets. A ∆class of equivalence allows us to characterize the distribution underlying the data: (i) when ∆ is large, there are only a few ∆-classes of equivalence whose elements are very stable, (ii) when ∆ is small, the number of ∆-classes increases while the related information is less stable and depends on a smaller number of elements. This allows us to study stable patterns which are robust against noise and also to dynamic changes in the data.

In the experiments, we compute the closure structure of a number of public datasets and show how ∆-closedness is taken into account. The closure structure may roughly be of three types, where the first levels of the structure are the most interesting, stable, and interpretable in standard ("typical") and almost standard ("intermediate") datasets, while in nonstandard datasets ("deterministic"), levels are more separated and less easily interpretable. For example, the closure structure allows us to determine the levels where the itemsets are the most diverse. It shows how the frequency and the stability are distributed among the levels, and when a search for interesting itemsets can be stopped without major loss. Moreover, the ∆closure structure is very stable for large ∆ and does not substantially depend on the data sampling used for the analysis. To the best of our knowledge, such kind of study and visualisation of data is rather unique. The proposed methodology and the related tools are useful for simple navigation in the pattern search space.

The paper has the following structure. Section II presents related work and motivation. In Section III we firstly introduce ∆-closure with the related equivalence classes and passkeys. We also provide a concrete example and we make precise an algorithm for discovering the ∆-closure and the related elements. Then in Section IV we present and discuss experiments about the computing of the ∆-closure structure and the stability of passkeys. We also show how this structure can be visualized and interpreted. Then we conclude the paper and give directions for future work.

II. RELATED WORK AND MOTIVATION

Following the research directions discussed in [START_REF] Boley | Efficient Discovery of Interesting Patterns Based on Strong Closedness[END_REF] and [START_REF] Makhalova | Introducing the closure structure and the GDPM algorithm for mining and understanding a tabular dataset[END_REF], we are interested in discovering closed and stable patterns together with their related equivalence classes and minimum generators (passkeys). As in [START_REF] Boley | Efficient Discovery of Interesting Patterns Based on Strong Closedness[END_REF], we introduce a ∆ parameter for studying ∆-closure and the robustness of the discovered patterns against noise, and as well the stability of the associated equivalence classes and their characteristic elements such as passkeys. Moreover, as in [START_REF] Makhalova | Introducing the closure structure and the GDPM algorithm for mining and understanding a tabular dataset[END_REF] we are also interested in investigating the data distribution through the closure structure. The objective is to build a kind of "internal picture" of the datasets under study and to check at which level -related to closure and ∆-closureone can discover the most interesting patterns, and as well interesting implications or association rules. In addition, in the present paper, we introduce Δclosure structure based on Δ-closure operator which subsumes ordinary closure.

The interest in data distribution is frequently appearing in the literature. Indeed, even the standard frequent itemset mining indirectly focuses on properties of the distribution. Such interest in the distribution can be stated more formally if the itemset frequency is verified as a statistical hypothesis [START_REF] Riondato | Finding the True Frequent Itemsets[END_REF]. Some papers study how to take noise into account in pattern mining, which can also be considered as indirect characterization of the distribution underlying the data. In particular, δ-free sets [START_REF] Boulicaut | Free-sets: a condensed representation of boolean data for the approximation of frequency queries[END_REF] are special classes of patterns that encode association rules with a high confidence -only few counterexamples are accepted-for allowing identification of noiseresistant patterns [START_REF] Pensa | Towards Fault-Tolerant Formal Concept Analysis[END_REF]. Other studies were also performed to find out classes of patterns resistant to noise [START_REF] Klimushkin | Approaches to the Selection of Relevant Concepts in the Case of Noisy Data[END_REF], [START_REF] Yang | Efficient discovery of error-tolerant frequent itemsets in high dimensions[END_REF].

Still about data distribution, in [START_REF] Ramesh | Feasible itemset distributions in data mining: theory and application[END_REF] the authors study distribution of sizes of frequent and maximal frequent itemsets in a database. In [START_REF] Flouvat | A new classification of datasets for frequent itemsets[END_REF] the authors are mostly interested in the distributions related to the border between frequent and infrequent itemsets together with the distribution of three concise representations: frequent closed, frequent free, and frequent essential itemsets. In [START_REF] Van | Fast Estimation of the Pattern Frequency Spectrum[END_REF] the authors propose an efficient approach to estimation of the number of frequent patterns for arbitrary minsup thresholds. In our study we assume that "important" patterns enabling better interpretation are not necessarily frequent, but are certainly stable to noise in data. So, here we study empirical distribution of noise tolerance in terms of ∆-stability introduced below.

Computational complexity remains an important challenge in itemset mining. The number of itemsets can be exponential w.r.t. the dataset size. Focusing on closed itemsets allows for a substantial reduction of this number by replacing an equivalence class, i.e., a whole class of itemsets having the same support, by the largest one which is the closed itemset [START_REF] Pasquier | Efficient mining of association rules using closed itemset lattices[END_REF]. For dealing with large amounts of closed itemsets, it is possible to evaluate itemsets thanks to quality metrics [START_REF] Kuznetsov | On interestingness measures of formal concepts[END_REF], [START_REF] Vreeken | Interesting Patterns[END_REF]. For example ∆-closedness evaluates the robustness of patterns and the corresponding closure operation [START_REF] Boley | Efficient Discovery of Interesting Patterns Based on Strong Closedness[END_REF]. Roughly speaking, a ∆-closed set cannot be augmented by any further item without decreasing its support by at least ∆. Actually, ∆closed sets allow to capture interesting itemsets that are stable, i.e., robust to noise and changes in the data. A very similar notion is δ-tolerance closed itemsets which are defined w.r.t. itemset support rather than an absolute change in support for Δ-closure [START_REF] Cheng | δ-Tolerance Closed Frequent Itemsets[END_REF]. An alternative to exhaustive enumeration of itemsets is based on "sampling" [START_REF] Dzyuba | Flexible constrained sampling with guarantees for pattern mining[END_REF] and on a gradual search for itemsets according to an interestingness measure or a set of constraints [START_REF] Smets | Slim: Directly Mining Descriptive Patterns[END_REF]. Such algorithms usually output a rather small set of itemsets while they may provide only an approximate solution.

The papers discussed above do not address the problem of discovering an intrinsic structure underlying the set of closed itemsets, and -more generally-the dataset itself. By contrast, in this paper, we define such a level-wise structure for representing this intrinsic structure of the datasets. Moreover, we build this structure by means of ∆-closedness operator and we study the stability of the related ∆-closure levels w.r.t. the size of ∆. The Δ-closure reveals the internal structure of a dataset and as well relates this structure to the distribution underlying the dataset.

III. FORMALISM

A dataset is modeled as a formal context [START_REF] Ganter | Formal Concept Analysis -Mathematical Foundations[END_REF], which is a triple K = (G, M, I), where G is a set of objects, M is a set of attributes, and I ⊆ G × M is an incidence relation such that (g, m) ∈ I if object g has attribute m. Two derivation operators, both denoted by (•) , are defined for A ⊆ G and B ⊆ M as follows:

A = {m ∈ M | ∀g ∈ A : gIm} (1) B = {g ∈ G | ∀m ∈ B : gIm} (2)
Intuitively, A is the set of attributes common to objects in A, while B is the set of objects which have all attributes in B. Sets A ⊆ G, B ⊆ M , such that A = A and B = B , are closed sets. Moreover, (•) is a closure operator equivalent to instance counting closure in data mining as introduced in [START_REF] Pasquier | Efficient mining of association rules using closed itemset lattices[END_REF]. For A ⊆ G, B ⊆ M , a pair (A, B) such that A = B and B = A, is called a formal concept, then A and B are closed sets and called extent and intent, respectively.

m 1 m 2 m 3 m 4 m 5 m 6 m 7 m 8 m 9 a b c d e f g h i g 1 × × × × × g 2 × × × × × g 3 × × × × g 4 × × × × g 5 × × × × g 6 × × × × g 7 × × × × × g 8 × × × × × g 9 × g 10 ×
Fig. 1: A toy binary dataset where 'x' in cell (i, j) indicates that object g i is related to attribute m j .

A. Δ-classes of equivalence

Definition III.1. An itemset B is called Δ-closed if for any attribute m ∈ M : |B | -|(B ∪ {m}) | ≥ ∆ ≥ 1. (3)
If ∆ = 1, then the itemset is just closed w.r.t. object counting [START_REF] Boley | Efficient Discovery of Interesting Patterns Based on Strong Closedness[END_REF], [START_REF] Pasquier | Efficient mining of association rules using closed itemset lattices[END_REF]. In [START_REF] Boley | Efficient Discovery of Interesting Patterns Based on Strong Closedness[END_REF] it was also shown that Δ-closedness is associated with a closure operator that we call Δ-closure and that works as follows. Definition III.2. Given an itemset X, the equivalence class Equiv ∆ (X) of X is the set of all itemsets whose ∆-closure is equal to the ∆-closure of X, i.e.,

Equiv ∆ (X) = {Y ⊆ M | φ ∆ (Y ) = φ ∆ (X)}. (4)
Some elements of a class of equivalence can be highlighted since they have special properties [START_REF] Makhalova | Introducing the closure structure and the GDPM algorithm for mining and understanding a tabular dataset[END_REF]. A concrete example about Δ-closure, Δ-keys, and Δpasskeys is given in Section III-B. It can be checked that Δclasses of equivalence with higher ∆ make coarser partitions of the set of concept intents (i.e., ∆-closed itemsets for ∆ = 1)) than Δ-classes of equivalence with smaller Δ, i.e., (∀∆ ≥ 2)(∀X ⊆ M )(∃X 1 , . . . , ∃X k )Equiv ∆ (X) = Equiv ∆-1 (X i ). Moreover, we can associate every Δ-class of equivalence with the size of its passkeys. This size corresponds to a level in the class of equivalence. Then for any Δ the complexity of the dataset can be captured by means of the distribution of the classes of equivalence within the levels, called here the level structure. Now, given an itemset X, let us define the following measures:

Definition III.3. A Δ-key X ∈ Equiv ∆ (B) is any minimal (w.r.t. subset relation) itemset in Equiv ∆ (B).
∆(X) = max{0 ≤ ∆ ≤ |G| | X is Δ-closed} (5) ∆ key (X) = max{0 ≤ ∆ ≤ |G| | X is a Δ-key} (6) ∆ passkey (X) = max{0 ≤ ∆ ≤ |G| | X is a Δ-passkey} (7)
All these measures capture the maximal Δ for which the itemset X preserves a certain property, i.e., being closed, being a key, or being a passkey. We will call these values ∆-values of closed itemsets, keys, and passkeys, respectively. Since any passkey is a key, ∆ passkey (X) ≤ ∆ key (X). It is clear from the definitions that for all 1 ≤ ∆ ≤ ∆ key (X), X is a Δ-key, and for all 1 ≤ ∆ ≤ ∆ passkey (X), X is a Δ-passkey.

Let us now relate these measures to the variations existing in a dataset. More precisely, is it possible to discover the itemsets that are likely to be preserved when a dataset is changed? Proposition III.1. Let ∆ passkey (X) = δ. Then at least δ objects should be removed from the dataset K to obtain a subdataset K s ⊆ K such that X is not a passkey in K s .

Proposition III.2. Let ∆ key (X) = δ. Then at least δ objects should be removed from the dataset K to obtain a subdataset K s ⊆ K such that X is not a key in K s .

Proposition III.3. Let ∆(X) = δ. Then at least δ objects should be removed from the dataset K to obtain a subdataset

K s ⊆ K such that X is not closed in K s .
The proofs of these propositions can be found in a technical report at http://arxiv.org/abs/2210.06926. These propositions show that, the higher is Δ, the deeper should the dataset be modified for removing a certain property of the dataset. Thus, if the structure of the dataset is captured with itemsets of high Δ, then the discovered itemsets are stable w.r.t. object removal. This fact is related to the distribution underlying the dataset the objects are taken from.

As a final remark, let us notice that an itemset X is a Δ-key

only if ∀Y ⊂ X (|Y | -|X | ≥ ∆)
. This is only a necessary condition. For example, in the following context with 3 objects G = {g 1 , g 2 , g 3 }, {g 1 } = ∅, {g 2 } = {m 1 }, and {g 3 } = {m 1 , m 2 }, the aforementioned condition is satisfied for {m 2 } and ∆ = 2. However, {m 2 } is not a Δ-key, but a δ-free set [START_REF] Boulicaut | Free-sets: a condensed representation of boolean data for the approximation of frequency queries[END_REF]. This shows that not every δ-free set is a Δ-key.

B. Example

Figure 1 displays a toy formal context. Every object except objects 1, 2, 3, and 4, has its own attribute set and thus the size of their passkeys is 1. The closed itemsets (intents) along with their supports (extents) are shown in Fig. 2. For brevity sake, we write 123 instead of {1, 2, 3} and so on. The concepts with extents 123456 and 12478 have passkeys equal to c and e, respectively. The passkey for concept with extent 12345678 is ab and the passkey for the top concept is ∅. The concept with extent 12 has a passkey larger than 2, namely cde. In Fig. 2, the Δ of the intent or Δ-measure is indicated in the square brackets for every concept when it is different from 1. The concept intents are the maximal elements in the Δclasses of equivalence for Δ not larger than their Δ-measure. Now consider the Δ-classes of equivalence for ∆ = 2. There are only 4 concepts with a Δ-measure not smaller than 2, thus, there are only 4 classes of equivalence for ∆ = 2 plus the "technical" class of equivalence at the bottom. These classes of equivalence with ∆ = 2 are formed by joining smaller classes of equivalence, i.e., the classes of equivalence with ∆ = 1 which are related to single concepts. Then in Fig. 2, the classes of equivalence with ∆ = 2 are shown in green areas, while the simple concepts, i.e., closed sets, correspond to classes of equivalence with ∆ = 1.

Now, if we consider Δ-classes of equivalence for ∆ = 3, there are just two closures abc and M . These equivalence classes are surrounded by the dotted line in the diagram. The Δ-class containing the closed set abc contains the concepts with the intents abc, ab, a, b, ∅. The other concepts belong to the Δ-class of equivalence whose maximal element is M . If we further increase Δ, the whole lattice will collapse into a single class of equivalence.

C. Computational considerations

The algorithm for computing Δ-equivalence classes is sketched in Algorithm 1 while details are provided in Algorithm 2. Firstly, we compute the concepts and the associated concept levels using GDPM algorithm [START_REF] Makhalova | Introducing the closure structure and the GDPM algorithm for mining and understanding a tabular dataset[END_REF], which iterates over keys searching for passkeys. Then, for every concept (A, B) we compute the value of ∆(B) as given by equation [START_REF] Dzyuba | Flexible constrained sampling with guarantees for pattern mining[END_REF]. Finally, for every concept (A, B) we compute the value of ∆ passkey (B) as given by equation ( 7). This procedure is made precise in Algorithm 2.

At the beginning the value of ∆ passkey (B) is initialized to 1 for every concept (A, B). Then, since any (∆ + 1)-class of Algorithm 1: A general algorithm for finding Δ-values for all concept intents and for their passkeys.

Data: K is the dataset Result: A set of concepts associated with their levels and the values of Δ for the concept intent and for the concept passkey

1 {c i ⇒ level i } = GDPM(K); 2 {c i ⇒ ∆ cls (c i ), Ref (c i ) } = ComputeDeltas({c i }); 3 {c i ⇒ ∆ passkey (c i )} = ComputePKDeltas({c i ⇒ level i , ∆ cls (c i ), Ref (c i ) });
equivalence is a union of Δ-classes of equivalence, we close every Δ-closed element with (∆ + 1)-closure (lines 10-19).

Finally, if two Δ-classes of equivalence have the same (∆+1)closure, we select one passkey among the possible (∆ + 1)passkeys (lines 20-29). The computational time complexity of Algorithm 2 is O(|G| • log(|G|) • |{c i }|), since we need to iterate over all possible thresholds of Δ and for any given threshold we need to find the Δ-closure of every concept, which can be found in at most log(|G|) time, since on every new iteration of loop 12-17 the distance to the closure from the concept is either doubled or the closure is found.

IV. EXPERIMENTS

The experiments are carried out using a machine with Intel Core i5 CPU, 16GB of RAM and Nvidea GeForce RTX 3080 video card operated under Ubuntu 20.04 operating system. Table I shows datasets with total computation time larger than 10 seconds. It also shows the number of objects, attributes, and closed itemsets. The computation time is divided into the computation time of the level structure with GDPM and the time for computing Δ-classes of equivalence. We can observe that the most important factor in the computation time is the number of closed itemsets in the dataset, the second factor being the size of the dataset.

A. Stability of passkeys

The goal of the experiments is to verify that Δ-closure structures are becoming more robust when Δ is increasing. In particular, Δ-closure structures when Δ is greater than 1 are more robust than closure structures or 1-closure structures.

For any Δ-equivalence class we distinguish its minimum elements called the Δ-passkeys and the maximum one called the Δ-closed itemset. They are both characterized by Δ, where the higher value for Δ, the more stable are passkeys/closed itemsets w.r.t. noise. A passkey can be considered as a most concise "definition" of a Δ-closed itemset.

In this section we study how stable Δ-equivalence classes are distributed in the dataset. To analyze the datasets we use the following observations:

(i) the equivalence classes with small passkeys are easier to mine by the algorithm that computes the closure Algorithm 2: An algorithm for computing Δ of passkeys for every concept. Data: {c i } is a set of concepts with associated levels, Δ-value for their intent and the most closest child Result: The value of Δ for passkeys of every concept 1 def ComputePKDeltas({c i }): are more stable w.r.t. noise. We analyzed 25 datasets from LUCS-KDD repository (their parameters are given in the supplementary material) paying attention to the properties mentioned above. Among all the datasets under analysis we discovered three main types of behaviors -typical, deterministic, and two kinds of intermediates-which are displayed in the columns of Fig. 3. the attribute frequency distribution and does not contain any information about interaction between attributes. By contrast, starting from the second level, the Δ-closure levels include information about interactions between attributes. For example, the second level contains Δ-closed itemsets generated from passkeys of size 2. Datasets may have mainly three main types of behaviors. A majority of the datasets which were analyzed, namely "auto", "cylinder bands", "dermatology", "ecoli", "glass", "heart diseases", "hepatitis", "horse colic", "pen digit", "soybean", "wine", "zoo", have a closure structure which is similar to the closure structure of "adult", as displayed in Fig. 3a. We call this type of behavior "typical". For all these datasets, the ratio of the equivalence classes having passkeys with a low Δ increases exponentially with the level number. It means that most of the closed itemsets are not stable and can easily become not closed if new objects are added. The equivalence classes which are the most stable to noise are located at the first levels.

B. Visualizing the closure structure

The second type of behavior is given in Fig. 3b and termed as "deterministic". It is observed for "car evaluation" and "nursery" datasets, where each level mostly contains Δequivalence classes with passkeys having the same value. Actually these datasets were built from hierarchical decision models. Then each entry in such a dataset usually has a rather "determinististic" behavior, i.e., a much less random nature while higher order interactions between attributes are unusual.

The "intermediate" behavior can be different. In some cases the number of "unstable" Δ-equivalence classes grows less rapidly than for the majority of "typical" datasets. This means that a larger number of "stable" Δ-equivalence classes are present in the upper levels. The other datasets with similar behavior are "iris", "led7", "mushroom", and "tic tac toe".

Another type of intermediate behavior was observed for the datasets "ionosphere", "page blocks", and "pima", and is presented in Fig. 3d, where the Δ-closure structure of "pima" is shown. For these datasets, the ratio of the Δ-equivalence classes including passkeys with larger Δ-values may increase with levels, i.e., as the levels increase, the ratio of more robust Δ-equivalence classes increases as well. Intuitively, this kind of behavior is related to datasets where the 1st level of the Δclosure structure contains a substantial number of attributes of very low frequency. This induces unreliable or unstable itemsets in the upper levels of the Δ-closure structure. Finally, the patterns with high Δ-values are even more interesting if they are located at high levels of the closure structure, since they reflect high order interactions between attributes and are very stable to noise.

V. DISCUSSION AND CONCLUSION

In this paper we introduced the "Δ-closure structure" of a dataset, along with their theoretical and practical properties. The closure structure provides a view of the content of a dataset and may be used for guiding data mining. We have found three main types of "dataset behaviors" which are showing the stability and robustness of the corresponding closed itemsets. Here after we synthesize the main aspects of the closure structures:

• The Δ-closure structure is a partition of the whole set of itemsets. It is determined by the Δ-closed itemsets and their equivalence classes. The latter are represented by delta-passkeys, i.e., i.e., the smallest keys in an equivalence class. Accordingly, passkeys can represent any closed itemset of attributes without loss. • The closure structure is organized around a set of levels depending on the size of the passkeys, which is related to the distribution of the closed itemsets and their passkeys. Moreover, the Δ-closure structure allows us to study the stability and robustness of the closed itemsets. In future work, further studies are needed for a deeper understanding of the benefits of the Δ-closure structure for dealing with noisy data. More specifically, we should investigate the advantages brought by the Δ-closure structure for taking noise into account (i) in practical applications, and (ii) in interpreting the resulting patterns.

  Given a non Δ-closed itemset B, i.e., ∃m ∈ M (|B | -|(B ∪ {m}) | < ∆), B can be closed by iteratively changing B to B ∪ {m}, for any m violating (3) until no such attribute is found. The result is the Δ-closure of B. The corresponding closure operator is denoted by φ ∆ . Moreover any closure operator is associated with a class of equivalence.

  We denote the set of Δ-keys of B, or Δ-key set, by Key ∆ (B). Definition III.4. An itemset X ∈ Key ∆ (B) is called a Δpasskey if it has the smallest size among all keys in Key ∆ (B). For a Δ-closed itemset B, the Δ-passkey set, denoted by pKey ∆ (B), is given by pKey ∆ (B) = {X | X ∈ Key ∆ (B), |X| = min Y ∈Key∆(B) |Y |}, and it verifies pKey ∆ (B) ⊆ Key ∆ (B).

Fig. 2 :

 2 Fig. 2: The hierarchical structure of Δ-equivalence classes.

2 foreach c ∈ {c i } do 3 ∆ 4 Cls 10 while

 23410 pk (c) = 1; hasU pdates is T rue do 11 hasU pdates = F alse; 12 foreach c ∈ {c i } do 13 cls = Cls(c); 14 if |Ext(Ref (cls))| -|Ext(cls)| < d then 15 Cls(c) = Cls(cls); 16 hasU pdates = T rue; {c i } do 21 cls = Cls(c); 22 level ∆ (cls) = min(level ∆ (cls), level(c)); 23 end 24 foreach c ∈ {c i } do 25 cls = Cls(c); 26 if level(c) == level ∆ (cls) then 27 ∆ pk (c) = d; ], and easier to interpret, since Δ-passkeys are the smallest elements in the Δ-equivalence classes, (ii) the Δ-equivalence classes where passkeys have high Δ

Fig. 3

 3 Fig. 3 visualizes Δ-closure-structures for different types of data. The first level at the bottom of each subfigure shows the distribution of passkeys constituted of only one attribute w.r.t. their values of Δ. Actually this first level corresponds to

Fig. 3 :

 3 Fig. 3: The distribution of Δ-passkeys by levels based on 10 frequency bins distributed within the Δ-closure levels. The horizontal bars represent closure levels, where the level number is given in parentheses. The leftmost value is the percentage of closed itemsets |C k |/|C| at level k. The width of a color within a bar is proportional to the ratio of Δ-passkeys for Δ in (v 1 , v 2 ] at level C k .

TABLE I :

 I Dataset characteristics and computation time

	Datasets	|G|	|M |	# closed	runtime, sec GDPM Δ-passkeys
	adult	48842	95	359141	984	247
	chess kr k	28056	40	84636	146	52
	cyl. bands	540	120	39829537	2404	8803
	horse col.	368	81	173866	11	21
	ionosphere	351	155	23202541	2467	3090
	mushroom	8124	88	181945	164	27
	nursery	12960	27	115200	46	24
	pen digits	10992	76	3605507	12863	5020
	soybean	683	99	2874252	379	393