Anne-Laure Fougeres

Cecile Mercadier
email: mercadier@math.univ-lyon1.fr

Risk measures and multivariate extensions of Breiman's Theorem

Keywords: Discrete time model, Ruin probability, Value-at-Risk, Multivariate regular variation, Dependent risks

Modeling insurance risks is a task that received an increasing attention because of Solvency Capital Requirements. The ruin probability has become a standard risk measure to assess regulatory capital. In this paper we focus on discrete time models for finite time horizon. Several results are available in the literature allowing to calibrate the ruin probability by means of the sum of the tail probabilities of individual claim amounts. The aim of this work is to obtain asymptotics for such probabilities under multivariate regularly variation and, more precisely, to derive them from Breiman's Theorem extensions. We thus exhibit new situations where the ruin probability admits computable equivalents. Consequences are also derived in terms of the Value-at-Risk.

Introduction

Let x be a nonnegative real number, (ρ i ) i∈N * and (X i ) i∈N * be two sequences of random variables (r.v.), and define the sequence (R i ) i∈N by the recursive equation

R 0 = x, R i = R i-1 (1 + ρ i ) -X i , i ∈ N * .
(1)

This model has been used to model insurance risks. In such a context, x stands for the initial capital, ρ i ∈ (-1, ∞) for the random interest rate for the ith period and X i for the total claim amount minus the total premium incomes, so that R i represents the discounted surplus computed from period 0 to i. In this paper we focus on discrete time models for fixed finite time horizon d.

The ruin probability has become a standard risk measure to assess regulatory capital in insurance. For the model (1), the ruin probability within finite time horizon [0, d] and initial capital reserve x is defined by

ψ(x, d) = P min 1≤k≤d R k < 0 | R 0 = x .
Introducing the so-called discount factor from period i to 0, denoted by Y i = i j=1 (1 + ρ j ) -1 , one can re-express the discrete time risk model as

R 0 = x, R i = Y -1 i   x - i j=1 X j Y j   , i ∈ N * .
(2)

Hence, the ruin probability can be written as

ψ(x, d) = P max 1≤k≤d k i=1 X i Y i > x . (3)
Several results exist on the limiting behavior of ψ(x, d), especially in the case of an initial reserve x tending to infinity and when the distribution function of the X ′ i s is subexponential and sometimes heavy tailed. Most of them are stated under independence between the components of X = (X 1 , . . . , X d ): When Y = (Y 1 , . . . , Y d ) is deterministic, we know these asymptotics from [START_REF] Embrechts | Estimates for the probability of ruin with special emphasis on the possibility of large claims[END_REF], [START_REF] Sgibnev | On the distribution of the maxima of partial sums[END_REF], [START_REF] Embrechts | Modelling extremal events[END_REF], [START_REF] Ng | Maxima of sums of heavy-tailed random variables[END_REF] and [START_REF] Zhu | The uniform approximation of the tail probability of the randomly weighted sums of subexponential random variables[END_REF] among others. The case where the Y ′ i s are bounded r.v. has been studied by Tang and Tsitsiashvili [2003a,b]. Models governed by a specific dependence structure of the random vector Y have been proposed by [START_REF] Nyrhinen | On the ruin probabilities in a general economic environment[END_REF], Cai [2002] and [START_REF] Chen | Finite time ruin probability with heavy-tailed insurance and financial risks[END_REF]. Recent extensions to any dependence structure of Y can be found in [START_REF] Goovaerts | The tail probability of discounted sums of Pareto-like losses in insurance[END_REF], [START_REF] Wang | Uniform estimate for maximum of randomly weighted sums with applications to insurance risk theory[END_REF] and [START_REF] Wang | Tail probabilities of randomly weighted sums of random variables with dominated variation[END_REF] where the results are proven under a moment condition.

In some situations, the assumption of independence between the claim amounts X i might be unrealistic. [START_REF] Cossette | The discrete-time risk model with correlated classes of business[END_REF] considered special models of dependency. Allowing general dependence among the X ′ i s but keeping them independent from the nonnegative weights Y ′ i s, [START_REF] Zhang | Approximation of the tail probability of randomly weighted sums and applications[END_REF] derived approximations for ψ(x, d) under the assumptions that the X ′ i s have (extended) regularly varying tail and are asymptotically independent. Roughly speaking, the latter notion means that the joint upper tail of two claim amounts is negligible compared with each univariate tail. More precisely, [START_REF] Zhang | Approximation of the tail probability of randomly weighted sums and applications[END_REF] established the following equivalences, as x tends to infinity:

ψ(x, d) ∼ P d i=1 X i Y i > x ∼ d i=1 P (X i Y i > x) .
(4)

In their context, the first equivalence in ( 4) is a consequence of the assumption that the left tails of the X ′ i s are lighter than their right tails (see ( 9) below for a precise definition) and of the positivity of the weights Y i . The second equivalence in (4), which follows from the asymptotic independence of the X ′ i s, has also its own interest. Indeed, several papers are concerned with how the tail of the marginal distribution of an individual summand influences the asymptotic behavior of the sum. [START_REF] Barbe | On the tail behaviour of sums of dependent risks[END_REF] extended the equivalence between the tail of a sum and the sum of the tails, already obtained by [START_REF] Wüthrich | Asymptotic value-at-risk estimates for sums of dependent random variables[END_REF], [START_REF] Alink | Diversification of aggregate dependent risks[END_REF][START_REF] Alink | Analysis of the expected shortfall of aggregate dependent risks[END_REF], to a broader class of dependence structures using multivariate extreme value theory. Recently, [START_REF] Kortschak | Asymptotic results for the sum of dependent non-identically distributed random variables[END_REF] treated the case of non-identically distributed and not necessarily positive random variables.

The aim of this paper is to establish, in various contexts of dependence, asymptotics of the ruin probability and the Value-at-Risk. Our generalizations of (4) are in two directions: we allow dependence between the claim amounts X i and the weights Y i ; and we relax the assumption of asymptotic independence of the claim amounts. These results are obtained under the key assumption that X is Multivariate Regularly Varying (MRV) at infinity, see for instance Resnick [2007, chapter 6]. More precisely, the paper highlights the fact that our extensions are valid as soon as the random vector (X 1 Y 1 , . . . , X d Y d ) is also MRV, which is known as a Breiman's type result.

The rest of the paper is organized as follows: Section 2 contains notation and preliminary results. Several generalizations of Breiman's Theorem are stated in Section 3. Section 4 takes benefit of these extensions to derive asymptotics of risk measures as ruin probability and Value-at-Risk. A brief conclusion is given in Section 5. The proofs are postponed until Section 6.

Notation, definitions and preliminary results

We start this section with some notation. Let x • y define the componentwise product of two vectors

x = (x 1 , . . . , x d ) T and y = (y 1 , . . . , y d ) T ∈ R d , i.e. x • y = (x 1 y 1 , . . . , x d y d ) T and let z -1 = (z -1 1 , . . . , z -1 d ) T be the componentwise inverse of z ∈ [0, ∞] d . We define the inverse image of a set A for y ∈ [0, ∞] d by y -1 • A = {x ∈ R d | y • x ∈ A} = {(y -1 1 x 1 , . . . , y -1 d x d ) T , x ∈ A} .
Let X be a d-dimensional random vector. We say that X is multivariate regularly varying if there exists a non-null Radon measure ν on B d , the Borel σ-field of

E d = [-∞, ∞] d \ {0}, such that ν((-∞, ∞) d \ {0}) > 0 and a normalizing function a (with a(t) → ∞) such that tP(X/a(t) ∈ •) v → ν , (5)
as t tends to infinity, where v → refers to vague convergence on B d . Recall that a Radon measure on B d is a measure that is finite on each compact set of E d , and that a set A ⊂ E d is relatively compact if it is bounded away from zero. Recall also that a sequence of measures ν n converges vaguely to ν on

B d if E d f dν n converges to E d f dν for any function f compactly supported on E d , or equivalently if ν n (K) converges to ν(K) for each relatively compact set K of E d such that ν(∂K) = 0.
The limit measure ν is necessarily homogeneous, i.e. ν(tK) = t -α ν(K) for some α > 0 and all relatively compact Borel set K of E d . The function a is regularly varying with index 1/α. We shall write X ∈ MRV(α, a, ν) if (5) holds or simply MRV if no confusion can arise.

In the convergence (5) one may choose the normalizing function a such that lim

t→∞ t P( X > a(t)) = 1 , (6)
where • denotes some norm on R d . Given the choice of a norm on R d , a polar representation of the measure ν can be obtained. This result is due to de Haan and [START_REF] De Haan | Limit theory for multivariate sample extremes[END_REF], refer for instance to Resnick [2007, Theorem 6.1 p. 173 and Section 6.5.5 p. 201]. Let S d-1

• denote the unit sphere of R d relatively to the norm • . Then X is MRV if and only if there exists a measure H • (•) on S d-1

• , a positive real α, and a normalizing function a (with

a(t) → ∞) such that tP X a(t) , X X ∈ • v → (η α × H • )(•) , (7)
as t tends to infinity, where vague convergence holds on (0, ∞] × S d-1

• and η α denotes the Radon measure on (0, ∞] defined by η α (x, ∞] = x -α . The measure H • is called the spectral measure, or angular measure, and the choice of the normalizing function in (6) implies that H • is a probability measure on S d-1

• . The link between the limit measure ν and the spectral measure H • can be explicited via the following decomposition, for any A ∈ B d ,

ν(A) = T (A) αdr r α+1 H • (dw) , (8)
where T (z) = ( z , z/ z ) for any z ∈ E d and T (A) is the usual image of A by T .

We will use the following property in the sequel, and sketch the proof in Section 5.

Lemma 1. If X is MRV(α, a, ν) on R d and if the left tail of each component X i is lighter than its right tail, i.e.

lim x→∞ P(X i < -x) P(X i > x) = 0 , 1 ≤ i ≤ d , (9)
then the measure ν is concentrated on [0, ∞] d \ {0}, the spectral measure

H • of X is con- centrated on S d-1 • ∩ [0, ∞] d and lim x→∞ P(max 1≤j≤d j i=1 X i > x) P( d i=1 X i > x) = lim x→∞ P( d i=1 X + i > x) P( d i=1 X i > x) = 1 .
Remark 1. Following Kortschak and Albrecher [2009, Lemma 3.1] the condition (9) may be replaced by

P(X i > a, X j > b) ≥ P(X i > a)P(X j > b) for all(a, b) ∈ R 2 and 1 ≤ i < j < n. ( 10
)
Throughout this paper, we will choose the ℓ 1 -norm • 1 defined by

x 1 = d j=1 |x i | ,
and we will denote S d-1

1 and H 1 the corresponding unit sphere and spectral measure. The motivation of these choices for the normalizing function and the norm on R d is that if X is MRV and satisfies the assumptions of Lemma 1, then on one hand, one has for all x > 0, lim

t→∞ tP d i=1 X i > a(t)x = x -α . ( 11
)
For any d-dimensional vector X, let Q(X) denote the limit, when it exists,

Q(X) = lim x→∞ P( d i=1 X i > x) d j=1 P(X j > x)
.

Similar limits to Q(X) have been considered in the literature. The following result can be found in several references and under different assumptions. We refer for instance to [START_REF] Alink | Diversification of aggregate dependent risks[END_REF], [START_REF] Barbe | On the tail behaviour of sums of dependent risks[END_REF], [START_REF] Kortschak | Asymptotic results for the sum of dependent non-identically distributed random variables[END_REF], Embrechts et al. [2009b,c].

The unique assumption needed is multivariate regular variation, as proven in Section 6. The components of a regularly varying random vector are called asymptotically independent if its spectral measure (with respect to any norm) is concentrated on the axes.

Lemma 2. If X is a d-dimensional multivariate regularly varying random vector with index α > 0 such that (9) holds, then Q(X) exists and

d -(1-α) + ≤ Q(X) ≤ d (α-1) + . ( 12
)
If the components of X are asymptotically independent, then Q(X) = 1.

Remark 2. Note that these bounds are universal and (obviously) does not depend on any particular choice of norm on R d , even though the proof we give below makes use of the ℓ 1norm. For analogous comments on the choice of norm, we refer to the remark in [START_REF] Mainik | On optimal portfolio diversi cation with respect to extreme risks[END_REF] after the representation (13) therein or to [Embrechts et al., 2009a, Proposition 4.6].

Remark 3. Note that contrary to [START_REF] Barbe | On the tail behaviour of sums of dependent risks[END_REF] and [Embrechts et al., 2009a, Corollary 4.2], we do not assume that the random variables X i are positive nor that they have the same marginal distribution. This implies in particular that the usual standardization w i H 1 (dw) = 1/d (Barbe et al. [2006, Equation (11)] or Beirlant et al. [2004, p. 260]) does not hold here, but is not needed. [START_REF] Kortschak | Asymptotic results for the sum of dependent non-identically distributed random variables[END_REF] also considered such generalizations to non positive or non identically distributed marginals. However, the bounds given in Lemma ( 12) are new in this context.

Remark 4. If X is multivariate regularly varying with index α > 0, for any fixed norm on R d , Q(X) depends only on α and on the associated spectral measure H. Thus, with an abuse of notation, we can denote Q(X) = Q(α, H). The properties of Q(α, H) have first been investigated by [Barbe et al., 2006, Proposition 2.2]. For a given α, the upper and lower bounds in ( 12) are achieved by independent or fully dependent components, and for a given spectral measure H, Q(α, H) is increasing in α. See also [START_REF] Kortschak | Asymptotic results for the sum of dependent non-identically distributed random variables[END_REF] and Embrechts et al. [2009b,c].

Multivariate extensions of Breiman's Theorem

Let X and Y be two random vectors in R d . A Breiman's type result consists in obtaining sufficient conditions for the vector Y • X to be MRV. We start by recalling a particular case of Basrak et al. [2002, Proposition A.1] who prove the multivariate regular variation of MX where M is a random matrix independent of X. Denote by ν Y the measure defined on E d by

ν Y (K) = E[ν(Y -1 • K)] = [-∞,∞] d ν(y -1 • K) P Y (dy) . (13)
Theorem 3. Let X ∈ MRV(α, a, ν). Let Y be a random vector independent of X. Assume that there exists a positive ε such that 0

< E[|Y i | α+ǫ ] < ∞ for each i in {1, . . . , d}. Then the random vector Y • X ∈ MRV(α, a, ν Y ), i.e. tP(Y • X ∈ a(t)•) v → ν Y ,
as t tends to infinity.

The hypothesis of independence between X and Y might be restrictive. In the sequel, we generalize Theorem 3 on this point. The rest of the section is divided according to the type of dependence considered between X and Y. We investigate three situations that seem meaningful in an actuarial context:

• X is MRV and asymptotically independent of Y;

• X has independent and identically distributed (i.i.d.) regularly varying components and Y is predictable with respect to X;

• X, Y are jointly MRV and asymptotically dependent.

In each of these situations, an extension of Breiman's Theorem is stated.

Case of asymptotic independence

The first generalization of Theorem 3 is done under the condition that X and Y are asymptotically independent in the following sense. We assume that

tP X a(t) , Y ∈ • v → (ν × L)(•) (14)
on the Borel sets of

E d × [-∞, ∞] d
where ν is a Radon measure on E d not concentrated at infinity and L is a probability measure on [-∞, ∞] d . Note that the roles of X and Y are not symmetric in this definition since it specifies their independence when X is large only. Condition ( 14) implies in particular that X ∈ MRV(α, a, ν) for some positive α. This case obviously contains the case of stochastic independence between X and Y, L being then the distribution of Y.

Following [START_REF] Maulik | Asymptotic independent and a network traffic model[END_REF], we make the following asymptotic negligibility assumption. Assume that for some δ > 0 and any i in {1, . . . , d}

lim ε→0 lim sup t→∞ tE |X||Y i | a(t) δ ✶ |X|/a(t)≤ε = 0 , (15)
and also that

[-∞,∞] d y α L(dy) < ∞ . ( 16
)
Let ν L be the measure defined on E d by

ν L (A) = (ν × L) {(x, y), x ∈ y -1 • A} = [-∞,∞] d ν(y -1 • A)L(dy) . Note that if Y * is a random vector with distribution L, then ν L = ν Y * .
Theorem 4. Assume that Assumptions ( 14), ( 15) and ( 16) hold. Then the random vector

Y • X ∈ MRV(α, a, ν L ), i.e. tP(Y • X ∈ a(t)•) v → ν L ,
as t tends to infinity.

The proofs of Theorem 4 and subsequent results are postponed until Section 6. We establish in the proof of Theorem 4 a more general result on a random vector MX, where M is a random matrix of size q × d.

Case of predictable weights

The following generalization uses the predictable framework introduced by Hult and Samorodnitsky [2008]. In this context, we assume that the components of the random vector X are independent (non necessarily identically distributed) and regularly varying at infinity. Then the vector X is multivariate regularly varying. More precisely, the sequence tP(X/a(t) ∈ •) converges vaguely to a measure ν, which is concentrated on the axes

ν = d j=1 ν j , where ν j (A) = ν(A ∩ δ j ) and δ j = {z ∈ R d | z j = 0, z i = 0, i = j} is the j-th punctured coordinate axis.
The predictable framework of [START_REF] Hult | Tail probabilities for infinite series of regularly varying random vectors[END_REF] consists in assuming that there exists a filtration with respect to which the X ′ j s are measurable and the Y ′ j s are predictable. This implies in particular that for each j, X j and Y j are independent. This framework is of interest in time series. An example is the EGARCH process of [START_REF] Nelson | Conditional heteroskedasticity in asset returns: a new approach[END_REF]. This process, say {ζ j }, can be expressed as

ζ j = X j Y j , Y j = exp ∞ i=1 c i η j-i ,
and the relevant filtration is F j = σ(X i , η i , i ≤ j -1), {X j } and {η j } are two i.i.d. sequences, not necessarily independent of each other, E[η j ] = 0, Var(η j ) = 1 and ∞ i=1 c 2 i < ∞. The process {Y j } is called the volatility. The extremal properties of this process have been studied by [START_REF] Davis | Point process convergence of stochastic volatility processes with application to sample autocorrelation[END_REF] in the case of independence between the processes {X j } and {Y j }.

Theorem 5. Let X and Y be two random vectors of R d . Assume that the components of X are independent (non necessarily identically distributed) and X is MRV(α, a, ν). Suppose that there exists a filtration {F j } such that X j is F j+1 -measurable and Y j is F j -measurable. Assume moreover that there exists ε > 0 such that

E[|Y j | α+ε ] < ∞ ∀j = 1, . . . , d . (17)
Then the random vector Y • X is multivariate regularly varying and

tP(Y • X ∈ a(t)•) v → d j=1 E[Y α j ]ν j , (18)
as t tends to infinity.

Remark 5. The assumptions of the previous theorem are weaker than those of [START_REF] Hult | Tail probabilities for infinite series of regularly varying random vectors[END_REF]. Their mean condition (3.1) and their assumptions (3.7)-(3.9) here reduce to (17) thanks to the finite time horizon d.

Case of joint multivariate regular variation

When neither independence nor asymptotic independence is a relevant assumption, one might be interested in extensions of Breiman's Theorem under asymptotic dependence. Specifically, we assume that X and Y are jointly multivariate regularly varying, i.e. there exist a and b such that

tP X i a(t) , Y i b(t) i=1,...,d ∈ • v → ν X,Y , (19)
as t tends to infinity, for a non-null Radon measure ν X,Y on E 2d . We assume that the measure ν X,Y • Π -1 is not identically zero. It means that there exists one index i such that the pair (X i , Y i ) is asymptotically dependent, i.e.

∃ i ∈ {1, . . . , d} : ν

(i) X,Y ((0, ∞] 2 ) > 0 , (20)
where ν

(i)
X,Y is defined as the restriction of ν to the i-th coordinates in x and y. For instance

when i = 1, ν (1) (A) = ν 2d (A×[-∞, ∞] 2d-2 ) for any A ∈ B 2 .
The normalizing functions a and b are regularly varying with indices respectively denoted by 1/α and 1/β for some positive real numbers α and β.

Define the map Π

: R d × R d → R d by Π(x, y) = y • x.
Theorem 6. Assume that Assumptions ( 19) and ( 20) hold. Then the random vector Y

• X ∈ MRV(αβ/(α + β), ab, ν X,Y • Π -1 ), i.e. tP(Y • X ∈ a(t)b(t)•) v → ν X,Y • Π -1 ,
as t tends to infinity.

In the particular case Y 1 = • • • = Y d , this result can be found in Resnick [2007, Proposition 7.6].

Remark 6. Traditional models include an assumption of regular variation on the exp(Y i )'s instead of the Y i 's; for instance when the Y i 's are log-normal. As a consequence, assuming that the discount factors have heavy tails might be considered as a strong hypothesis in a usual context. However, this model yields to quite light tails for large values of β and allows dependence between the X and the Y components.

Application to risk measures

The link between a Breiman's type result and asymptotics for the ruin probability is highlighted in the first part of this section. The equivalents that we obtain depend on the ratio Q(Y • X) which is not always explicitly known. Since numerical approximations can be done, we show some patterns in the second part of this section in the case of the asymmetric logistic measure. Finally, in the third part of this section, we derive the consequences of our results on another risk measure: the Value-at-Risk.

Ruin probability

Recall that the ruin probability associated to the model ( 1) is defined by

ψ(x, d) = P max 1≤k≤d k i=1 X i Y i > x
and that for any random vector Z we have defined Q(Z) by

Q(Z) = lim x→∞ P d i=1 Z i > x d i=1 P(Z i > x)
when the limit exists.

The next result establishes that the first equivalence of (4) holds if the random vector (X 1 Y 1 , . . . , X d Y d ) is MRV and if all the random variables X i have lower tails lighter than their upper tails.

Theorem 7. Under the assumptions of any one of Theorems 3, 4, 5 or 6, if the components of Y are nonnegative and (9) holds, then one has, when x tends to infinity,

ψ(x, d) ∼ P d i=1 X i Y i > x . (21)
The second equivalence of ( 4) is specific to the asymptotic independence context considered by [START_REF] Zhang | Approximation of the tail probability of randomly weighted sums and applications[END_REF] and does not hold in general. However, as soon as a Breiman's type result is valid, one can derive other asymptotics for (21) depending explicitly on the margins of X and on Q(Y • X). We make these relations explicit in the following corollaries, which proofs are deferred to Section 6. Let us introduce the two sets

T = {z ∈ Rd | d i=1 z i > 1} and T i = {z ∈ Rd | z i > 1}.
Corollary 8. Under the assumptions of Theorem 3, if the components of Y are nonnegative and (9) holds, then when x tends to infinity,

ψ(x, d) ∼ Q(Y • X) d i=1 E[Y α i ] P(X i > x) , (22)
where

Q(Y • X) = ν Y (T ) d i=1 ν Y (T i ) = E[ν(Y -1 • T )] d i=1 E[Y α i ]ν(T i ) . ( 23
)
Note that for any deterministic vector Y, the numerator ν Y (T ) of the middle term in (23) corresponds to the extreme risk index of the portfolio denoted by γ Y . See [START_REF] Mainik | On optimal portfolio diversi cation with respect to extreme risks[END_REF] and references therein for complements on this risk measure.

Corollary 9. Under the assumptions of Theorem 4, if the components of Y are nonnegative and (9) holds, then when x tends to infinity,

ψ(x, d) ∼ Q(Y • X) d i=1 P(X i > x) ∞ 0 y α i dL(y) , (24)
where

Q(Y • X) = ν L (T ) d i=1 ν L (T i ) = ν(y -1 • T ) dL(y) d i=1 ν(T i ) y α i dL(y)
.

Remark 7. Corollary 8 recovers the results (for finite time horizon) of [START_REF] Zhang | Approximation of the tail probability of randomly weighted sums and applications[END_REF]. Indeed, under Corollary 8, as well as under Corollary 9, it can easily be proven that the components of Y • X are asymptotically independent as soon as those of X are asymptotically independent, so that Q(Y • X) = 1. This comes from linearity properties combined with the fact that for any vector y with positive components, any i ∈ {1, . . . , d} and any subset A of

E d one has (y -1 • A) ∩ Re i = y -1 • (A ∩ Re i ) .
Corollary 10. Under the assumptions of Theorem 5, if the components of Y are nonnegative and (9) holds, then when x tends to infinity,

ψ(x, d) ∼ d i=1 E[Y α i ] P(X i > x) . (25)
In the latter corollary, the independence of the claim amounts X i yields Q(Y • X) = 1, so that actuarial and financial informations are separated in (25).

Corollary 11. Under the assumptions of Theorem 6, one has, when x tends to infinity,

ψ(x, d) ∼ Q(Y • X) d i=1 P(X i Y i > x) , (26)
where

Q(Y • X) = ν X,Y • Π -1 (T ) d i=1 ν X,Y • Π -1 (T i )
.

Let us make some comments on the constant Q(Y • X). We know from Lemma 2 that under the assumptions of any of the previous corollaries, one has

d (1-δ) + ≤ Q(Y • X) ≤ d (δ-1) +
with δ = α under Corollary 8 and 9 and δ = αβ/(α + β) under Corollary 11. Consequently, when δ ∈ (0, 1] the control is accurate. If δ > 1, these bounds fail to be sharp as the dimension d grows up. However, it is possible to derive bounds depending on moments of Y. Let Y (1) and Y (d) be the smaller and greater components of Y. For instance, Corollary 8 yields

d i=1 ν(T i ) 1 E[Y α (1) ] d i=1 E[Y α i ]ν(T i ) Q(X) ≤ Q(Y • X) ≤ d i=1 ν(T i ) 1 E[Y α (d) ] d i=1 E[Y α i ]ν(T i ) Q(X) .
These bounds are particularly interesting in the fact that they allow to separate the financial information contained in Y and the actuarial part given by X and ν. The previous upper and lower bounds can also be relaxed and lead to

E[Y α (1) ] max i=1,...,d E[Y α i ] Q(X) ≤ Q(Y • X) ≤ E[Y α (d) ] min i=1,...,d E[Y α i ] Q(X) .
Analogous bounds may be obtained under the assumptions of Corollary 9:

[0,∞] d y α (1) dL(y) max i=1,...,d [0,∞] d y α i dL(y) Q(X) ≤ Q(Y • X) ≤ [0,∞] d y α (d) dL(y) min i=1,...,d [0,∞] d y α i dL(y) Q(X) .
When the limit measure ν of X is given, numerical bounds for Q(Y • X) can thus be obtained from the preceding inequalities. This relies on the computation of the value of Q(X), that is illustrated in the following section for the bivariate logistic dependence case.

Illustrative features of Q(X)

In this part, we illustrate the behaviour of the term Q(X) in the bivariate asymmetric logistic case. For other examples, refer e.g. to [START_REF] Alink | Diversification of aggregate dependent risks[END_REF], [START_REF] Barbe | On the tail behaviour of sums of dependent risks[END_REF], [START_REF] Kortschak | Asymptotic results for the sum of dependent non-identically distributed random variables[END_REF], Embrechts et al. [2009b,c]. Recall that when d = 2, Q(X) may be written in terms of the limit measure ν as follows:

Q(X) = ν(T ) ν(T 1 ) + ν(T 2 )
,

where

T = {z ∈ R2 | z 1 + z 2 > 1} and T i = {z ∈ R2 | z i > 1}.
We consider a random vector X = (X 1 , X 2 ) with the following bivariate extreme value distribution: -both margins are identically distributed from the Fréchet(µ = 0, σ = 1, ξ = α) distribution, which means that there exists some positive α such that for any positive x one has F X i (x) = P(X i ≤ x) = exp (-x -α ) ; -the dependence between the margins of X is characterized by a function ℓ satisfying the conditions of a stable tail dependence function (see e.g. L1, L2 and L3 conditions of [Beirlant et al., 2004, page 257]), so that the distribution function of X can be written as

F X (x) = P (X ≤ x) = exp -ℓ x -α 1 , x -α 2 .
We fix the dependence tail function to be the asymmetric logistic defined by

ℓ ψ 1 ,ψ 2 ,r (x 1 , x 2 ) = (1 -ψ 1 )x 1 + (1 -ψ 2 )x 2 + (ψ 1 x 1 ) 1/r + (ψ 2 x 2 ) 1/r r ,
for any 0 < r ≤ 1 and 0 ≤ ψ 1 , ψ 2 ≤ 1. This parametric model has been widely used, and is e.g. presented in Section 9.2.2 of [START_REF] Beirlant | Ruin probabilities with dependent rates of interest[END_REF].

We compute Q(X) by standard upper and lower Riemann approximations. In dimension two, this procedure is very accurate: lower and upper curves are so closed that they are indistinguishable, see Figures 1 and2.

In the special case where ψ 1 = ψ 2 = 1, the distribution of X is the bivariate symmetric logistic model with Fréchet margins. The strength of dependence between the components of X is a decreasing function of r. In particular, the independence (resp. the total positive dependence) corresponds to r = 1 (resp. r → 0). The values of Q(X) for this model are given on Figure 1 as functions of r and α. Several aspects of these plots were already observed in the literature for other models: Q(X) tends to 1 as r tends to 1; Q(X) tends to 2 α-1 as r tends to 0; For α < 1 (resp. α > 1), Q(X) is strictly increasing (resp. decreasing) in r; For all r, Q(X) is strictly increasing in α.

ψ ψ 1 = = 1, , ψ ψ 2 = = 1 α α Q( ( X
More generally, the patterns of Q(X) when ψ 1 and ψ 2 are in (0, 1) are presented in Figure2. For specific values of α, ψ 1 and r, we plot Q(X) as a function of the parameter ψ 2 . Again we observe several elements on these plots: Q(X) is strictly decreasing (resp. increasing) in ψ 1 or ψ 2 when α < 1 (resp. α > 1); For any values of ψ 1 and ψ 2 , Q(X) is strictly increasing (resp. decreasing) in r when α < 1 (resp. α > 1).

The computation of numerical approximations of Q(X) can also be performed for values of d larger than two. Note however that this procedure becomes quickly time consuming when d is large.

Value-at-Risk

Another classical risk measure is the Value-at-Risk, defined as the following: Given p ∈ (0, 1) and a random variable Z, the Value-at-Risk of Z at level p, denoted by VaR p (Z), is the (1 -p)-th quantile of Z, i.e.:

VaR p (Z) = inf{z | P(Z > z) ≤ p} .
An important feature of risk measures is the property of subadditivity (resp. superadditivity), which can be written for the Value-at-Risk in terms of the ratio

R (p) Z = VaR p d i=1 Z i d i=1 VaR p (Z i ) , (27)
where Z denotes the vector (Z 1 , . . . , Z d ) T . We say that subadditivity (resp. superadditivity)

holds for VaR p (Z) if R (p)
Z ≤ 1 (resp. > 1). Analogously, asymptotic subadditivity (resp. asymptotic superadditivity) holds for VaR p (Z) if:

lim p→0 R (p) Z ≤ 1 (resp. > 1) .
The following corollary summarizes the consequences of the asymptotics obtained in the previous subsection on the Value-at-Risk, and allows in particular to discuss explicit cases where asymptotic subadditivity (resp. asymptotic superadditivity) of VaR p (Y • X) holds. Such kind of comments are given by Embrechts et al. [2009a] in the case of elliptical claim amounts.

Corollary 12. Assume that the components of Y are nonnegative and (9) holds. Then, under the assumptions of Theorem 3,

lim p→0 R (p) Y•X = Q(Y • X) 1/α d i=1 E[Y α i ]ν(T i ) 1/α d i=1 (E[Y α i ]ν(T i )) 1/α .
Under the assumptions of Theorem 4,

lim p→0 R (p) Y•X = Q(Y • X) 1/α d i=1 ∞ 0 y α i L(dy)ν(T i ) 1/α d i=1 ∞ 0 y α i L(dy)ν(T i ) 1/α .
Under the assumptions of Theorem 5,

lim p→0 R (p) Y•X = d i=1 E[Y α i ]ν(T i ) 1/α d i=1 (E[Y α i ]ν(T i )) 1/α .
Under the assumptions of Theorem 6,

lim p→0 R (p) Y•X = Q(Y • X) 1/α d i=1 ν X,Y • Π -1 (T i ) 1/α d i=1 (ν X,Y • Π -1 (T i )) 1/α .
The proof of Corollary 12 is postponed until Section 6.

Concluding comments and discussion

In this paper, discrete time risk models with finite time horizon have been considered, and asymptotics for risk measures, as ruin probability or Value-at-Risk, have been obtained under different dependence settings for the claim amounts and discount factors. The fruitful role of the multivariate regularly varying setting has been highlighted, as well as the usefulness of Breiman's type results. This allowed to generalize (for finite time horizon) [START_REF] Zhang | Approximation of the tail probability of randomly weighted sums and applications[END_REF] result outside the asymptotic independence of the claim amounts, and outside the independence of claim amounts and discount factors.

A specific parameter Q(Y • X) arises from these asymptotics, for which explicit bounds have been provided in terms of the limit measure of X and the characteristics of the discount factors These bounds are easily numerically computable as far as the time horizon d is not too large. Dealing with very high dimensions still represents a challenging numerical problem. Ideas coming from the algorithms developed in [START_REF] Arbenz | The aep algorithm for the fast computation of the distribution of the sum of dependent random variables[END_REF] could be promising for this task.

Another important issue is to measure the accuracy of the approximations stated in Theorem 7 and its corollaries. Such a problem requires MRV hypotheses with second order conditions, as formulated in the univariate case by [START_REF] Degen | Risk concentration and diversification -second-order properties[END_REF]. This will be the scope of a future work.

Proofs

Proof of Lemma 1. We prove the first assertion of Lemma 1 by induction. Let K be a relatively compact set of E d-1 such that ν(∂K) = 0 and let ǫ > 0. By assumption, the convergence (5) implies

lim t→∞ tP(a(t) -1 X ∈ K × (-∞, -ǫ]) = ν(K × (-∞, -ǫ]) . Besides, tP(a(t) -1 X ∈ K × (-∞, -ǫ]) ≤ tP(X d < -a(t)ǫ) → 0 , therefore ν(K × (-∞, -ǫ]) = 0, so the support of ν is included in R d-1 × [0, ∞) \ {0}.
The rest of the induction argument is along the same lines. It is then obvious that the measure H • is concentrated on S d-1

• ∩ [0, ∞] d . Now for simplicity, we consider the following sets

T max = {z ∈ Rd | max 1≤k≤d k i=1 z i > 1} , T = {z ∈ Rd | d i=1 z i > 1} , T + = {z ∈ [0, ∞] d | d i=1 z i > 1} .
The last assertion of Lemma 1 may be then written as the following:

ν(T max ) ν(T ) = ν(T + ) ν(T ) = 1,
which is obvious when ν is concentrated on the nonnegative quadrant.

Proof of Lemma 2. Recall that H 1 denotes the spectral measure of X with respect to the ℓ 1 -norm. This is a probability measure concentrated on {w ∈ S d-1

1 | w i ≥ 0 , 1 ≤ i ≤ d}.
The choice of the ℓ 1 -norm yields that lim t→∞ tP( d i=1 X i > a(t)) = 1. Besides, one gets from (8), for any y ∈ [0, ∞] d , that

ν z ∈ [0, ∞] d \ {0} | d i=1 y i z i > 1 = S d-1 1 d i=1 y i w i α H 1 (dw) , which implies in particular that lim t→∞ tP(X i > a(t)) = S d-1 1 w α i H 1 (dw) .
As a consequence, one gets

lim t→∞ t d i=1 P(X i > a(t)) = d i=1 S d-1 1 w α i H 1 (dw) , so that Q(X) = 1 d i=1 S d-1 1 w α i H 1 (dw)
.

If α > 1, then the function x → x α is convex. By Jensen's inequality, we have

1 d d i=1 S d-1 1 w α i H 1 (dw) ≥ 1 d d i=1 S d-1 1 w i H 1 (dw) α = 1 d S d-1 1 d i=1 w i H 1 (dw) α = d -α .
If α < 1, then again by Jensen's inequality, the reverse bound holds:

d -α = 1 d d i=1 S d-1 1 w i H 1 (dw) α ≥ 1 d d i=1 S d-1 1 w α i H 1 (dw) .
On the other hand, if α > 1, then w α ≤ w for w ∈ [0, 1] so

d i=1 S d-1 1 w α i H 1 (dw) ≤ d i=1 S d-1 1 w i H 1 (dw) = 1 .
The reverse inequality obviously holds for α < 1. Gathering these bounds yields (12). The equality Q(X) = 1 when the components of X are asymptotically independent follows from the fact that H 1 is then concentrated on the axes.

Proof of Theorem 4

In order to simplify notations, we denote by • a given norm on an Euclidean space and for any q × d matrix M , the induced matrix norm is M = sup x =1 M x . For a q × d matrix M and a set K ⊂ R d , we define

M -1 • K = {x ∈ R d | M x ∈ K} .
In order to prove Theorem 4, we state and prove a more general result. We use the concept of asymptotic independence introduced by [START_REF] Maulik | Asymptotic independent and a network traffic model[END_REF]. Assume that M is a random matrix of size q × d and X ∈ R d is a random vector satisfying

X a(t) , M ∈ • v → (ν × G)(•) (28)
as t tends to infinity on E d × [-∞, ∞] qd , where ν is a Radon measure on E d not concentrated at infinity and G is a probability measure on [-∞, ∞] qd . This implies that ν is homogeneous with positive index α, say, and we can still choose the normalizing function a such that lim t→∞ tP( X > a(t)) = 1 .

Assume also that there exists δ > 0 such that lim ǫ→0 lim sup

t→∞ tE M X a(t) δ ✶ { X ≤ǫa(t)} = 0 , (29) [-∞,∞] qd M α G(dM ) < ∞ . ( 30
)
Let ν G denote the measure defined on

E d by ν G (K) = ν ⊗ G ({(x, M ) | M x ∈ K}) = [-∞,∞] qd ν(M -1 • K) G(dM ) = E[ν(M * -1 • K)] ,
where M * is a random matrix with distribution G.

Theorem 13. Assume that Assumptions (28), ( 29) and ( 30) hold. Then

tP(a(t) -1 MX ∈ •) v → ν G
as t tends to infinity.

Proof of Theorem 13. Let K be relatively compact in E d such that ν G (∂K) = 0. Fix some real number s > 0, write

P (MX ∈ a(t)K) = P (MX ∈ a(t)K , M ≤ s) + P ∈ a(t)K , M > s) .
Since K is relatively compact in E d , there exists ǫ > 0 such that x ≥ ǫ for all x ∈ K. Thus M ≤ s and M x ∈ K imply that x ≥ s -1 ǫ, thus Assumption (28) implies that

lim t→∞ tP (MX ∈ a(t)K , M ≤ s) = E ν(M * -1 • K)✶ { M ≤s} .
By homogeneity of ν, Condition (30) and the bounded convergence theorem, it holds that

lim s→∞ E[ν(M * -1 • K)✶ { M ≤s} ] = E[ν(M * -1 • K)] .
Next, since y ∈ K implies y ≥ ǫ, we have, for s ≥ 1,

P(MX ∈ a(t)K , M > s) ≤ P ( M X > ǫa(t) , M > s) ≤ P ( X > ǫa(t) , M s) + P M X a(t) > ǫ , X a(t) ≤ ǫ .
By Assumption (28),

lim t→∞ tP ( X > ǫa(t) , M > s) = ǫ -α P( M * > s) ,
which can be made arbitrarily small by choosing s large enough. By Markov inequality and Assumption (29

) lim sup t→∞ tP M X a(t) > ǫ , X a(t) ≤ ǫ ≤ ǫ -δ lim sup t→∞ tE M X a(t) δ ✶ { X ≤ǫa(t)} = 0 .
Since s can be chosen arbitrarily large, we have proven that lim t→∞ tP(MX ∈ a(t)K) = E[ν(M * -1 • K)] and this concludes the proof.

Proof of Theorem 5. It is a straightforward consequence of Hult and Samorodnitsky [2008, Theorem 3.1] for p = 1; A 0 = 0 and Z 0 = 0; A j = Y j e j and Z j = X j for j = 1, . . . , d; A j = 0 and Z j = 0 for j > d.

Proof of Theorem 6. Let K be a relatively compact set in E d such that ν X,Y • Π -1 (∂K) = 0. Then ν X,Y (∂Π -1 K) = 0, since Π is continuous, which implies that ∂Π -1 (K) ⊂ Π -1 (∂K). Moreover, Π -1 (K) is relatively compact in E 2d . To see this, we can choose some arbitrary norm and prove that for some ξ > 0 the set

K = {x ∈ R d | x > ξ} is such that Π -1 (K) is bounded away from zero in E 2d .
Choose for instance the euclidean norm in both E d and E 2d . Then

x • y 2 = d i=1 x 2 i y 2 i ≤ d i=1 (x 2 i + y 2 i )y 2 i ≤ d i=1 (x 2 i + y 2 i ) 2 = (x, y) 2 . Thus, if (x, y) ∈ Π -1 (K), then (x, y) ≥ x•y > ξ and this proves that Π -1 (K) is relatively compact in E 2d . Denote c(t) = a(t)b(t). Then, lim t→∞ tP(Y • X/c(t) ∈ K) = lim t→∞ tP X i a(t) , Y i b(t) i=1,...,d ∈ Π -1 (K) = ν X,Y • Π -1 (K) .
Proof of Theorem 7. We want to state (21). Under the assumptions of any of the theorems of Section 3, Y • X is multivariate regularly varying. Let µ denote the Radon measure associated to Y • X and a the normalizing function. Then ( 21) is equivalent to

µ(T ) = µ(T max ) = µ(T + ) .
This holds as soon as the support of µ is included in [0, ∞] d \ {0}, which we prove now. By Lemma 1, it suffices to prove that (9) holds for the vector Y • X, i.e. for each i = 1, . . . , d,

lim x→∞ P(X i Y i < -x) P(X i Y i > x) = 0 . (31)
The latter equivalence holds under the assumptions of Theorems 3, 4 and 5 since X i and Y i are asymptotically independent, so that

P(X i Y i > x) ∼ c i P(X i > x) and P(X i Y i < -x) ∼ c i P(X i < -x)
in each of these cases, where c i = E[Y α i ] under the assumptions of Theorems 3 and 5 and c i = ∞ 0 y α i L(dy) under the assumptions of Theorem 4. Under the assumptions of Theorem 6, we must prove that if X and Y are two jointly regularly random variables such that Y is non negative and X satisfies (9), then XY also satisfies (9). This clearly holds, since for any x > 0, 0 = lim t→∞ tP(X < -a(t)x) = ν X,Y ((-∞, -x) × [0, ∞)) .

Thus the support of ν X,Y is included in [0, ∞] 2 . This proves (31) and concludes the proof of ( 21) under the assumptions of Theorem 6, and therefore also the proof of Theorem 7.

Proof of Corollaries 8 -9 -10 -11. We need to prove respectively ( 22) -( 24) -( 25) and ( 26). Note that (26) follows directly from the definition of Q(X • Y) and is valid in all contexts. Thus, we only need to give further equivalents of d i=1 P(X i Y i > x).

• Under the assumptions of Theorems 3 and 5, X i and Y i are independent, and Breiman's Theorem applies: P(X i Y i > x) ∼ E[Y α i ] P(X i > x). Under the assumptions of Theorem 5, the components of Y • X are asymptotically independent as shown by ( 18), so that a consequence of Lemma 2 is that Q(Y • X) = 1.

• Under the assumptions of Theorem 4, it follows from conditions ( 14), ( 15) and ( 16) that P(X i Y i > x) ∼ ∞ 0 y α i L(dy)P(X i > x).

Proof of Corollary 12. The following lemma is useful to prove Corollary 12, and is a straightforward consequence of usual properties of inverses of regularly varying functions, refer e.g. to Resnick [1987, Proposition 0.8 (vi)]. We state it here (without proof) for completeness. Under the assumptions of any of the theorems 3, 4, 5 or 6, we want to obtain the following common expression:

lim p→0 VaR p d i=1 X i Y i d i=1 VaR p (X i Y i ) = {Q(Y • X)D} 1/α , ( 32
)
where D is defined by This limit D can be explicited using the change of variable x = γ(t), where γ(t) = a(t) under the assumptions of Theorems 3-4-5, and γ(t) = a(t)b(t) under the assumptions of Theorem 6. One gets then, under any of these assumptions:

D = d i=1 µ(T i ) d i=1 µ(T i ) 1/α α , (33)
lim t→∞ tP(X i Y i > γ(t)) = µ(T i ) . ( 34
)
This gives the numerator announced in (33). The denominator is obtained applying twice Lemma 14. Indeed, we get first from (34) that, for any i, j ∈ {1, . . . , d}, which is the expected denominator in (33). Checking that µ(T j ) has the different expressions announced under the assumptions of Theorems 3-4-5-6 is straightforward and ends the proof of Corollary 12.

Figure 1 :

 1 Figure1: Values of Q(X) as a function of r or α when X has the symmetric logistic dependence structure ℓ 1,1,r and Fréchet(α) margins. Left plot gives Q(X) as a function of the dependence parameter r ∈ (0, 1), for values of α among {0.1, 0.8, 1.5, 2.2, 3}, and right plot gives Q(X) as a function of α ∈ (0.1, 3) for values of r among {0.1, 0.3, 0.5, 0.7, 0.9}.

Figure 2 :

 2 Figure2: Values of Q(X) as a function of ψ 2 when X has the asymmetric dependence structure ℓ ψ 1 ,ψ 2 ,r and Fréchet(α) margins. Left plot (resp. right plot) gives the values of Q(X) for α = 1/2 (resp. α = 2) as a function of ψ 2 for different choices of ψ 1 , from top to bottom (resp. from bottom to top) ψ 1 = 0.1, . . . , 0.9, and the values r = 1/4, 3/4.

Lemma 14 .

 14 Let 0 < α < ∞. If X and Y are two random variables with regularly varying upper tails with index -α, then for 0 ≤ a ≤ ∞,

  in terms ofT i = {z ∈ Rd | z i > 1}, µ denoting the limit measure of Y • X. Define the random variable Z via the relationship VaR p Z = d i=1 VaR p (X i Y i ). Lemma 14 says that (32) is equivalent to lim x→∞ P( d i=1 X i Y i > x) P(Z > x) = Q(Y • X)D , which is also equivalent, by definition of Q(Y • X), to lim x→∞ d i=1 P (X i Y i > x) P(Z > x) = D .

  any i ∈ {1, . . . , d},lim p→0 VaR p (Z) VaR p (X i Y i ) = lim p→0 d j=1 VaR p (X j Y j ) VaR p (X i Y i ) =

Acknowledgements

This research is supported in part by the Agence Nationale de la Recherche through the AST&RISK project (ANR-08-BLAN-0314-01). The authors would like to express a special thank to Philippe Soulier for many fruitful discussions and advices during this work. They are also grateful to Stéphane Loisel who gave several constructive comments towards the final version of the paper.


Jean-Daniel Boissonnat
email: daniel.boissonnat@sophia.inria.fr

Arijit Ghosh
email: arijit.ghosh@sophia.inria.fr

Manifold Reconstruction using Tangential Delaunay Complexes *

Keywords: F.2.2 [Nonnumerical Algorithms and problems]: Geometrical problems and computations; I.3.5 [Computational Geometry and Object Modeling]: Curve, surface, solid, and object representations Algorithms, Theory Tangential Delaunay complex, Weighted Delaunay triangulation, manifold reconstruction, manifold learning, sampling conditions, sliver exudation

. Differently from previous methods, we do not construct any subdivision of the embedding d-dimensional space. As a result, the running time of our algorithm depends only linearly on the extrinsic dimension d while it depends quadratically on the size of the input sample, and exponentially on the intrinsic dimension k. To the best of our knowledge, this is the first certified algorithm for manifold reconstruction whose complexity depends linearly on the ambient dimension. We also prove that for a dense enough sample the output of our algorithm is isotopic to the manifold and a close geometric approximation of the manifold.

INTRODUCTION

Manifold reconstruction consists in computing a PL approximation of an unknown manifold M ⊂ R d from a finite sample of unorganized points P lying on M or close to M. When the manifold is a two-dimensional surface embedded in R 3 , the problem is known as the surface reconstruction problem. Surface reconstruction is a problem of major practical interest which has been extensively studied in the fields of Computational Geometry, Computer Graphics and Computer Vision. In the last decade, solid foundations have been established and the problem is now pretty well understood. Refer to Dey's book [START_REF] Dey | Curve and Surface Reconstruction: Algorithms with Mathematical Analysis[END_REF], and the survey by Cazals and Giesen in [START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF] for recent results. The output of those methods is a triangulated surface that approximates M. This triangulated surface is usually extracted from a 3-dimensional subdivision of the ambient space (typically a grid or a triangulation). Although rather inoffensive in 3-dimensional space, such data structures depend exponentially on the dimension of the ambient space, and all attempts to extend those geometric approaches to more general manifolds has led to algorithms whose complexities depend exponentially on d [START_REF] Niyogi | Finding the Homology of Submanifolds with High Confidence from Random Samples[END_REF][START_REF] Chazal | Smooth Manifold Reconstruction from Noisy and Non uniform Approximation with Guarantees[END_REF][START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF].

The problem in higher dimensions is also of great practical interest in data analysis and machine learning. In those fields, the general assumption is that, even if the data are represented as points in a very high dimensional space R d , they in fact live on a manifold of much smaller intrinsic dimension [START_REF] Seung | The manifold ways of perception[END_REF]. If the manifold is linear, well-known global techniques like principal component analysis (PCA) or multi-dimensional scaling (MDS) can be efficiently applied. When the manifold is highly nonlinear, several more local techniques have attracted much attention in visual perception and many other areas of science. Among the prominent algorithms are Isomap [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], Laplacian eigenmaps [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF], Hessian eigenmaps [START_REF] Donohu | Hessian eigenmaps: new locally linear embedding techniques for high dimensional data[END_REF], diffusion maps [START_REF] Lafon | Diffusion Maps and Coarse-Graining: A Unified Framework for Dimensionality Reduction, Graph Partitioning, and Data Set Parameterization[END_REF][START_REF] Nadler | Diffusion maps, spectral clustering and eigenfunctions of fokker-planck operators[END_REF], principal manifolds [START_REF] Zhang | Principal manifolds and nonlinear dimension reduction via local tangent space alignment[END_REF]. Most of those methods reduces to computing an eigendecomposition of some connection matrix. In all cases, the output is a mapping of the original data points into R k where k is the estimated intrinsic dimension of M. Those methods come with no or very limited guarantees. For example, Isomap provides a correct embedding only if M is isometric to a convex open set of R k . To be able to better approximate the sampled manifold, another route is to extend the work on surface reconstruction and to construct a PL approximation of M from the sample in such a way that, under appropriate sampling conditions, the quality of the approximation can be guaranteed. First investigations along this line can be found in the work of Cheng, Dey and Ramos [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF], and Boissonnat, Guibas and Oudot [START_REF] Boissonnat | Manifold reconstruction in arbitary dimensions using witness complexes[END_REF]. In both cases, however, the complexity of the algorithms is exponential in the ambient dimension d, which highly reduces their practical relevance.

In this paper, we extend the geometric techniques developped in small dimensions and propose a way to avoid computing data structures in the ambient space. We assume that M is a smooth manifold of known dimension k and that we can compute the tangent space to M at any sample point. Under those conditions, we propose a provably correct algorithm that allows to construct a simplicial complex of dimension k that approximates M. The complexity of the algorithm is linear in d, quadratic in the size n of the sample, and exponential in k. Our work builds on [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF] and [START_REF] Boissonnat | Manifold reconstruction in arbitary dimensions using witness complexes[END_REF] but dramatically reduces the dependance on d. To the best of our knowledge, this is the first certified algorithm for manifold reconstruction whose complexity depends only linearly on the ambient dimension. In the same spirit, Chazal and Oudot [START_REF] Chazal | Towards Persistence-Based Reconstruction in Euclidean Spaces[END_REF] have devised an algorithm of intrinsic complexity to solve the easier problem of computing the homology of a manifold from a sample.

Our approach is based on two main ideas : the notion of tangential Delaunay complex defined in [START_REF] Freedman | Efficient simplicial reconstructions of manifolds from their samples[END_REF][START_REF] Boissonnat | A coordinate system associated with points scattered on a surface[END_REF][START_REF] Flötotto | A coordinate system associated to a point cloud issued from a manifold: definition, properties and applications[END_REF], and the technique of sliver removal by weighting the sample points [START_REF] Cheng | Sliver Exudation[END_REF]. The tangential complex is obtained by gluing local (Delaunay) triangulations around each sample point. The tangential complex is a subcomplex of the d-dimensional Delaunay triangulation of the sample points but it can be computed using mostly operations in the k-dimensional tangent spaces at the sample points. Hence the dependence on k rather than d in the complexity. However, due to the presence of so-called inconsistencies, the local triangulations may not form a triangulated manifold. Although this problem has already been reported [START_REF] Freedman | Efficient simplicial reconstructions of manifolds from their samples[END_REF], no solution was known except for the case of curves (k = 1) [START_REF] Flötotto | A coordinate system associated to a point cloud issued from a manifold: definition, properties and applications[END_REF]. We show that we can remove inconsistencies by weighting the sample points under appropriate sample conditions. We can then prove that the approximation returned by our algorithm is isotopic to M, and a close geometric approximation of M.

Our algorithm can be seen as a local version of the cocone algorithm of Cheng et al. [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF]. By local, we mean that we do not compute any d-dimensional data structure like a grid or a triangulation of the ambient space. Still, the tangential complex is a subcomplex of the d-dimensional Delaunay triangulation of the data points and therefore implicitly relies on a global partition of the ambient space. This is key to our analysis and makes our method depart from other local algorithms that have been proposed in the surface reconstruction literature [START_REF] Cohen-Steiner | A greedy Delaunay Based Surface Reconstruction Algorithm[END_REF][START_REF] Gopi | Surface Reconstruction based on Lower Dimensional Localized Delaunay Triangulation[END_REF].

Notations In the rest of the paper, we assume that M is a smooth manifold of dimension k embedded in R d . We call P = {p1, . . . , pn} a finite sample of points from M. We denote by Tp the k-dimensional tangent space at point p ∈ M. We write B(c, r) for the d-dimensional ball centered at c of radius r. We define the angle between two vector spaces U and V as

U V = max u∈U min v∈V ∠uv. ( 1
)
If τ is a j-simplex, the (dj)-dimensional normal space of aff(τ ) is denoted by Nτ .

DEFINITIONS AND PRELIMINARIES

Weighted Delaunay triangulation

Weighted points

A weighted point is a pair consisting of a point p of R d , called the center of the weighted point, and a non-negative real number ω(p), called the weight of the weighted point. It might be convenient to visualize the weighted point (p, ω(p)) as the hypersphere (we will simply say sphere in the sequel) centered at p of radius ω(p).

Two weighted points (or spheres) (p, ω(p)) and (q, ω(q)) are called orthogonal when pq 2 = ω(p) 2 + ω(q) 2 , further than orthogonal when pq 2 > ω(p) 2 + ω(q) 2 , and closer than orthogonal when pq 2 < ω(p) 2 + ω(q) 2 .

Given a point set P = {p1, . . . , pn} ⊆ R d , a weight function on P is a non-negative real-valued function ω : P → [0, ∞). Write p ω i = (pi, ω(pi)) and P ω = {p ω 1 , . . . , p ω n }. We define the relative amplitude of ω, denoted as ω, as max p∈P,q∈P\{p} ω(p) ||p-q|| . In the paper, we assume that ω ≤ ω0 < 1/2, for some constant ω0 to be fixed later.

Given a subset τ of d + 1 weighted points whose centers are affinely independent, there exists a unique sphere orthogonal to the weighted points of τ . The sphere is called the orthosphere of τ and its center and radius are called the orthocenter and the orthoradius of τ . If τ is a j-simplex, j < d, the orthosphere of τ is the smallest sphere that is orthogonal to the (weighted) vertices of τ . Plainly, its center oτ lies in aff(τ ). The radius of the orthosphere of τ is denoted by R ′ τ . A finite set of weighted points P ω is said to be in general position if there exists no sphere orthogonal to d+2 weighted points of P ω .

Weighted Voronoi diagram and Delaunay triangulation

Let ω be a weight function defined over P. We define the weighted Voronoi cell of p ∈ P as

Vor ω (p) = {x ∈ R d : ||p -x|| 2 -ω 2 (p) ≤ ||q -x|| 2 -ω 2 (q), ∀q ∈ P}.
The weighted Voronoi cells and their k-dimensional faces, 0 ≤ k ≤ d, form a cell complex, called the weighted Voronoi diagram of P, that decomposes R d into convex polyhedral cells. See [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF]. Let τ be a subset of points of P and write Vor ω (τ ) = ∩x ∈ τ Vor ω (x). If the points of P are in general position, Vor ω (τ ) = ∅ when |τ | > d + 1. The collection of all simplices conv(τ ) such that Vor ω (τ ) = ∅ constitutes the weighted Delaunay triangulation Del ω (P). The mapping that associates to the face Vor ω (τ ) of Vor ω (P) the face conv(τ ) of Del ω (P) is a duality, i.e. a bijection that reverses the inclusion relation.

Alternatively, a d-simplex τ is in Del ω (P) if the orthosphere of τ is further than orthogonal from all weighted points in P ω \ {τ ω }.

The weighted Delaunay triangulation of a set of weighted points can be computed efficiently in small dimensions and has found many applications, see e.g. [START_REF] Aurenhammer | An optimal algorithm for constructing the weighted voronoi diagram in the plane[END_REF][START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF]. In this paper, we use weighted Delaunay triangulations for two main reasons. The first one is that the restriction of a d-dimensional weighted Voronoi diagram to an affine space of dimension k is a k-dimensional weighted Voronoi diagram that can be computed without computing the d-dimensional diagram (see Lemma 1). The other main reason is that some flat simplices named slivers can be removed from a Delaunay triangulation by weighting the vertices (see [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF][START_REF] Boissonnat | Manifold reconstruction in arbitary dimensions using witness complexes[END_REF][START_REF] Cheng | Sliver Exudation[END_REF] and Section 4).

p i p j x p ′ i p ′ j H Vor(p i p j )
Lemma 1. Let H be a k-dimensional affine space of R d .
The restriction of the weighted Voronoi diagram of P to H is the k-dimensional weighted Voronoi diagram of P ′ where P ′ is the orthogonal projection of P onto H and the squared weight of

p ′ i is ω 2 (pi) -pi -p ′ i 2 .
Proof. By Pythagoras theorem, we have

∀ x ∈ H ∩ Vor ω (pi), x -pi 2 -ω 2 (pi) ≤ x -pj 2 -ω 2 (pj) ⇔ x - p ′ i 2 + pi -p ′ i 2 -ω 2 (pi) ≤ x -p ′ j 2 + pj -p ′ j 2 -ω 2 (pj),
where p ′ i denotes the orthogonal projection of pi ∈ P onto H. Hence the restriction of Vor ω (P) to H is the weighted Voronoi diagram of the weighted points (p ′ i , ωi) ∈ H where

ω 2 i = -pi -p ′ i 2 + ω 2 (pi).

Lemma 2 ([13]

). If τ is a simplex of Del ω (P), then 1. ∀z ∈ aff(Vor ω (τ )) and ∀ p, q ∈ τ we have qz ≤

p-z √ 1-4ω 2 0 . 2. Rτ ≤ R ′ τ √ 1-4ω 2 0 . 3. ∀z ∈ aff(Vor ω (τ )) and ∀p ∈ τ , rz = p p -z 2 -ω 2 (p) ≥ R ′ τ .

Sampling conditions

Local feature size

The medial axis of M is the closure of the sets of points of R d that have more than one nearest neighbor on M. The local feature size of x ∈ M, lfs(x), is the distance of x to the medial axis of M. As is well known and can be easily proved, lfs is Lipschitz continuous i.e, lfs(x) ≤ lfs(y) + xy .

(ε, δ)-sample

The point sample P is said to be an (ε, δ)-sample (where 0 < δ < ε < 1) if (1) for any point x ∈ M there exists a point p ∈ P such that ||x -p|| ≤ ε lfs(x), and (2) for any two distinct points p, q ∈ P, ||p -q|| ≥ δ lfs(p). 1 The ratio ε/δ is called the relative density of P. 1 Observe that the sparsity condition ( 2) is mandatory if one wants to infer the dimension of M from a sample [START_REF] Giesen | Shape dimension and intrinsic metric from samples of manifolds[END_REF].

We will use the following results from [START_REF] Giesen | Shape dimension and intrinsic metric from samples of manifolds[END_REF]. We write lp for the distance between p ∈ P and its nearest neighbor in P \ {p}. Lemma 3. Given an (ε, δ)-sample P of M, we have 1. δ lfs(p) ≤ lp ≤ 2ε 1-ε lfs(p). 2. For any two points p, q ∈ M such that ||p-q|| = t lfs(p), 0 < t < 1, sin ∠(pq, Tp) ≤ t/2.

3.

Let p be a point in M. Let x be a point in Tp such that ||p -x|| ≤ t lfs(p) for some 0 < t ≤ 1/4. Let x ′ be the point on M closest to x. Then ||xx ′ || ≤ 2t2 lfs(p).

Slivers and good simplices

Consider a j-simplex τ , where 1 ≤ j ≤ k + 1. We denote by Rτ , Lτ , Vτ and ρ(τ ) = Rτ /Lτ the circumradius, the shortest edge length, the volume, and the radius-edge ratio of τ respectively. We further define σ(τ ) = `Vτ /L j τ ´1/j , as the sliverity measure of τ . The orthocenter of τ is denoted by oτ and its orthoradius by R ′ τ . If p ∈ τ , we define τp = τ \ {p} to be the (j -1)-face of τ opposite to p. We also write Dτ (p) for the distance from p to the affine hull of τp, and Hτ (p) for the signed distance 2 from oτ to aff(τp). We state without proof the following easy lemma. Lemma 4. Let τ be a simplex of Del ω (P) and p ∈ τ s.t.

1. There exists z ∈ aff(Vor ω (τ )) s.t. zp ≤ L1.

2. pq ≤ L2 for all vertices q of τ .

3. R ′ τ ≤ γ0Lτ . Then |Hτ (x)| = dist(oτ , aff(τx)) ≤ L1 + (1 + γ0 + ω0)L2 for all vertices x of τ .
Lemma 5 ([13]). Let τ be a simplex of Del ω (P). Let p be any vertex of τ and write Hτ (ω(p)) (instead of Hτ (p)) for the signed distance of oτ to aff(τp) parametrized by the weight of p. We have Hτ (ω(p)) = Hτ (0) -ω 2 (p) 2Dτ (p) . Slivers are a special type of flat simplices. The property of being a sliver is measured in terms of a parameter σ0, called the sliverity bound, to be fixed later in Section 4.

Definition 1 (Sliver). Given a positive parameter σ0, slivers are defined by induction on the dimension : (1) a simplex of dimension less than 3 is not a σ0-sliver, and (2) for j ≥ 3, a j-simplex τ is a σ0-sliver if σ(τ ) < σ0 and, ∀τ ′ ⊂ τ , σ(τ ′ ) ≥ σ0. Lemma 6. If a j-simplex τ is a σ0-sliver, then Dτ (p) < jσ0Lτ for all vertices p ∈ τ .

Proof. The volume of τ is

Vτ = Vτ p . Dτ (p) j = σ j-1 (τp)L j-1 τp . Dτ (p) j
and it is also equal to σ j (τ )L j τ . Since τ is a σ0-sliver, we have σ(τ ) < σ0 and σ(τp) ≥ σ0. Therefore we get Definition 2 (Good simplex). Given two positive constants ρ0 and σ0, a simplex τ is called a (ρ0, σ0)-good simplex if ρ(τ ) ≤ ρ0 and τ nor its subsimplices are σ0-slivers.

Dτ (p) = j σ j (τ ) σ j-1 (τp) × L j τ L j-1 τp < jσ0Lτ .
The following important lemma is known (see e.g. [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF]).

Lemma 7 (Normal approximation). Let τ be a (ρ0, σ0)good j-simplex for j ≤ k with vertices on a k-dimensional smooth manifold M, and p ∈ τ s.t. the lengths of the edges of τ that are incident to p are less than c3εlfs(p) for c3ε < 1. Then, for any normal vector np of M at p, τ has a normal nτ such that ∠npnτ ≤ αj (σ0)ε where αj (σ0) depends on j, σ0, ρ0 and c3, and 1/αj (σ0) vanishes with σ0.

In the sequel, we will simply write α(σ0) for α k (σ0).

Tangential Delaunay complex and inconsistent configurations

Let Del ω p i (P) be the weighted Delaunay triangulation of P restricted to the tangent space Tp i . Equivalently, the simplices of Del ω p i (P) are the simplices of Del ω (P) whose Voronoi dual faces intersect Tp i , i.e. τ ∈ Del ω p i (P) iff Vor ω (τ )∩ Tp i = ∅. Observe that Del ω p i (P) is in general a k-dimensional triangulation. Since this situation can always be ensured by applying some infinitesimal perturbation on P, we will assume, in the rest of the paper, that all Del ω p i (P) are kdimensional triangulations. Finally, write star(pi) for the star of pi in Del ω p i (P), i.e. the set of simplices that are incident to pi in Del ω p i (P). We call tangential Delaunay complex or tangential complex for short, the simplicial complex {τ, τ ∈ star(p), p ∈ P}. We denote it by Del ω T M (P). By our assumption above, Del ω T M (P) is a k-dimensional complex contained in Del ω (P). By duality, computing star(pi) is equivalent to computing the restriction of the (weighted) Voronoi cell of pi to Tp i , which, by Lemma 1, reduces to computing a cell in a kdimensional weighted Voronoi diagram embedded in Tp i . It follows that the tangential complex can be computed without constructing any data structure of dimension higher than k, the intrinsic dimension of M.

The tangential Delaunay complex is not in general a triangulated manifold and therefore not a good approximation of M. This is due to the presence of so-called inconsistencies. Consider a k-simplex τ of Del ω T M (P) with two vertices pi and pj such that τ is in star(pi) but not in star(pj) (refer to Fig. 3). We write Bi(τ ) for the open ball centered on Tp i that is orthogonal to the (weighted) vertices of τ ω , and denote by cp i and rp i its center and its radius. According to our definition, τ is inconsistent iff Bi(τ ) is further than orthogonal from all weighted points in P ω \ τ ω while there exists a weighted point of P ω \ τ ω , say p ω l , that is closer than orthogonal from Bj (τ ). We deduce from the above discussion that the line segment [cp i cp j ] has to penetrate the interior of Vor ω (p l ).

We formally define an inconsistent configuration as follows.

Definition 3 (Inconsistent configuration). φ = [p1 , . . . , p k+2 ] is called an inconsistent configuration of Del ω T M (P) witnessed by the triplet (pi, pj, p l ) if

• The k-simplex τ = φ \ {p l } is in star(pi) but not in star(pj).
• τ is a (ρ0, σ0)-good simplex.

• Vor ω (p l ) is the first cell of Vor ω (P) whose interior is intersected (at a point denoted by i φ ) by the line segment [cp i cp j ] oriented from cp i to cp j . Here cp i = Tp i ∩ Vor ω (τ ) and cp j = Tp j ∩ aff(Vor ω (τ )).

Note that i φ is the center of a sphere that is orthogonal to the weighted vertices of τ and also to p ω l , and further than orthogonal from all the other weighted points of P ω . Equivalenty, i φ is the point on [cp i cp j ] that belongs to Vor ω (φ).

An inconsistent configuration is thus a (k + 1)-simplex of Del ω (P). Such a configuration or its subfaces do not necessarily belong to Del ω T M (P) 3 We write In ω (P) for the subcomplex of Del(P) consisting of all the inconsistent configurations of Del ω T M (P) and their subfaces.

STRUCTURAL RESULTS

In the rest of the paper, we will need several constants : ω0, ρ0 and σ0 (to define slivers, good simplices and inconsistent configurations), and c3 (to be able to use lemma 7). 3 In fact, as already noted, no (k + 1)-simplex belongs to Del ω T M (P) when the points are in general position. These constants will be fixed in Section 3 and in Section 4. For simplicity, we will write sliver for σ0-sliver and inconsistent configuration for (ρ0, σ0)-inconsistent configuration.

p a q t x ′ x L/2 L/2 L/2 < L/2 > π/2
We give now an hypothesis which is assumed to be satisfied in the following results.

Hypothesis 1. P is an (ε, δ)-sample of M of bounded relative density, i.e. ε/δ ≤ η0 for some positive constant η0. We assume further that ε is small enough and that 2α(σ0)ε < 1. Finally, we assume that ω ≤ ω0 where ω is the relative amplitude of the weight assignement ω and ω0 is a positive constant less than 1/2.

Properties of the tangential Delaunay complex

We now give the following two lemmas which are slight variants of results of [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF].

Lemma 8. For all x ∈ Tp∩Vor ω (p), there exists a positive constant c1 such that ||p -x|| ≤ c1εlfs(p).

Proof. Assume for a contradiction that there exists a point x ∈ Vor ω (p) ∩ Tp s.t. ||p -x|| > c1εlfs(p) with c1(1 -c1ε) > 2 + c1ε(1 + c1ε) (*). Let q be a point on the line segment [px] such that ||p -q|| = c1εlfs(p)/2. Let q ′ be the point nearest to q on M. From Lemma 3, we have ||q -q

′ || ≤ c 2 1 ε 2 lfs(p)/2. Hence, ||p -q ′ || ≤ ||p -q|| + ||q -q ′ || < c 1 2 ε(1 + c1ε)lfs(p). From the 1-Lipschitz property, lfs(q ′ ) ≤ lfs(p) + ||p -q ′ || < (1 + c 1 2 ε(1 + c1ε)) lfs(p) < c 1 2 (1 -c1ε) lfs(p), which yields ||p -q ′ || ≥ ||p -q|| -||q -q ′ || > c1 2 ε(1 -c1ε) lfs(p)
Since P is an ε-sample, there exists a point t ∈ P, s.t.

||q ′ -t|| ≤ εlfs(q ′ ) < c 1 2 ε(1 -c1ε) lfs(p). We thus have ||q -t|| ≤ ||q -q ′ || + ||q ′ -t|| < c1 2 εlfs(p).
From Fig. 4, we can see that ∠ptx > π/2. This implies that ||x -

p|| 2 -||x -t|| 2 -||p -t|| 2 > 0. Hence, ||x -p|| 2 -||x -t|| 2 -ω 2 (p) + ω 2 (t) ≥ ||p -t|| 2 -ω 2 (p) ≥ ||p -t|| 2 -ω 2 0 .||p -t|| 2 > 0 (as ω0 < 1 2 )
This implies x ∈ Vor ω (p), which contradicts our initial assumption. We conclude that Vor

p-x √ 1-4ω 2 0 ≤ c 1 εlfs(p) √ 1-4ω 2 0 . Hence, p -q ≤ c ′ 2 εlfs(p) where c ′ 2 = c1(1 + 1/ p 1 -4ω 2 0 ) > 2c1. 1b.
From the definition of Del ω T M (P), there exists a vertex r of τ such that [pq] ∈ star(r). From 1a, rp and rq are at most c ′ 2 εlfs(r). From the 1-Lipschitz property of lfs and assuming that 2c ′ 2 ε < 1, we have lfs(r) ≤ lfs(p) + pr ≤ lfs(p)

1-c ′ 2 ε ≤ 2 lfs(p). It follows that p -q ≤ p -r + r -q ≤ 4c ′ 2 εlfs(p).
= p ||z -p|| 2 -ω 2 (p) ≤ ||z -p|| ≤ c1εlfs(p).
For any vertex q of τ , we have pq ≤ c2εlfs(p) (By 1.). Assuming 2c2ε ≤ 1 and using the fact that lfs is 1-Lipschitz, lfs(p) ≤ 2lfs(q). Therefore, taking for q a vertex of the shortest edge of τ , we have, using Lemma 3,

rz ≤ c1 " ε δ " δ lfs(p) ≤ c1 " ε δ " δ × 2lfs(q) ≤ 2c1η0 Lτ = ρ ′′ 0 Lτ .
From Lemma 2, we have Rτ ≤

R ′ τ √ 1-4ω 2 0 . Therefore ρ(τ ) ≤ ρ ′ 0 = ρ ′′ 0 √ 1-4ω 2 0 .

Properties of inconsistent configurations

We now give the lemmas on inconsistent configurations which are central to the proof of correctness of the reconstruction algorithm given later in the paper.

Lemma 10. Let φ ∈ In ω (P) be an inconsistent configuration witnessed by (pi, pj, p l ). Then, there exists positive constants c3 > c2 and ρ0 > ρ ′ 0 s.t.

1. pi φ ≤ c 3 2 εlfs(p) for all vertices p of φ.

2. If pq is an edge of φ then pq ≤ c3ε lfs(p).

3. If τ is a j-dimensional face of φ (j ≤ k + 1) and R ′ τ is the orthosphere of τ , then R ′ τ ≤ ρ0Lτ and Rτ /Lτ ≤ ρ0.

Proof. From the definition, τ = φ \ {p l } belongs to Del ω p i (P). We first bound dist(i φ , aff(τ )) = oτi φ where oτ is the orthocenter of τ . By Lemma 8, picp i ≤ c1εlfs(pi) and, by Lemma 2, we have

pj -cp i ≤ pi -cp i p 1 -4ω 2 0 ≤ c1ε lfs(pi) p 1 -4ω 2 0 Let c3 = c2/ p 1 -4ω 2 0 > c2. By Lemma 7, ∠(aff(τ ), Tp i ) ≤ α(σ0) ε, which implies that sin ∠(aff(τ ), Tp i ) ≤ α(σ0)ε and tan 2 ∠(aff(τ ), Tp i ) ≤ α 2 (σ0) ε 2 1 -α 2 (σ0) ε 2 < 4 α 2 (σ0) ε 2 ,
since 2α(σ0)ε < 1 (Hypothesis 1). Observing that ||pioτ || ≤ ||pi -cp i ||, since oτ belongs to aff(Vor ω (τ )) and therefore is the closest point to pi in aff(Vor ω (τ )), we deduce

cp i -oτ ≤ cp i -pi sin ∠(aff(τ ), Tp i ) ≤ α(σ0)c1ε 2 lfs(pi).
Moreover, ||pjoτ || ≤ ||pjcp i || as oτ is the closest point to pj in aff(Vor ω (τ )). Hence we have,

cp j -oτ ≤ pj -oτ tan ∠(aff(τ ), Tp j ) < 2 α(σ0) c1 ε 2 p 1 -4ω 2 0 lfs(pi) As i φ ∈ [cp i cp j ], we conclude that oτ -i φ ≤ 2 α(σ0) c1 ε 2 p 1 -4ω 2 0 lfs(pi).
1. Assuming that 2α(σ0) ε < 1 and using the facts that pioτ ≤ pii φ , we get

pi -i φ ≤ pi -oτ + oτ -i φ ≤ pi -cp i + oτ -i φ ≤ c1ε + 2α(σ0)c1ε 2 p 1 -4ω 2 0 ! lfs(pi) ≤ c2 4 εlfs(pi)
From Lemma 2, we have pi φ ≤

p i -i φ √ 1-4ω 2 0 = c 3 4
εlfs(pi) for all vertices p of φ.

2. Using part 1 of this lemma, we have

p -q ≤ p -i φ + q -i φ ≤ c3 2 εlfs(pi)
Moreover, lfs(pi) ≤ 2lfs(p) since lfs is a 1-Lipschitz function and by taking c3ε < 1. Hence pi φ ≤ c 3 2 εlfs(p) and pq ≤ c3ε lfs(p).

3. Let r φ = p i φpi 2ω 2 (pi). Since i φ ∈ Vor ω (τ ), the ball centered at i φ with radius r φ is orthogonal to the weighted vertices of τ . From Lemma 2, we have

r φ ≥ R ′ τ . Hence it suffices to show that r φ ≤ ρ0Lτ . Using pi -i φ ≤ c 3 4 εlfs(pi) we get r φ = p i φ -pi -ω 2 (pi) ≤ i φ -pi ≤ c3 4 εlfs(pi).
Let q be a vertex of a shortest edge of τ . We have, from part 2 of this lemma, piq ≤ c3 εlfs(q). Therefore

R ′ τ ≤ r φ ≤ c3 εlfs(q) ≤ c3 " ε δ " δ × lfs(q) ≤ c3 η0 Lτ
From Lemma 2, we have Rτ ≤ R ′ τ / p 1 -4ω 2 0 . Therefore , Rτ /Lτ ≤ ρ0 = c3η0/ p 1 -4ω 2 0 .

Lemma 11. Let τ be a j-simplex of Del ω T M (P) ∪ In ω (P), j ≤ k + 1. For any vertex p of τ , there exists a constant c4 s.t. the distance between the orthocenter o(τ ) of τ and aff(τp) is at most c4εlfs(p).

Proof. 1. We first consider the case where τ ∈ In ω (P). Then there exists an inconsistent configuration φ witnessed by points (pi, pj, p l ) s.t. τ is a j-dimensional subsimplex of φ. From Lemma 10, we have i φp ≤ c 3 2 εlfs(p), qp ≤ c3lfs(p) for all q ∈ φ and, since τ ∈ In ω (P), ρ(τ ) ≤ ρ0. Using the above facts and Lemma 4, we get dist(oτ , aff(τp)) ≤ dist(i φ , aff(τp)) 2. Consider now the case where τ ∈ Del ω T M (P). By definition, τ ∈ Del ω q (P) for some vertex q of τ . From Lemma 8, we have qcq ≤ c1εlfs(q) ≤ 2c1εlfs(p) where cq = Vor ω (τ )∩Tq and p is any vertex of τ . The last inequality follows from the facts that lfs is 1-Lipschitz, pq ≤ c2εlfs(p) (Lemma 9) and by taking c2ε ≤ 1. Therefore, using Lemma 4, we get dist(oτ , aff(τp)) ≤ 2c1εlfs(p) p 1 -4ω 2 0 + (1 + ρ ′ 0 + ω0)c2 εlfs(p).

The next lemma shows that inconsistent configurations are slivers provided that σ0 is sufficiently large wrt ε.

Lemma 12. Let φ be an inconsistent configuration and ε < f(σ0) = (k + 1) σ0 ρ0 (c3 + 2 α(σ0))

.

If the subfaces of φ are not slivers, then φ is a sliver.

Proof. Let pq be the smallest edge of φ and let r be a vertex in φ \ {p, q}. Since φr = φ \ {r} is not a sliver (as all subfaces of φ are not slivers) and rx ≤ c3lfs(x) (Lemma 10), we have sin ∠(aff(φr), Tx) ≤ α(σ0) ε for all vertices x of φr by Lemma 7. From Lemma 3, sin ∠(pr, Tp) ≤ c

Figure 1 :

 1 Figure 1: Refer to Lemma 1. The red line denotes the k-dimensional plane H and the black line denotes Vor ω (pipj ).

Figure 2 :

 2 Figure 2: [pq] and [pr] are edges of the star of p in Del ω T M (P) since their dual Voronoi edges intersect the tangent space Tp at p.

∈Figure 3 :

 3 Figure 3: An inconsistent configuration in the unweighted case. Edge [pipj] is in Delp i (P) but not in Delp j (P) since Vor(pipj) intersects Tp i but not Tp j . This happens because [cp i cp j ] penetrates (at i φ ) the Voronoi cell of a point p l = pi, pj, therefore creating an inconsistent configuration φ = [pi, pj , p l ].

Figure 4 :

 4 Figure 4: Refer to Lemma 8. x ′ is a point on the line segment such that ||px ′ || = c1ε lfs(p), L = c1ε lfs(p), ∠pax ′ = π/2 and ∠ptx ≥ ∠ptx ′ > π/2.

  ω0)c3εlfs(p) = c4εlfs(p),

  3 2 ε. ThereforeD φ (r) = sin ∠(pr, aff(φr)) × pr ≤ (sin ∠(pr, Tp) + sin ∠(aff(φr), Tp)) × pr ≤ ( c3 2 + α(σ0)) ε pr .Using the facts that vol(φr) = σ k (φr) L k φr , σ(φr) ≥ σ0, pr ≤ 2R φ , ρ(φr) ≤ ρ0 and ε < f(σ0), we get vol(φ) = D φ (r) × vol(φr)

  Lemma 9. There exists positive constants c2 and ρ ′ 0 s.t. 1. If pq is an edge of Del ω T M (P), then ||p -q|| ≤ c2εlfs(p). Consider first the case where pq is an edge of Del ω p (P). Then Tp ∩ Vor ω (pq) = ∅. Let x ∈ Tp ∩ Vor ω (pq). From Lemma 8, we have ||p -x|| ≤ c1εlfs(p). By Lemma 2, ||q -x|| ≤

	ω (p)∩Tp ⊆ B(p, c1εlfs(p)) if Inequality (*) is satisfied, which is true for all c1 = 3 + √ 2 ≈ 4.41 and ε < 0.09.
	2. If τ is a simplex of Del ω T M (P), then R ′ τ ≤ ρ ′ 0 Lτ and ρ(τ ) = Rτ /Lτ ≤ ρ ′ 0 .
	Proof. 1a.

  The first part of the lemma is proved by taking c2 = 4c ′ 2 > 8c1. 2. Assume that τ ∈ star(p). Let z ∈ Vor ω (τ ) ∩ Tp, and rz = p ||z -p|| 2ω 2 (p). By definition, the ball centered at z with radius rz is orthogonal to the weighted vertices of τ . From Lemma 2, we have rz ≥ R ′ τ . Hence it suffices to prove rz ≤ ρ ′ 0 Lτ . Since z ∈ Vor ω (τ ) ∩ Tp, we deduce from Lemma 8 that ||z -p|| ≤ c1εlfs(p). Therefore

	rz

H(p) is positive if oτ ∈ aff(τ ) and p lie on the same side of aff(τp), and it is negative if they lie on different sides of aff(τp).

ACKNOWLEDGEMENTS

The authors thank the anonymous referees for their insightful comments, and Mariette Yvinec for helpful discussions.

This work was partially supported by the ANR project GAIA. For the full version of the paper refer to [7].

Number of local neighbors

Lemmas 9 and 10 show that, in order to construct star(p) and search for inconsistencies involving p, it is enough to consider the points of P that lie in ball Bp = B(p, c3εlfs(p)). Since ε and lfs(p) are not known in practice, we will consider instead the ball B ′ p = B(p, c3η0 lp) where lx = min q∈P,q =x

xq .

It is easily seen that lx : M → R is 1-Lipschitz and, by Lemma 3, we have δlfs(p) ≤ lp ≤ 2ε 1-ε lfs(p). It follows that B ′ p contains Bp if ε/δ ≤ η0. We call LNp = B ′ p ∩ P the local neighborhood of p.

Lemma 13. The number of points of LNp is less than N < 2 O(k) .

Proof. For convenience, write ν = 4c3η0ε and assume that ν ≤ 1/2. We observe that LNp ⊂ B ′′ p = B(p, νlfs(p)) since, by Lemma 3, lp ≤ 2ε 1-ε lfs(p) ≤ 4εlfs(p). We will count the number of points in B ′′ p ∩ P. Let x and y be two points of B ′′ p ∩ P. We have lfs(x), lfs(y) ≥ lfs(p)(1ν) ≥ lfs(p)/2, since lfs is 1-Lipschitz. The balls B(x, lx/2) and B(y, ly/2) are disjoint, and, since lx ≥ δ lfs(x) ≥ δ 2 lfs(p) (and similarly for ly), the balls Bx = B(x, δ 4 lfs(p)) and By = B(y, δ 4 lfs(p)) are also disjoint. Observe that both balls Bx and By are contained in B + p = B(p, µ εlfs(p)) where µ = ν ε + δ 4ε ≤ 4c3η0 + 1 4 . A packing argument now allows to conclude. Specifically, by Lemma 15, we have that vol

, where φ k is the volume of the k-dimensional unit ball. We conclude that the number of points of P ∩ B ′′ p is less than `8µε

MANIFOLD RECONSTRUCTION

In this section, we will show how to find a weight assignment for the point set P so that we can remove all inconsistent configurations. Once this is done, all the stars become coherent and the resulting weighted tangential Delaunay complex is a simplicial k-manifold.

Algorithm

The input to the algorithm is a point sample P = {p1, . . . , pn} together with a bound on the relative density η0 of P. We assume in addition that the dimension k of M is given and that we know the tangent space Tp at each point p ∈ P.

The algorithm fixes ω0 in (0, 1/2) and then computes a weight assignment ω ∈ [0, ω0] such that no inconsistent configuration remains in the weighted tangential Delaunay complex. More precisely, we will assign a weight to each point p ∈ P in turn so as to remove all j-dimensional slivers incident to p in Del ω T M (P) ∪ In ω (P) for j ∈ {3, . . . , k + 1}. By Lemma 12, we know that, if σ0 is large enough, removing inconsistent configurations from the tangential Delaunay complex reduces to removing slivers.

For a given j-simplex τ = [p0, . . . , pj] ∈ Del ω T M (P) ∪ In ω (P) we have

the last inequalities follow from the facts that Lτ ≤ 2 Rτ , and the radius-edge ratio of the simplices in Del ω T M (P) ∪ In ω (P) is ≤ ρ0 (Lemmas 9 and 10). This implies that σ(τ ) ∈ (0, 2 ρ0]. In the first step of our reconstruction algorithm, we pick a random value of σ0 from (0, 2ρ0], and once σ0 is selected, we try to remove all slivers from Del ω T M (P) ∪ In ω (P). If we fail to remove slivers or if we still have inconsistencies, then we go back and select a new value of σ0.

Manifold-Reconstruction(P = {p1, . . . , pn}) S0 Calculate LNp i ∀pi ∈ P and ρ1 S1 Select σ0 at random from (0, 2 ρ1]

else output Del ω T M (P) Before we give the details of the function weight() we first define the subroutine skyline(p, S, σ0) that will be used in the function. Let S be a set of simplices incident on p.

Let τ be a simplex in S whose subfaces are not σ0-slivers. Such a simplex is called a candidate simplex. We associate to τ interval W (τ ) that consists of all squared weights ω 2 (p) for which τ appears as a simplex in Del ω T M (P) ∪ In ω (P). We define the skyline of p as the lower boundary of the union of all rectangles R(τ 1. We first detect all possible j-simplices for all 3 ≤ j ≤ k + 1 of Del ω T M (P) incident on p for all possible ω(p). This is done in the following way: (1) we vary the weight of p from 0 to to ω0lp, keeping the weights of the other points constant; (2) for each new weight assignment to p, we modify the stars of the points in LNp and detect from the stars the new j-simplices incident

We determine the next weight assignment of p in the following way. For each new simplex τ currently incident to p, we keep it in a priority queue ordered by the weight of p at which τ will disappear for the first time.

Hence the minimum weight in the priority queue gives the next weight assignment for p. Since the number of points in LNp is bounded, the number of simplices incident to p is also bounded, as well as the number of times we have to change the weight of p.

2. Once we have detected all possible j-simplices, for all j ∈ {1, . . . , k + 1}, that can be incident on p in the weighted tangential Delaunay complex, we then detect all possible inconsistent configurations incident on p, by calling the function detect inconsistentconfiguration(p, σ0).

Function detect inconsistent-configurations(p, σ0)

1. We vary the weight of p from 0 to ω0lp, keeping the weight of the rest of the points constant. Once we have assigned a new weight to p we modify the stars of the points in LNp.

2. Detecting the inconsistent configurations incident to p is more complicated than detecting the simplices incident to p. We consider all points pi in LNp. Let τ be a k-simplex in the star of pi, and let pj be a vertex of τ such that τ is not in the star of pj. We calculate the Voronoi diagram of the points in LNp restricted to the line segment [cp i cp j ], where cp i = Tp i ∩ Vor ω (τ ) and cp j = Tp j ∩ aff(Vor ω (τ )). From the restricted Voronoi diagram, we find a point p whose Voronoi cell intersects for the first time the line segment [cp i cp j ] oriented from cp i to cp j . If p ∈ φ = τ ∪ {l}, then we report φ.

3. As in the detect simplices function the weight of p is changed only a finite number of times. For each current inconsistent configuration φ incident to p, we keep in a priority queue the weight of p for which φ will disappear for the first time. The minimum weight in the priority queue gives the next weight assignment of p.

Analysis of the algorithm

Definition 4 (Sliverity range). Let ω be a weight assignment of relative amplitude at most ω0 we keep fixed except for ω(p). The sliverity range of a point p ∈ P is the measure of the set of all squared weights ω 2 (p) for which p is incident to a sliver in Del ω T M (P) ∪ In ω (P).

Lemma 14. Under Hypothesis 1, the sliverity range of p is at most Σ(p) = N k+2 (k + 1)c5 σ0 l 2 p for some constant c5. Proof. Let τ be a simplex incident on p. We call the sliverity range of τ the measure of the set of squared weights for which τ is a sliver in Del ω T M (P) ∪ In ω (P). If ω(p) is the weight of p, we write H(ω(p)) for the signed distance of the orthocenter of τ to aff(τp). From Lemma 11, we have |H(ω(p))| ≤ c4εlfs(p), for all τ ∈ Del ω T M (P) ∪ In ω (P). Moreover, using Lemma 5 and the fact that τ is a sliver,

2Dp ≤ H(0) -ω 2 (p) 2jσ 0 Lτ . It follows that the sliverity range of τ is at most 4jσ0Lτ c4εlfs(p). Using the facts that Lτ ≤ c3εlfs(p) (from Lemmas 9 and 10), lfs(p) ≤ lp/δ and ε/δ ≤ η0, the sliverity range of τ is less than 4jc3c4σ0η 2 0 l 2 p = j c5σ0 l 2 p . By Lemma 13, the number of j-simplices that are incident to p is at most N j . Hence, the sliverity range of p is less that

Theorem 1. Under Hypothesis 1 and if

then, for any σ0 ∈ [σmin, σmax], the above algorithm outputs Del ω T M (P) without any slivers nor inconsistent configuration. Proof. As in subroutine skyline(), let Sp denotes the set of all possible simplices that can be incident on p in the complex Del ω T M (P) ∪ In ω (P) for all possible weight assignments ω of relative amplitude ω ≤ ω0. By Lemma 14, the sliverity range of p is less than N k+2 (k + 1)c5 σ0 l 2 p . If the sliverity range of p is less than ω 2 0 l 2 p , the total range of all possible squared weights, or, equivalently, if

then we will be able to remove all slivers incident on p by selecting the highest point on the skyline. If we select a value of σ0 in the interval (σmin, σmax], Lemma 14 ensures that removing all j-dimensional slivers for all j ∈ {3, . . . , k + 1} in Del ω T M (P) ∪ In ω (P) will also result in removing inconsistent configurations from Del ω T M (P), i.e. In ω (P) = ∅.

Time and space complexity

Theorem 2. The time complexity of the algorithm is

where λ = (σmaxσmin)/(2ρ0), and its space complexity is

Proof. We only sketch the complexity analysis. See [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF] for a detailed discussion. Step S0 of Manifold-Reconstruction can easily be performed in O(d) |P| 2 time.

We show now that the expected number of times Step S1 of ManifoldReconstruction is repeated is 1 λ . Indeed, for any simplex τ ∈ Del ω T M (P)∪In ω (P), we have σ(τ ) ∈ (0, 2ρ0]. Hence, the probability that, for the selected value of σ0, the algorithm removes all slivers and inconsistencies is at least λ = (σmaxσmin)/(2ρ0). It follows that the expected number of times S1 is performed is less than

The time complexities of functions update(Del ω T M (P), p), detect simplices (p, σ0), and detect inconsistent-configurations(p, σ0) are 2 O(k(k+log d)) . Indeed, we need to project |LNp| < N < 2 O(k) points onto Tp, which costs O(kd) × |LNp| = O(d) 2 O(k) time. The Delaunay simplices that are computed have their vertices in LNp and are of dimension at most k + 1. Hence their number is 2 O(k 2 ) . The basic operations (mainly in-sphere predicates) amount to evaluating signs of determinants of k × k matrices. The cost of such a basic operation is O(k 3 ). The total cost of both functions is thus

We deduce that the expected time complexity of Manifold-Reconstruction is

We easily deduce from the above discussion that the total space complexity of the algorithm is

Topological and Geometric guarantees

We assume that the conditions of Theorem 1 are satisfied and denote by Del ω T M (P) the tangential complex output by our algorithm. Let π be the mapping that maps any point of Del ω T M (P) to its closest point on M. The proof of topological correctness (see [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF]) uses ideas from [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF][START_REF] Boissonnat | Manifold reconstruction in arbitary dimensions using witness complexes[END_REF][START_REF] Chazal | Smooth Manifold Reconstruction from Noisy and Non uniform Approximation with Guarantees[END_REF][START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF]. • Pointwise approximation : ∀x ∈ M, dist(x, π -1 (x))

= O(ε 2 lfs(x));

where τ is a k-simplex of Del ω T M (P) containing the point π -1 (x),

• Topological correctness : π defines an isotopy between Del ω T M (P) and M.

Proof sketch. By Lemma 7, one can show that the maximum distance from a point of a simplex τ ∈ Del ω T M (P) to the closest point on M is O(ε 2 lfs(x)). It follows that the projection π that maps every point of τ to its closest point on M is injective: if we extend an open segment of length lfs(y) from every manifold point y in all normal directions to M, these segments do not intersect, and they can be used as the fibers of a tubular neighborhood M of M. Each point of such a segment has y as its unique closest neighbor on M. For small enough ε, the simplex τ is contained in M. Thus, the mapping π defines an isotopy between τ and a corresponding manifold patch.

One can show that two k-simplices of Del ω T M (P) that share a subface have normal spaces that differ by at most O(ε), and that the mapping π extends continuously across the subfaces. It follows that the projection π restricted to the two adjacent k-simplices is a homeomorphism that is invertible locally. (In topological terms, π : Del ω T M (P) → M is a covering map, if we can establish that it is surjective.) By assumption, on every component, there is at least one vertex of a simplex of Del ω T M (P). This ensures that π(Del ω T M (P)) contains that vertex, and since the mapping can be continued locally, it follows that every component of S is covered at least once. It is now still possible that some component is covered more than once by π. This would imply that some sample point p ∈ P is covered more than once. However, one can show quite easily that no point p of Del ω T M (P) except p itself has p as its closest neighbor on M.

It follows that the mapping π defines an ambient isotopy between Del ω T M (P) and M.

CONCLUSION

We have given the first algorithm that is able to reconstruct a smooth manifold in a time that depends only linearly on the dimension of the ambient space. We believe that our algorithm is of practical interest when the dimension of the manifold is small, even if it is embedded in a space of very high dimension. This situation is quite common in practical applications in machine learning.

The algorithm is rather simple. The basic ingredients we need are data structures for constructing weighted Delaunay triangulations in k-flats. We will report experimental results in a forthcoming paper.

We have assumed that dimension of M is known. If not, we can use algorithms given in [START_REF] Giesen | Shape dimension and intrinsic metric from samples of manifolds[END_REF][START_REF] Cheng | Provable Dimension Detection using Principle Component Analysis[END_REF] to estimate the dimension of M and the tangent space at each sample point. One interesting feature of our approach is that it is pretty robust and still works if we only have approximate tangent spaces at the sample points.

We have also assumed that we know an upper bound on the relative density η0 of the input sample P. Ideas from [START_REF] Funke | Smooth-surface reconstruction in near-linear time[END_REF][START_REF] Boissonnat | Manifold reconstruction in arbitary dimensions using witness complexes[END_REF] may be useful to convert a sample to a subsample with a bounded relative density.

We forsee other applications of the tangential complex and of our construction each time computations in the tangent space of a manifold are required, e.g. for dimensionality reduction and approximating the Laplace Beltrami operator [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF].

Finally, let us mention that removing inconsistencies among stars that have been computed independently is a useful paradigm that has already been used for maintaining dynamic meshes [START_REF] Shewchuk | Star splaying: An algorithm for repairing delaunay triangulations and convex hulls[END_REF] and generating anisotropic meshes [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. We hope that this paper will motivate further applications. Lemma 15 ( [26]). Let A = B(p, r) ∩ M where r ≤ εlfs(M). Then, φ k r k ≥ vol(A)φ k r k /2 k , where φ k is volume of the k-dimensional unit ball.

APPENDIX