The empirical Christoffel function with applications in data analysis

TL;DR

This paper explores the empirical Christoffel function's applications in machine learning, providing methods for support estimation, consistency results, and demonstrating its effectiveness in density estimation, outlier detection, and affine matching.

Contribution

It introduces a thresholding scheme for support approximation, establishes a consistency result between empirical and population Christoffel functions, and demonstrates practical applications in data analysis.

Findings

Effective support estimation from finite samples.

Consistency between empirical and population Christoffel functions.

Successful applications in density estimation, outlier detection, and affine matching.

Abstract

We illustrate the potential applications in machine learning of the Christoffel function, or more precisely, its empirical counterpart associated with a counting measure uniformly supported on a finite set of points. Firstly, we provide a thresholding scheme which allows to approximate the support of a measure from a finite subset of its moments with strong asymptotic guaranties. Secondly, we provide a consistency result which relates the empirical Christoffel function and its population counterpart in the limit of large samples. Finally, we illustrate the relevance of our results on simulated and real world datasets for several applications in statistics and machine learning: (a) density and support estimation from finite samples, (b) outlier and novelty detection and (c) affine matching.

Tables

Table 1. Table 1: AUPR scores ( × 100 absent 100 \times 100 ) for the network intrusion detection task for the three different methods considered in this paper. 1SVM corresponds to one-class SVM with Gaussian kernel and varying kernel scale parameter σ 𝜎 \sigma and SVM parameter ν 𝜈 \nu . Christoffel corresponds to the empirical Christoffel function with varying degree d 𝑑 d . KDE corresponds to kernel density estimation with Gaussian kernel and varying scale parameter σ 𝜎 \sigma . The best scores are higlighted in boldface font.

1SVM
	0.005	0.01	0.02	0.05	0.1	0.2
0.01	10	17	17	17	15	11
0.02	2	17	18	17	15	12
0.05	8	1	14	18	15	11
0.1	9	8	12	17	14	11
0.2	7	9	8	17	14	13
0.5	3	5	9	15	17	16
1	3	6	9	14	19	18
2	4	4	5	1	18	18
5	4	3	4	9	12	16
10	5	4	4	7	10	15

Equations

$$\mathrm{vol}({\mathbf{B}}_{\delta }(\mathbf{x}))=\frac{{\pi }^{\frac{p}{2}}}{\Gamma \left(\frac{p}{2}+1\right)}{\delta }^p,\forall \mathbf{x}\in {\mathbb{R}}^p.$$

$$\mathbf{M}_{2}(\mu):\quad\begin{array}[]{rccccccccc}&&&1&X_{1}&X_{2}&X_{1}^{2}&X_{1}X_{2}&X_{2}^{2}\\ &&&\\ 1&&&1&y_{10}&y_{01}&y_{20}&y_{11}&y_{02}\\ X_{1}&&&y_{10}&y_{20}&y_{11}&y_{30}&y_{21}&y_{12}\\ X_{2}&&&y_{01}&y_{11}&y_{02}&y_{21}&y_{12}&y_{03}\\ X_{1}^{2}&&&y_{20}&y_{30}&y_{21}&y_{40}&y_{31}&y_{22}\\ X_{1}X_{2}&&&y_{11}&y_{21}&y_{12}&y_{31}&y_{22}&y_{13}\\ X_{2}^{2}&&&y_{02}&y_{12}&y_{03}&y_{22}&y_{13}&y_{04}\\ \end{array}.$$

$$\exists Q\in {\mathbb{R}}^{s(m)\times s(m)},Q=Q^T,Q⪰0,P(\mathbf{x})={\mathbf{v}}_m(\mathbf{x}{)}^{T}Q{\mathbf{v}}_m(\mathbf{x}),\forall \mathbf{x}\in {\mathbb{R}}^p.$$

$$⟨P_{\alpha },P_{\beta }{⟩}_{\mu }={\delta }_{\alpha =\beta },⟨P_{\alpha },X^{\beta }{⟩}_{\mu }=0,\mathrm{if}\beta {<}_{gl}\alpha ,⟨P_{\alpha },X^{\alpha }{⟩}_{\mu }>0.$$

$$(\mathbf{x},\mathbf{y})\mapsto {\kappa }_{\mu ,d}(\mathbf{x},\mathbf{y}):=\sum _{\alpha \in {\mathbb{N}}_d^p}{P}_{\alpha }(\mathbf{x})P_{\alpha }(\mathbf{y})={\mathbf{v}}_d(\mathbf{x}{)}^T{\mathbf{M}}_d(\mu {)}^{-1}{\mathbf{v}}_d(\mathbf{y}),$$

$$P_{\alpha }(X)=\int K(X,\mathbf{y})P_{\alpha }(\mathbf{y})d\mu (\mathbf{y}),\forall \alpha \in {\mathbb{N}}^p,$$

$$p(X)=\int K(X,\mathbf{y})p(\mathbf{y})d\mu (\mathbf{y}),\forall p\in \mathbb{R}[X].$$

$${\Lambda }_{\mu ,d}(\text{boldmath}\xi )=\underset{P\in \mathbb{R}[X{]}_d}{\min }\left\{{\int }_{{\mathbb{R}}^p}P(\mathbf{x}{)}^{2}d\mu (\mathbf{x}):P(\text{boldmath}\xi )=1\right\}.$$

$$P_d^{\ast }(X)=\frac{{\kappa }_{\mu ,d}(X,\text{boldmath}\xi )}{{\kappa }_{\mu ,d}(\text{boldmath}\xi ,\text{boldmath}\xi )}={\Lambda }_{\mu ,d}(\xi ){\kappa }_{\mu ,d}(X,\text{boldmath}\xi ).$$

$${\Lambda }_{\mu ,d}(\text{boldmath}\xi )$$

$${\Lambda }_{\mu ,d}(\text{boldmath}\xi ){\text{boldmath}\xi }^{\alpha }$$

$$\displaystyle\min_{P\in\mathbb{R}[\mathbf{x}]_{d}}\,\left\{\,\int P^{2}(\mathbf{x})\,d\mu(\mathbf{x}):\>P(\mathbf{x})=\mathbf{x}^{\alpha}+\sum_{\beta<_{gl}\,\alpha}\theta_{\beta}\,\mathbf{x}^{\beta}\,\quad\mbox{for some $\{\theta_{\beta}\}_{\beta<_{gl}\alpha}$}\right\}.$$

$${\mathbf{v}}_d(\mathbf{x}{)}^T{\mathbf{M}}_d(\mu {)}^{-1}{\mathbf{v}}_d(\mathbf{y})$$

$$=(A{\mathbf{v}}_d(\mathbf{x}){)}^T{\left({\int }_{{\mathbb{R}}^p}(A{\mathbf{v}}_d(\mathbf{z}))(A{\mathbf{v}}_d(\mathbf{z}){)}^{T}d\mu (\mathbf{z})\right)}^{-1}(A{\mathbf{v}}_d(\mathbf{y}))$$

$$2^{3-\frac{{\delta }_k{d}_k}{{\delta }_k+\mathrm{diam}(S)}}{d}_k^p{\left(\frac{e}{p}\right)}^p\exp \left(\frac{p^2}{d_k}\right)\le {\alpha }_k$$

$${\alpha }_k:=\frac{{\delta }_k^p{\omega }_p}{\lambda (S)}\frac{(d_k+1)(d_k+2)(d_k+3)}{(d_k+p+1)(d_k+p+2)(2d_k+p+6)}$$

$$d_H(X,Y)=\max \left\{\underset{\mathbf{x}\in X}{\sup }\underset{\mathbf{y}\in Y}{\inf }\mathrm{dist}(\mathbf{x},\mathbf{y}),\underset{\mathbf{y}\in Y}{\sup }\underset{\mathbf{x}\in X}{\inf }\mathrm{dist}(\mathbf{x},\mathbf{y})\right\}.$$

$$S_k:=\left\{\mathbf{x}\in {\mathbb{R}}^p:s(d_k){\Lambda }_{{\mu }_S,d_k}(\mathbf{x})\ge {\alpha }_k\right\}.$$

$$d_H(S_k,S)$$

$$d_H(\partial S_k,\partial S)$$

$$2^{3-\frac{{\delta }_k{d}_k}{{\delta }_k+\mathrm{diam}(S)}}{d}_k^p{\left(\frac{e}{p}\right)}^p\exp \left(\frac{p^2}{d_k}\right)\le {\alpha }_k$$

$${\alpha }_k:=w_{-}{\delta }_k^p{\omega }_p\frac{(d_k+1)(d_k+2)(d_k+3)}{(d_k+p+1)(d_k+p+2)(2d_k+p+6)}.$$

$$S_k:=\left\{\mathbf{x}\in {\mathbb{R}}^p:s(d_k){\Lambda }_{{\mu }_S,d_k}(\mathbf{x})\ge {\alpha }_k\right\}.$$

$$d_H(S_k,S)$$

$$d_H(\partial S_k,\partial S)$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^p}{\sup }\left\{|{\Lambda }_{{\mu }_n,d}(\mathbf{x})-{\Lambda }_{\mu ,d}(\mathbf{x})|\right\}\stackrel{a.s.}{\underset{n\to \infty }{⟶}}0.$$

$$\Vert {\Lambda }_{{\mu }_n,d}-{\Lambda }_{\mu ,d}{\Vert }_{\infty }\stackrel{a.s.}{\underset{n\to \infty }{⟶}}0.$$

$$d_H(\partial L_n,\partial L)\stackrel{a.s.}{⟶}0,$$

$$\underset{d\to \infty }{\lim }d{\Lambda }_{\mu ,d}(x)\to \pi h(x)\sqrt{1-x^2}.$$

$$\int P^{2}d\mu =P^T{\mathbf{M}}_d(\mu )P\text{mbox}and{\mathbf{M}}_d(\mu )\succ 0,$$

$$P_d^{\ast }(\text{boldmath}\xi )=1;2{\mathbf{M}}_d(\mu )P_d^{\ast }=\theta {\mathbf{v}}_d(\xi ),$$

$$2{\kappa }_{\mu ,d}(\text{boldmath}\xi ,\text{boldmath}\xi {)}^{-1}=2(P_d^{\ast }{)}^T{\mathbf{M}}_d(\mu )P_d^{\ast }=\theta .$$

Full text

1 Introduction

The main claim of this paper is that the Christoffel function (a tool from Approximation Theory) can prove to be very useful machine learning applications. The Christoffel function is associated with a finite measure and a degree parameter $d$. It has an important history of research with strong connection to orthogonal polynomials [Szegö (1974), Dunkl and Xu (2001)], interpolation and approximation theory [Nevai (1986), De Marchi et al. (2014)]. Its typical asymptotic behavior as $d$ increases is of particular interest because, in some specific settings, it provides very relevant information on the support and the density of the associated input measure. Important references include [Máté and Nevai (1980), Máté et al. (1991), Totik (2000), Gustafsson et al. (2009)] in a single dimension, [Bos (1994), Bos et al. (1998), Xu (1999), Berman (2009), Kroo and Lubinsky (2013)] for specific multivariate settings and [Kroo and Lubinsky (2012)] for ratii of mutually absolutely continuous measures. The topic is still a subject of active research, but regarding properties of the Christoffel function, a lot of information is already available.

The present work shows how properties of the Christoffel function can be used successfully in some machine learning applications. To the best of our our knowledge, this is the first attempt in such a context with the recent work of [Lasserre and Pauwels (2016)] and [Malyshkin (2015)]. More precisely, we consider the empirical Christoffel function, a specific case where the input measure is a scaled counting measure uniformly supported on a set (a cloud) of datapoints. This methodology has three distinguishing features: (i) It is extremely simple and involves no optimization procedure, (ii) it scales linearly with the number of observations (one pass over the data is sufficient), and (iii) it is affine invariant. These three features prove to be especially important in all the applications that we consider.

In [Lasserre and Pauwels (2016)] we have exhibited a striking property of some distinguished family of sum-of-squares (SOS) polynomials $(Q_{d})_{d\in\mathbb{N}}$, indexed by their degree ($2d\in\mathbb{N}$), and easily computed from empirical moments associated with a cloud of $n$ points in $\mathbb{R}^{p}$ which we call $\mathbf{X}$. The associated family of sublevel sets $S_{\alpha,d}=\{\mathbf{x}:Q_{d}(\mathbf{x})\leq\alpha\}$, for various values of $\alpha>0$, approximates the global shape of original cloud of points $\mathbf{X}$. The degree index $d$ can be used as a tuning parameter, trading off regularity of the polynomial $Q_{d}$ with the fitness of the approximation of the shape of $\mathbf{X}$ (as long as the cloud contains sufficiently many points). Remarkably, even with relatively low degree $d$, the sets $S_{\alpha,d}$ capture acurately the shape of $\mathbf{X}$, and so provides a compact (algebraic) encoding of the cloud.

In fact the reciprocal function $\mathbf{x}\mapsto Q_{d}(\mathbf{x})^{-1}$ is precisely the Christoffel function $\Lambda_{\mu_{n},d}$ associated to the empirical counting measure supported on $\mathbf{X}$ and the degree index $d$. Some properties of the Christoffel function stemming from approximation theory suggest that it could be exploited in a statistical learning context by considering its empirical counterpart. The purpose of this work is to push this idea further. In particular we investigate (a) further properties of the Christoffel function which prove to be relevant in some machine learning applications, (b) statistical properties of the empirical Christoffel function as well as (c) further applications to well known machine learning tasks.

Contributions

This paper significantly extends [Lasserre and Pauwels (2016)] in several directions. Indeed our contribution is threefold:

I. We first provide a thresholding scheme which allows to approximate the compact support $S$ of a measure with strong asymptotic guarantees. This result rigorously establishes the property that, as $d$ increases, the scaled Christoffel function decreases to zero outside $S$ and remains positive in the interior of $S$.

II. In view of potential applications in machine learning we provide a rationale for using the empirical Christoffel function in place of its population counterpart in the limit of large sample size. We consider a compactly supported population measure $\mu$ as well as an empirical measure $\mu_{n}$ uniformly supported on a sample of $n$ vectors in $\mathbb{R}^{p}$, drawn independently from $\mu$. For each fixed $d$ we show a highly desirable strong asymptotic property as $n$ increases. Namely, the empirical Christofell function $\Lambda_{\mu_{n},d}(\cdot)$ converges, uniformly in $\mathbf{x}\in\mathbb{R}^{p}$, to $\Lambda_{\mu,d}(\cdot)$, almost-surely with respect to the draw of the random sample.

III. We illustrate the benefits of the empirical Christoffel function in some important applications, mainly in machine learning. The rationale for such benefits builds on approximation properties of the Christoffel function combined with our consistency result. In particular, we first show on simulated data that the Christoffel function can be useful for density estimation and support inference. In [Lasserre and Pauwels (2016)] we have described how the Christoffel function yields a simple procedure for intrusion detection in networks, and here we extend these results by performing a numerical comparison with well established methods for novelty detection on a real world dataset. Finally we show that the Christoffel function is also very useful to perform affine matching and inverse affine shuffling of a dataset.

Comparison with existing literature on set estimation

Support estimation and more generaly set estimation has a long history in statistics and we intend to give a nonexhaustive overview in this section. The main question of interest is that of inferering a set (support, level sets of the density function …) based on independants samples from an unknown distribution. Pioneering works include [Rényi and Sulanke (1963), Geffroy (1964)] followed by [Chevalier (1976), Devroye and Wise (1980)] and resulted in the introduction and first analyses for estimators based on convex hull for convex domains or union of balls for nonconvex sets. This motivated the development of minimax statistical analysis for the set estimation problem [Hardle et al. (1995), Mammen and Tsybakov (1995), Tsybakov (1997)] and the introduction of more sophisticated optimal estimators, such as the excess mass estimator [Polonik (1995)]. Strong relations between set estimation and density estimation lead to the development of the plugin approach for support and density level set estimation [Cuevas and Fraiman (1997), Molchanov (1998)] with futher generalization proposed in [Cuevas et al. (2006)] and a precise minimax analysis described in [Rigollet and Vert (2009)].

These works provide a rich statitical analysis of the main estimation approaches currently available. The topic is still active with more precise questions ranging from inference of topological properties [Aaron and Bodart (2016)], new geometric conditions [Cholaquidis et al. (2014)], adaptivity to local properties of the underlying density [Patschkowski and Rohde (2016), Singh et al. (2009)].

One of the goals of our work is the introduction of the Christoffel function as a tool to solve similar problems. This approach has several advantages

• •

The Christoffel function allows to encode the global shape of a cloud of points in any finite dimension using a polynomial level set. This kind of encoding is relatively simple and compact. This has clear advantages, for example, the evaluation of a polynomial has a complexity which does not depend on the size of the sample used to qualibrate its coefficients and the boundary of the corresponding sublevel set as a very compact representation as an algebraic set. Furthermore, it turns out that the estimation of the empirical Christoffel function has a computational cost which is linear in the sample size. This is in contrast with distance based approaches for which membership evaluation requires to query all the sample points. As pointed out in [Baíllo et al. (2000)], the practical use of multidimensional set estimation techniques involves formidable computational difficulties so that simplicity arises as a major advantage in this context.

• •

The proposed approach is specific in the sense that it relies on tools which were not considered before for support estimation such as orthogonal polynomials. Topological properties of the support of the distribution or its boundary arise as major questions beyond minimax analysis [Aaron and Bodart (2016)]. In this realm, the objects which we manipulate have a simple algebraic description and could be coupled with computational real algebraic geometry tools to infer topological properties such as, for example, Betty numbers [Basu et al. (2005)]. This strong algebraic structure could in principle allow to push further the statistical settings which could be handled, with, for example, notions such as singular measures and intrinsic dimension.

We see these facts as potential advantages of the Christoffel function in the context of support estimation and relevant motivation to further study the potential of this procedure in modern data analysis contexts. However, we emphasize that this work constitutes only a first step in this direction. Indeed, we are not able to provide a complete statistical efficiency analysis as precisely described in the support and set estimation literature (e.g. [Cuevas et al. (2006)]). This would require further studies of precise properties of the Christoffel function itself which are not available given the state of knowledge for this object. We aim at providing a rationale for the proposed approach and motivation for future studies, among which a complete statistical analysis is a longer term goal.111In particular qualibration of the underlying polynomial degree as a function of the sample size is out of the scope of this paper and left for future research.

Organisation of the paper

Section 2 describes the notation and definitions which will be used throughout the paper. In Section 3 we introduce the Christoffel function, outline some of its known properties and describe our main theoretical results. Applications are presented in Section 4 where we consider both simulated and real world data as well as a comparison with well established machine learning methods. For clarity of exposition most proofs and technical details are postponed to the Appendix in Section 6.

2 Notation, definitions and Preliminary results

2.1 Notation and definitions

We fix the ambient dimension to be $p$ throughout the text. For example, we will manipulate vectors in $\mathbb{R}^{p}$ as well as $p$-variate polynomials with real coefficients. We denote by $X$ a set of $p$ variables $X_{1},\ldots,X_{p}$ which we will use in mathematical expressions defining polynomials. We identify monomials from the canonical basis of $p$-variate polynomials with their exponents in $\mathbb{N}^{p}$: we associate to $\alpha=(\alpha_{i})_{i=1\ldots p}\in\mathbb{N}^{p}$ the monomial $X^{\alpha}:=X_{1}^{\alpha_{1}}X_{2}^{\alpha_{2}}\ldots X_{p}^{\alpha_{p}}$ which degree is $\deg(\alpha):=\sum_{i=1}^{p}\alpha_{i}=|\alpha|$. We use the expressions $<_{gl}$ and $\leq_{gl}$ to denote the graded lexicographic order, a well ordering over $p$-variate monomials. This amounts to, first, use the canonical order on the degree and, second, break ties in monomials with the same degree using the lexicographic order with $X_{1}=a,X_{2}=b\ldots$ For example, the monomials in two variables $X_{1},X_{2}$, of degree less or equal to $3$ listed in this order are given by: $1,\,X_{1},\,X_{2},\,X_{1}^{2},\,X_{1}X_{2},\,X_{2}^{2},\,X_{1}^{3},\,X_{1}^{2}X_{2},\,X_{1}X_{2}^{2},\,X_{2}^{3}$.

We denote by $\mathbb{N}^{p}_{d}$, the set $\left\{\alpha\in\mathbb{N}^{p};\;\deg(\alpha)\leq d\right\}$ ordered by $\leq_{gl}$. $\mathbb{R}[X]$ denotes the set of $p$-variate polynomials: linear combinations of monomials with real coefficients. The degree of a polynomial is the highest of the degrees of its monomials with nonzero coefficients222For the null polynomial, we use the convention that its degree is [math] and it is $\leq_{gl}$ smaller than all other monomials.. We use the same notation, $\deg(\cdot)$, to denote the degree of a polynomial or of an element of $\mathbb{N}^{p}$. For $d\in\mathbb{N}$, $\mathbb{R}[X]_{d}$ denotes the set of $p$-variate polynomials of degree at most $d$. We set $s(d)={p+d\choose d}$, the number of monomials of degree less or equal to $d$.

We will denote by $\mathbf{v}_{d}(X)$ the vector of monomials of degree less or equal to $d$ sorted by $\leq_{gl}$, i.e., $\mathbf{v}_{d}(X):=\left(X^{\alpha}\right)_{\alpha\in\mathbb{N}^{p}_{d}}\in\mathbb{R}[X]^{s(d)}_{d}$. With this notation, we can write a polynomial $P\in\mathbb{R}[X]_{d}$ as $P(X)=\left\langle\mathbf{p},\mathbf{v}_{d}(X)\right\rangle$ for some real vector of coefficients $\mathbf{p}=\left(p_{\alpha}\right)_{\alpha\in\mathbb{N}_{d}^{p}}\in\mathbb{R}^{s(d)}$ ordered using $\leq_{gl}$. Given $\mathbf{x}=(x_{i})_{i=1\ldots p}\in\mathbb{R}^{p}$, $P(\mathbf{x})$ denotes the evaluation of $P$ with the assignments $X_{1}=x_{1},X_{2}=x_{2},\ldots X_{p}=x_{p}$. Given a Borel probability measure $\mu$ and $\alpha\in\mathbb{N}^{p}$, $y_{\alpha}(\mu)$ denotes the moment $\alpha$ of $\mu$, i.e., $y_{\alpha}(\mu)=\int_{\mathbb{R}^{p}}\mathbf{x}^{\alpha}d\mu(\mathbf{x})$. Finally for $\delta>0$ and every $\mathbf{x}\in\mathbb{R}^{p}$, let $\mathbf{B}_{\delta}(\mathbf{x}):=\{\mathbf{x}:\|\mathbf{x}\|\leq\delta\}$ be the closed Euclidean ball of radius $\delta$ and centered at $\mathbf{x}$. We use the shorthand notation $\mathbf{B}$ to denote the closed Euclidean unit ball. For a given subset of Euclidean space, $A$, $\partial A$ denotes the topological boundary of $A$. Recall that its Lebesgue volume ${\rm vol}(\mathbf{B}_{\delta}(\mathbf{x}))$ satisfies:

[TABLE]

Furthermore, let $\omega_{p}:=\frac{2\pi^{\frac{p+1}{2}}}{\Gamma\left(\frac{p+1}{2}\right)}$ denote the surface of the $p$ dimensional unit sphere in $\mathbb{R}^{p+1}$. Throughout the paper, we will only consider measures of which all moments are finite.

Moment matrix

For a finite Borel measure $\mu$ on $\mathbb{R}^{p}$ denote by ${\rm supp}(\mu)$ its support, i.e., the smallest closed set $\mathbf{\Omega}\subset\mathbb{R}^{p}$ such that $\mu(\mathbb{R}^{p}\setminus\mathbf{\Omega})=0$. The moment matrix of $\mu$, $\mathbf{M}_{d}(\mu)$, is a matrix indexed by monomials of degree at most $d$ ordered by $\leq_{gl}$. For $\alpha,\beta\in\mathbb{N}^{p}_{d}$, the corresponding entry in $\mathbf{M}_{d}(\mu)$ is defined by $\mathbf{M}_{d}(\mu)_{\alpha,\beta}:=y_{\alpha+\beta}(\mu)$, the moment $\int\mathbf{x}^{\alpha+\beta}d\mu$ of $\mu$. When $p=2$ and $d=2$, letting $y_{\alpha}=y_{\alpha}(\mu)$ for $\alpha\in\mathbb{N}_{4}^{2}$, we have

[TABLE]

The matrix $\mathbf{M}_{d}(\mu)$ is positive semidefinite for all $d\in\mathbb{N}$. Indeed, for any $\mathbf{p}\in\mathbb{R}^{s(d)}$, let $P\in\mathbb{R}[X]_{d}$ be the polynomial with vector of coefficients $\mathbf{p}$; then $\mathbf{p}^{T}\mathbf{M}_{d}(\mu)\mathbf{p}=\int_{\mathbb{R}^{p}}P(\mathbf{x})^{2}d\mu(\mathbf{x})\geq 0$. We also have the identity $\mathbf{M}_{d}(\mu)=\int_{\mathbb{R}^{p}}\mathbf{v}_{d}(\mathbf{x})\mathbf{v}_{d}(\mathbf{x})^{T}d\mu(\mathbf{x})$ where the integral is understood elementwise.

Sum of squares (SOS)

We denote by $\Sigma[X]\subset\mathbb{R}[X]$ (resp. $\Sigma[X]_{d}\subset\mathbb{R}[X]_{d}$), the set of polynomials (resp. polynomials of degree at most $d$) which can be written as a sum of squares of polynomials. Let $P\in\mathbb{R}[X]_{2m}$ for some $m\in\mathbb{N}$, then $P$ belongs to $\Sigma[X]_{2m}$ if there exists a finite $J\subset\mathbb{N}$ and a family of polynomials $P_{j}\in\mathbb{R}[X]_{m}$, $j\in J$, such that $P=\sum_{j\in J}P_{j}^{2}$. It is obvious that sum of squares polynomials are always nonnegative. A further interesting property is that this class of polynomials is connected with positive semidefiniteness. Indeed, $P$ belongs to $\Sigma[X]_{2m}$ if and only if

[TABLE]

As a consequence, every real symmetric positive semidefinite matrix $Q\in\mathbb{R}^{s(m)\times s(m)}$ defines a polynomial in $\Sigma[X]_{2m}$ by using the representation (2.1).

Orthonormal polynomials

We define a classical [Szegö (1974), Dunkl and Xu (2001)] family of orthonormal polynomials, $\left\{P_{\alpha}\right\}_{\alpha\in\mathbb{N}^{p}_{d}}$ ordered according to $\leq_{gl}$, which satisfies for all $\alpha\in\mathbb{N}^{p}_{d}$

[TABLE]

Existence and uniqueness of such a family is guaranteed by the Gram-Schmidt orthonormalization process following the $\leq_{gl}$ ordering on monomials and by the positivity of the moment matrix, see for instance [Dunkl and Xu (2001)] Theorem 3.1.11.

Let $\mathbf{D}_{d}(\mu)$ be the lower triangular matrix of which rows are the coefficients of the polynomials $P_{\alpha}$ defined in (2.2) ordered by $\leq_{gl}$. It can be shown that $\mathbf{D}_{d}(\mu)=\mathbf{L}_{d}(\mu)^{-T}$, where $\mathbf{L}_{d}(\mu)$ is the Cholesky factorization of $\mathbf{M}_{d}(\mu)$. Furthermore, there is a direct relation with the inverse moment matrix as $\mathbf{M}_{d}(\mu)^{-1}=\mathbf{D}_{d}(\mu)^{T}\mathbf{D}_{d}(\mu)$ ([Helton et al. (2008)] Proof of Theorem 3.1).

3 The Christoffel function and its empirical counterpart

3.1 The Christoffel function

Let $\mu$ be a finite Borel measure on $\mathbb{R}^{p}$ with all moments finite and such that its moment matrix $\mathbf{M}_{d}(\mu)$ is positive definite for every $d=0,1,\ldots$. For every $d$, define the function $\kappa_{\mu,d}:\mathbb{R}^{p}\times\mathbb{R}^{p}\to\mathbb{R}$ by:

[TABLE]

where the family of polynomials $\left\{P_{\alpha}\right\}_{\alpha\in\mathbb{N}_{d}^{p}}$ is defined in (2.2) and the last equality follows from properties of this family (see also [Lasserre and Pauwels (2016)]). The kernel $(\mathbf{x},\mathbf{y})\mapsto K(\mathbf{x},\mathbf{y}):=\,\sum_{\alpha\in\mathbb{N}^{p}}P_{\alpha}(\mathbf{x})\,P_{\alpha}(\mathbf{y})$ is a reproducing kernel on $L_{2}(\mu)$ because

[TABLE]

that is, the $(P_{\alpha})$ are eigenvectors of the associated operator on $L_{2}(\mu)$, and so

[TABLE]

The function $\mathbf{x}\mapsto\Lambda_{\mu,d}(\mathbf{x}):=\kappa_{\mu,d}(\mathbf{x},\mathbf{x})^{-1}$ is called the Christoffel function associated with $\mu$ and $d\in\mathbb{N}$. The following result states a fundamental extremal property of the Christoffel function.

Theorem 3.1** (see e.g. [Dunkl and Xu (2001), Nevai (1986)]).**

Let ${\boldmath\xi}\in\mathbb{R}^{p}$ be fixed, arbitrary. Then

[TABLE]

The Christoffel function plays an important role in orthogonal polynomials and the theory of interpolation and approximation, see e.g. [Szegö (1974), Dunkl and Xu (2001)]. One is particularly interested in the asymptotics of the normalized Christoffel function $\mathbf{x}\mapsto s(d)\Lambda_{\mu,d}(\mathbf{x})$ as $d\to\infty$. The subject has a very long history in the univariate case, see [Nevai (1986)] for a detailed historical account prior to the 80’s. The first quantitative asymptotic result was given in [Máté and Nevai (1980)] and was latter improved by [Máté et al. (1991)] and [Totik (2000)]. In the multivariate setting, precise results are known in some particular cases such as balls, spheres and simplices [Bos (1994), Bos et al. (1998), Xu (1996), Xu (1999), Kroo and Lubinsky (2013)] but much remains to be done for the general multivariate case. A typical example of asymptotic result is given under quite general (and technical) conditions in [Kroo and Lubinsky (2013), Kroo and Lubinsky (2012)]. This work shows that, as $d\to\infty$, the limit of the ratio, $\frac{\Lambda_{\nu,d}}{\Lambda_{\mu,d}}$, of two Christoffel functions associated to two mutually absolutely continuous measures $\mu$ and $\nu$, converges to the density $\frac{d\nu}{d\mu}(\mathbf{x})$ on the interior of their common support.

Remark 3.2**.**

Notice that Theorem 3.1 also provides a method to compute the numerical value $\Lambda_{\mu,d}(\mathbf{x})$ for $\mathbf{x}\in\mathbb{R}^{n}$, fixed, arbitrary. Then indeed (3.2) is a convex quadratic programming problem which can be solved efficiently, even in high dimension using first order methods such as projected gradient descent and its stochastic variants. This is particularly interesting when the nonsingular moment matrix $\mathbf{M}_{d}(\mu)$ is large.

We next provide additional insights on Theorem 3.1 and solutions of (3.2).

Theorem 3.3**.**

For any ${\boldmath\xi}\in\mathbb{R}^{p}$, the optimization problem in (3.2) is convex with a unique optimal solution $P^{*}_{d}\in\mathbb{R}[X]_{d}$ defined by

[TABLE]

In addition,

[TABLE]

The proof is postponed to Section 6. Interestingly, each of the orthonormal polynomials $(P_{\alpha})_{\alpha\in\mathbb{N}^{p}}$ also satisfies an important and well-known extremality property.

Theorem 3.4** (see e.g. [Dunkl and Xu (2001)]).**

Let $\alpha\in\mathbb{N}^{p}$ be fixed, arbitrary and let $d=|\alpha|$. Then up to a multiplicative positive constant, $P_{\alpha}$ is the unique optimal solution of

[TABLE]

Finally, we highlight the following important property which will be useful in the sequel.

Theorem 3.5** (See e.g. [Lasserre and Pauwels (2016)]).**

$\Lambda_{\mu,d}$* is invariant by change of polynomial basis $\mathbf{v}_{d}$, change of the origin of $\mathbb{R}^{p}$ or change of basis in $\mathbb{R}^{p}$.*

Remark 3.6**.**

All these statements can be deduced from identity (3.1). Indeed, we have, for any $\mathbf{x},\mathbf{y}\in\mathbb{R}^{p}$,

[TABLE]

for any invertible matrix $A$ of suitable size. All the proposed transformations induce a change of basis of polynomials up to degree $d$ which can be represented by such an $A$.

3.2 When $\mu$ is the Lebesgue measure

In this section we consider the important case of the Lebesgue measure on a compact set $S\subset\mathbb{R}^{p}$ such that ${\rm cl}({\rm int}(S))=S$. It is known that in this case the Christoffel function encodes information on the set $S$; see for example the discussion in Section 3.1. In particular, the scaled Christoffel function remains positive on the interior of $S$. We push this idea further and present a new result asserting that it is possible to recover the set $S$ with strong asymptotic guaranties by carefully thresholding the corresponding scaled Christoffel function.

For any measurable set $A$, denote by $\mu_{A}$ the uniform probability measure on $A$, that is $\mu_{A}=\lambda_{A}/\lambda(A)$ where $\lambda$ is the Lebesgue measure and $\lambda_{A}$ the measure consisting of the restriction of Lebesgue measure to $A$ which is defined by $\lambda_{A}(A^{\prime})=\lambda_{A\cap A^{\prime}}$ for any measurable set $A^{\prime}$.

Threshold and asymptotics

The main idea is to use quantitative lower bounds on the scaled Christoffel function, $s(d)\Lambda_{\mu_{S},d}$, on the interior of $S$ (Lemma 6.2) and upper bounds outside $S$ (Lemma 6.6). Recall that $\mu_{S}$ denotes the uniform measure on $S$. In combining these bounds one proves the existence of a sequence of thresholds of the scaled Christoffel function which estimate $S$ in a strongly consistent manner. Let us introduce the following notation and assumption.

Assumption 3.7**.**

• (a)

$S\subset\mathbb{R}^{p}$* is a compact set such that ${\rm cl}({\rm int}(S))=S$.*

• (b)

The sequence $\left(\delta_{k}\right)_{k\in\mathbb{N}}$ is a decreasing sequence of positive numbers converging to [math]. For every $k\in\mathbb{N}$, let $d_{k}$ be the smallest integer such that:

[TABLE]

where ${\rm diam}(S)$ denotes the diameter of the set $S$, and

[TABLE]

Remark 3.8** (On Assumption 3.7).**

• •

$d_{k}$* is well defined. Indeed, since $\delta_{k}$ is positive, the left hand side of (3.7) goes to [math] as $k\to\infty$ while the right hand side remains bounded for increasing values of $d_{k}$.*

• •

From the definition of $d_{k}$ and the fact that $\delta_{k}$ is decreasing, the sequence $\{d_{k}\}_{k\in\mathbb{N}}$ is non decreasing. Indeed, in (3.7) the right hand side is an increasing function of $\delta_{k}$ while the left hand side is decreasing so that if (3.7) is satisfied for a certain value of $d_{k}$ and $\delta_{k}$, it is also satisfied with the same $d_{k}$ and any value of $\delta\geq\delta_{k}$.

• •

Given $\left\{\delta_{k}\right\}_{k\in\mathbb{N}}$, computing $d_{k}$ can be done recursively and only requires the knowledge of ${\rm diam}(S)$ and $\lambda(S)$.

• •

A similar condition can be enforced if only upper bounds on ${\rm diam}(S)$ and on $\lambda(S)$ are available. In this case, replace these quantities by their upper bounds in (3.7) to obtain a similar result.

We are now ready to state the first main result of this section whose proof is postponed to Section 6 for sake of clarity of exposition. Recall the definition of the Hausdorff distance $d_{H}(X,Y)$ between two subsets $X,Y$ of $\mathbb{R}^{p}$:

[TABLE]

Theorem 3.9**.**

Let $S\subset\mathbb{R}^{p}$, $\{\delta_{k}\}_{k\in\mathbb{N}}$, $\{\alpha_{k}\}_{k\in\mathbb{N}}$ and $\{d_{k}\}_{k\in\mathbb{N}}$ satisfy Assumption 3.7. For every $k\in\mathbb{N}$ let $S_{k}\subset\mathbb{R}^{p}$ be the set defined by,

[TABLE]

Then, as $k\to\infty$,

[TABLE]

Remark 3.10**.**

The relevance of Hausdorff distance and the notion of distance between topological boundaries is discussed in [Cuevas et al. (2006)] and [Singh et al. (2009)].

Extension to more general probability measures

Theorem 3.9 can easily be extended to probability measures that are more general than uniform distributions, in which case we consider the following alternative assumption.

Assumption 3.11**.**

• (a)

$S\subset\mathbb{R}^{p}$* is a compact set such that ${\rm cl}({\rm int}(S))=S$.*

• (b)

The function $w\colon{\rm int}(S)\to[w_{-},+\infty)$ is integrable on ${\rm int}(S)$ with $w_{-}>0$. The measure $\mu$ is such that for any measurable set $A$, $\mu(A)=\int_{A\cap S}w(\mathbf{x})d\mathbf{x}$ and $\mu(S)=1$. $\left\{\delta_{k}\right\}_{k\in\mathbb{N}}$ is a decreasing sequence of positive numbers which converges to [math]. For every $k\in\mathbb{N}$, let $d_{k}$ be the smallest integer such that:

[TABLE]

where

[TABLE]

Under Assumption 3.11 we obtain the following analogue of Theorem 3.9.

Theorem 3.12**.**

Let $S\subset\mathbb{R}^{p}$, $w\colon S\to[w_{-},+\infty)$, $\{\delta_{k}\}_{k\in\mathbb{N}}$, $\{\alpha_{k}\}_{k\in\mathbb{N}}$ and $\{d_{k}\}_{k\in\mathbb{N}}$ satisfy Assumption 3.11. For every $k\in\mathbb{N}$, let

[TABLE]

Then, as $k\to\infty$,

[TABLE]

The proof is postponed to Section 6.

3.3 Discrete approximation via the empirical Christoffel function

In this section $\mu$ is a probability measure on $\mathbb{R}^{p}$ with compact support $S$. We focus on the statistical setting where information on $\mu$ is available only through a sample of points drawn independently from the given distribution $\mu$. In this setting, for every $n\in\mathbb{N}$, let $\mu_{n}$ denote the empirical measure uniformly supported on an independent sample of $n$ points distributed according to $\mu$. It is worth emphasizing that in principle $\Lambda_{\mu_{n},d}$ is easy to compute and requires the inversion of a square matrix of size $s(d)$, see (3.1). Note that the definition in (3.1) can only be used if the empirical moment matrix, $\mathbf{M}_{d}(\mu_{n})$ is invertible which is the case almost surely if $\mathbf{M}_{d}(\mu)$ is invertible and $n$ is large enough. Alternatively, the numerical evaluation of $\Lambda_{\mu_{n},d}(\mathbf{x})$ at $\mathbf{x}\in\mathbb{R}^{n}$, fixed arbitrary, reduces to solving the convex quadratic programming problem (3.2), which can in principle be done efficiently even in high dimension, see Remark 3.2.

Our second main result is for fixed $d\in\mathbb{N}$ and relates the population Christoffel function $\Lambda_{\mu,d}$ and its empirical version $\Lambda_{\mu_{n},d}$, as $n$ increases. We proceed by distinguishing what happens far from $S$ and close to $S$. First, Lemma 6.8 ensures that both Christoffel functions associated with $\mu$ and $\mu_{n}$ vanish far from $S$ so that the influence of this region can be neglected. Second, when closer to $S$ one remains in a compact set and the strong law of large numbers applies.

Theorem 3.13**.**

Let $\mu$ be a probability measure on $\mathbb{R}^{p}$ with compact support. Let $\left\{X_{i}\right\}_{i\in\mathbb{N}}$ be a sequence of i.i.d. $\mathbb{R}^{p}$-valued random variables with common distribution $\mu$. For $n=1,2,\ldots$, define the (random) empirical probability measure $\mu_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$. Then, for every $d\in\mathbb{N}$, $d>0$, such that the moment matrix $\mathbf{M}_{d}(\mu)$ is invertible, it holds that

[TABLE]

Equivalently

[TABLE]

where $\|\cdot\|_{\infty}$ denotes the usual “sup-norm”.

A detailed proof can be found in Section 6. Theorem 3.13 is a strong result which states a highly desirable property, namely that almost surely with respect to the random draw of the sample, the (random) function $\Lambda_{\mu_{n},d}(\cdot)$ converges to $\Lambda_{\mu,d}(\cdot)$ uniformly in $\mathbf{x}\in\mathbb{R}^{p}$ as $n$ increases. Since we manipulate polynomials, it can be checked that [Cuevas et al. (2006), Theorem 1] for general level sets can be applied in the setting of Christoffel level set estimation. We get the following consequence in terms of consistency of the boundary of plugin estimates for Christoffel level sets.

Theorem 3.14**.**

Let $\mu$ be a probability measure on $\mathbb{R}^{p}$ with compact support. Let $\left\{X_{i}\right\}_{i\in\mathbb{N}}$ be a sequence of i.i.d. $\mathbb{R}^{p}$-valued random variables with common distribution $\mu$. For $n=1,2,\ldots$, define the (random) empirical probability measure $\mu_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$. Then, for every $d\in\mathbb{N}$, $d>0$, such that the moment matrix $\mathbf{M}_{d}(\mu)$ is invertible and any $c\in(0,\sup_{\mathbf{x}\in\mathbb{R}^{p}}\{\Lambda_{\mu,d}(\mathbf{x})\})$, as $n$ increases, it holds that

[TABLE]

where $L=\left\{\mathbf{x}\in\mathbb{R}^{p},\Lambda_{\mu,d}(x)\geq c\right\}$ and $L_{n}=\left\{\mathbf{x}\in\mathbb{R}^{p},\Lambda_{\mu_{n},d}(x)\geq c\right\}$.

4 Applications

4.1 Rationale

In this section we describe some applications for which properties of the Christoffel function prove to be very useful in a statistical context. We only consider the case of bounded support. A relevant property of the scaled Christoffel function is that it encodes information on the support and the density of a population measure $\mu$:

• •

[Máté et al. (1991)], [Totik (2000)] and [Kroo and Lubinsky (2012)] provide asymptotic results involving the density of the input measure.

• •

Theorem 3.9 provides asymptotic results related to the support of the input measure.

The support and density of a measure is of interest in many statistical applications. However, the aforementioned results are limited to population measures which are not accessible in a statistical setting. In the context of empirical Christoffel functions, Theorem 3.13 suggests that these properties still hold (at least in the limit of large number of samples) when one uses the empirical measure $\mu_{n}$ in place of the population measure $\mu$. Combining these ideas suggests to use of the empirical Chritoffel function in statistical applications such as (a) density estimation, (b) support inference or (c) outlier detection. This is illustrated on simulated and real world data and we compare the performance with well established methods for the same purpose. Finally, we also describe another application, namely inversion of affine shuffling, whose links with statistics are less clear.

All results presented in this section are mainly for illustrative purposes. In particular, the choice of the degree $d$ as a function of the sample size $n$ was done empirically and a precise quantitative analysis is a topic of future research beyond the scope of the present paper.

4.2 Density estimation

Most asymptotic results regarding the scaled Christoffel function suggest that the limiting behaviour involves the product of a boundary effect term and a density term. Hence if one knows both the Christoffel function and the boundary effect term, one has access to the density term. Unfortunately, this boundary term is only known in specific situations, the most typical example being the Euclidean ball. Hence, in the present state of knowledge, one of the following is assumed to hold true.

• •

The support of the population measure $\mu$ is $S$ and $\lim_{d\to\infty}s(d)\Lambda_{\lambda_{S},d}$ exists (possibly unknown).

• •

The support is unknown but contains a set $S$ with the same property as above. In this case, we consider the restriction of the population measure $\mu$ to $S$. Note that a sample from the restriction is easily obtained from a sample from $\mu$ by rejection.

In both cases, assuming that $\mu$ has a density $h$ on $S$, it is expected that the ratii $\frac{\Lambda_{\mu,d}}{\Lambda_{\lambda_{S},d}}$ or $\frac{s(d)\Lambda_{\mu,d}}{\lim_{d\to\infty}s(d)\Lambda_{\lambda_{S},d}}$, converge to $h$. An example of such a result in the univariate setting is the following.

Theorem 4.1** (Theorem 5 [Máté et al. (1991)]).**

Suppose that $\mu$ is supported on $[-1,1]$ with density $h\geq a>0$. Then for almost every $x\in[-1,1]$,

[TABLE]

Extensions include [Totik (2000)] for general support and [Kroo and Lubinsky (2012)] for the multivariate setting. Combining Theorems 4.1 and 3.13 suggest that the empirical Christoffel function can be used for density estimation. For illustration purposes, we set $\mu$ to be the restriction of a Gaussian to $[-1,1]$. We perform the following experiment for given $n,d\in\mathbb{N}$.

• •

Generate $x_{1},\ldots,x_{n}\in[-1,1]$ sampled independently from $\mu$.

• •

Compute and plot $x\to\frac{\Lambda_{\mu_{n},d}(x)}{\Lambda_{\lambda_{S},d}(x)}$. Note that $\Lambda_{\lambda_{S},d}$ is easily derived from the moments of the uniform distribution on $[-1,1]$.

The result is presented in Figure 1 and a comparison is given with a classical technique, kernel density estimation [Rosenblatt (1956), Parzen (1962)] with the Gaussian kernel. The result suggest that empirical Christoffel based density estimation is competitive with kernel density estimation in this setting. It is worth noticing how simple the methodology is with a single parameter to tune.

4.3 Support inference

Combining Theorems 3.9 and 3.13 suggest that one may recover the unknown support of a population measure $\mu$ from $n$ independant samples by thresholding the scaled empirical Christoffel function. In this section we set $\mu$ to be the uniform probability measure over a star shaped domain in $\mathbb{R}^{2}$ (see Figure 2). For different values of the degree $d$ and sample size $n$, we plot in Figure 2 the corresponding sample and the associated level set $\left\{\mathbf{x}\in\mathbb{R}^{2}:\,s(d)\Lambda_{\mu_{n},d}(\mathbf{x})=\alpha(\delta)\right\}$, where $\delta=0.5$ and $\alpha$ is given in Assumption 3.7.

The results displayed in Figure 2 show that for well chosen values of $d$ and with $\alpha$ as in Assumption 3.7, the support of the population measure is rather well approximated from a finite independent sample. The results even suggest that a careful tuning of the degree $d$ and the threshold level set $\alpha$ allows to approximate the support extremely well for larger sample sizes. Of course the degree $d$ should be chosen to avoid a form of over-fitting as the results suggest for small sample sizes and large values of $d$. A precise analysis of this phenomenon is a topic of future research.

4.4 Outlier detection

In [Lasserre and Pauwels (2016)] we suggested that the empirical Christoffel function could be used for the purpose of detecting outliers and the claim was supported by some numerical experiments. The rationale for this is that the empirical Christoffel function encodes information about the population density and outliers can be seen as samples from low density areas. We follow the same line and consider the network intrusion detection task described in [Williams et al. (2002)] based on the KDD cup 99 dataset [Lichman (2013)]. Following the pre-processing described in [Williams et al. (2002), Lasserre and Pauwels (2016)], we build up five datasets consisting of network connections represented by labeled vectors in $\mathbb{R}^{3}$ where each label indicates wether the connection was an attack or not.

[TABLE]

All the experiments described in this section are performed on the “other” dataset which is the most heterogeneous. The main task is to recover attacks from the collection of points in $\mathbb{R}^{3}$, ignoring the labels, and then compare the predictions with the ground thruth (given by the labels). We compare different methods, each of them assign a score to an individual, the higher the score, the more likely the individual is to be an outlier, or an attack. The metric that we use to compare different methods is the area under the Precision Recall curve (AUPR); see for example [Davis and Goadrich (2006)]. We compare three different methods, each of them producing a score reflecting some degree of outlyingness.

• •

Empirical Christoffel function.

• •

Kernel density estimation [Rosenblatt (1956), Parzen (1962)] with Gaussian kernel. The value of the density estimated at each datapoint is used as an outlyingness score.

• •

One-class SVM [Schölkopf et al. (2001)] with Gaussian kernel. The value of the estimated decision function at each datapoint is used as an outlyingness score. We used the implementation provided in the kernlab package [Zeileis et al. (2004)].

The first two methods involve only a single parameter while the last method requires two parameters to be tuned. The results are given in Table 1 and the corresponding curves can be found in Section 6.1. Table 1 suggest that one-class SVM and the empirical Christoffel perform similarly and clearly outperform the kernel density estimation approach. It is worth noticing here that the one-class SVM provides slightly better performances but requires a precise tuning of the second parameter.

4.5 Inversion of affine shuffling

This last application has fewer connections with statistics. Suppose that we are given two matrices $\mathbf{X}\in\mathbb{R}^{p\times n}$ and $\mathbf{X}^{\prime}\in\mathbb{R}^{p\times n}$. Furthermore, we know that there exists an invertible affine mapping $\mathcal{A}\colon\mathbb{R}^{p}\to\mathbb{R}^{p}$ such that after a potential permutation of the columns, $\mathcal{A}$ defines a bijection between the columns of $\mathbf{X}^{\prime}$ and those of $\mathbf{X}$. The problem is to recover the correspondence between the columns of $\mathbf{X}$ and the columns of $\mathbf{X}^{\prime}$, whence the name “affine shuffling inversion”. Note that the columns may be shuffled in an arbitrary way and therefore the matching problem is not trivial. In this setting we can use the affine invariance property of the Christoffel function described in Theorem 3.5. This is based on the two following observations.

• •

The Christoffel function only depends on the empirical moments and hence is not sensitive to reshufling of the columns.

• •

Working with the affine image amounts to perform a change of basis and a change of origin. By Theorem 3.5, the evaluation of the Christoffel function does not change.

This suggests the following procedure.

• •

Compute $\Lambda_{\mathbf{X},d}$ and $\Lambda_{\mathbf{X}^{\prime},d}$ the Christoffel functions associated to the columns of $\mathbf{X}$ and the columns of $\mathbf{X}^{\prime}$ respectively.

• •

Set $A\in\mathbb{R}^{n}$ to be the vector with $\Lambda_{\mathbf{X},d}(\mathbf{X}_{i})$ as $i$-th entry where $\mathbf{X}_{i}$ is the $i$-th column of $\mathbf{X}$. Set $A^{\prime}$ similarly.

• •

Match the values in $A$ to the values in $A^{\prime}$ according to their rank.

The proposed procedure defines a unique permutation between columns of $\mathbf{X}$ and columns of $\mathbf{X}^{\prime}$ when there are no ties in the vectors $A$ and $A^{\prime}$. In this case, Theorem 3.5 ensures that we have found the correct correspondance. In case of ties, the procedure does not allow to elicit completely the correspondance matching. Overall, the method is not garanteed to work but allows to treat simple cases easily. Investigating the robustness of this procedure to noise or to matching mispecification is the subject of future research.

An illustration is given in Figure 3 where a moon shaped cloud of points in $\mathbb{R}^{2}$ is deformed by an affine transformation and the matching between the points between the two clouds is recovered by matching the corresponding Christoffel function values. The correspondance between Christoffel function values is illustrated on a quantile quantile plot.

5 Conclusion

In this paper we have investigated the potential of the empirical Christoffel function for some applications in statistics and machine learning. This question led us to investigate its theoretical properties as well as potential paths toward applications, mostly in a statistical framework.

On the theoretical side, we proposed two main contributions. The first one provides an explicit thresholding scheme which allows to use the Christoffel function to recover the support of a measure with strong asymptotic guarantees. Although this property finds its root in the long history of results regarding asymptotic properties of the Christoffel function, we have provided a systematic way to tune the threshold and the degree to ensure strong convergence guarantees. The second main contribution relates the empirical Christoffel function to its population counterpart in the limit of large samples. This type of results is new and paves the way toward a much more precise understanding of relations between these two objects in a small sample setting.

On the practical side, we have illustrated the relevance of the Christoffel function as a practical tool in a machine learning context. In particular, simulations and experiments on real world data support our claim that the empirical Christoffel function is potentially useful for density estimation, support inference and outlier detection. Finally, in another application outside the statistical framework (detection of affine matching between two clouds of points), we have again illustrated the potential of the Christoffel function as a tool in shape recognition and shape comparison.

Both theory and applications suggest a broad research program. As already mentioned, an important issue is to quantify the deviation of the empirical Christofel function from its population counterpart in a finite sample setting. Results in this direction could have both theoretical and practical impacts and would to compare more accurately the performance of Christoffel-based approaches with state-of-the art methods. Furthermore, the use of the Christoffel function in a statistical framework raises questions specific to each application considered in this paper and will be the subject of future investigations. Finally, there are still important open questions on the Christoffel function itself and works in the line of [Berman (2009)] are of great interest to address applications in statistics.

Acknowledgements

The research of the first author was funded by the European Research Council (ERC) under the European’s Union Horizon 2020 research and innovation program (grant agreement 666981 TAMING).

6 Appendix

6.1 Precision recall curves from section 4.4

This section displays the curves from which the AUPR scores were measured in Section 4.4. Christoffel function and kernel density estimation are presented in Figure 4 and the one-class SVM is presented in Figure 5. A detailed discussion the experiment is given in Section 4.4.

6.2 Proof of Theorem 3.3

Proof.

In the optimization problem (3.2) the objective function $P\mapsto\int P^{2}d\mu$ is strongly convex in the vector of coefficients of $P$ because

[TABLE]

and therefore (3.2) reads $\min\{P^{T}\mathbf{M}_{d}(\mu)\,P:P^{T}\mathbf{v}_{d}(\xi)=1\}$, which is a convex optimization problem with a strongly convex objective function. Slater’s condition holds (only one linear equality constraint) and so the Karush-Kuhn-Tucker (KKT) -optimality conditions are both necessary and sufficient. At an optimality solution $P^{*}_{d}$ they read:

[TABLE]

for some scalar $\theta$. Multiplying by $(P^{*}_{d})^{T}$ yields

[TABLE]

Hence necessarily

[TABLE]

which is (3.3). Next, let $\mathbf{e}_{\alpha}\in\mathbb{R}^{s(d)}$ be the vector with null coordinates except the entry $\alpha$ which is $1$. From the definition of the moment matrix $\mathbf{M}_{d}(\mu)$,

[TABLE]

which is (3.5). In particular with $\alpha:=0$, we recover (3.4),

[TABLE]

∎

6.3 Proof of Theorems 3.9 and 3.12

6.3.1 Lower bound on the Christoffel function inside $S$

We will heavily rely on results from [Bos (1994)] (note that similar results could be obtained on the box, see for example [Xu (1995)]). In particular, we have the following result.

Lemma 6.1**.**

We have for any $d\geq 2$

[TABLE]

Proof.

Combining Lemma 2 in [Bos (1994)] and the last equation of the proof of Lemma 3 in [Bos (1994)], we have

[TABLE]

The result follows by using the expression given for $s(d)$ and simplifying factorial terms. ∎

From this result, we deduce the following bound.

Lemma 6.2**.**

Let $\delta>0$ and $\mathbf{x}\in S$ such that ${\rm dist}(\mathbf{x},\partial S)\geq\delta$. Then

[TABLE]

Proof.

Decompose the measure $\mu_{S}$ into the sum,

[TABLE]

Hence, by monotonicity of the Christoffel function with respect to addition and closure under multiplication by a positive term (this follows directly from Theorem 3.1), we have

[TABLE]

Next, by affine invariance of the Christoffel function (Theorem 3.5),

[TABLE]

where $\mathbf{B}$ is the unit Euclidean ball in $\mathbb{R}^{p}$. The result follows by combining (6.1), (6.2), Lemma 6.1 and the fact that $\frac{\lambda(\mathbf{B}_{\delta}(\mathbf{x}))}{\lambda(\mathbf{B})}=\delta^{p}$. ∎

6.3.2 Upper bound on the Christoffel function outside $S$

We next exhibit an upper bound on the Christoffel function outside of $S$. We first provide a useful quantitative refinement of the “Needle polynomial” introduced in [Kroo and Lubinsky (2012)].

Lemma 6.3**.**

For any $d\in\mathbb{N}$, $d>0$, and any $\delta\in(0,1)$, there exists a $p$-variate polynomial of degree $2d$, $q$, such tha

[TABLE]

Proof.

Let $r$ be the univariate polynomial of degree $2d$, defined by

[TABLE]

where $T_{d}$ is the Chebyshev polynomial of the first kind. We have

[TABLE]

Furthermore, for $t\in[-1,1]$, we have $0\leq 1+\delta^{2}-t^{2}\leq 1+\delta^{2}$. $T_{d}$ has absolute value less than $1$ on $[-1,1]$ and is inceasing on $[1,\infty)$ with $T_{d}(1)=1$, so for $t\in[-1,1]$,

[TABLE]

For $|t|\in[\delta,1]$, we have $\delta^{2}\leq 1+\delta^{2}-t^{2}\leq 1$, so

[TABLE]

Let us bound the last quantity. Recall that for $t\geq 1$, we have the following explicit expression

[TABLE]

We have $1+\delta^{2}+\sqrt{(1+\delta^{2})^{2}-1}\geq 1+\sqrt{2}\delta$, which leads to

[TABLE]

where we have used concavity of the $\log$ and the fact that $1+\sqrt{2}\geq 2$. It follows by combining (6.3), (6.4), (6.5) and (6.6), that $q\colon y\to r(\|y-x\|_{2})$ satisfies the claimed properties. ∎

We recall the following well known bound for the factorial taken from [Robbins (1955)].

Lemma 6.4** ([Robbins (1955)]).**

For any $n\in\mathbb{N}$, we have

[TABLE]

We deduce the following Lemma.

Lemma 6.5**.**

For any $d\in\mathbb{N}$, $d>0$, we have

[TABLE]

Proof.

This follows from a direct computation using Lemma 6.4.

[TABLE]

which proves the result. ∎

Combining the last two Lemma, we get the following bound on the Christoffel function.

Lemma 6.6**.**

Let $\mathbf{x}\not\in S$ and $\delta$ be such that ${\rm dist}(\mathbf{x},S)\geq\delta$. Then, for any $d\in\mathbb{N}$, $d>0$, we have

[TABLE]

Proof.

We may translate the origin of $\mathbb{R}^{p}$ at $\mathbf{x}$ and scale the coordinates by $\delta+{\rm diam}(S)$, this results in $\mathbf{x}=0$ and distance from $\mathbf{x}$ to $S$ is at most $\delta^{\prime}=\frac{\delta}{\delta+{\rm diam}(S)}\leq 1$. Furthermore, $S$ is contained in the unit Euclidean ball $\mathbf{B}$. Using invariance of the Christoffel function with respect to change of origin and change of basis in $\mathbb{R}^{p}$, (Theorem 3.5), this affine transformation does not change the value of the Christoffel function. Now the polynomial described in Lemma 6.3 provides an upper bound on the Christoffel function. Indeed for any $d^{\prime}\in\mathbb{N}$, we have

[TABLE]

where we have used $\delta^{\prime}\leq 1$ to obtain the last inequality. Combining Lemma 6.5 and (6.7), we obtain for any $d^{\prime}\in\mathbb{N}$

[TABLE]

Since in (6.8) $d^{\prime}\in\mathbb{N}$ was arbitrary, we obtain in particular

[TABLE]

The result follows by from (6.9) by setting $\delta^{\prime}=\frac{\delta}{\delta+{\rm diam}(S)}$. ∎

6.3.3 Proof of Theorem 3.9

Proof.

Let us first prove that $\lim_{k\to\infty}d_{H}(S,S_{k})=0$. We take care of both expressions in the definition of $d_{H}$ separately. Fix an arbitrary $k\in\mathbb{N}$, from Assumption 3.7 and Lemma 6.6, for any $\mathbf{x}\in\mathbb{R}^{p}$ such that ${\rm dist}(\mathbf{x},S)>\delta_{k}$,

[TABLE]

From this we deduce that $\mathbb{R}^{p}\setminus S_{k}\supseteq\left\{\mathbf{x}\in\mathbb{R}^{p}:{\rm dist}(\mathbf{x},S)>\delta_{k}\right\}$ and thus $S_{k}\subseteq\left\{\mathbf{x}\in\mathbb{R}^{p}:{\rm dist}(\mathbf{x},S)\leq\delta_{k}\right\}$. Since $k$ was arbitrary, for any $k\in\mathbb{N}$,

[TABLE]

Inequality (6.10) allows to take care of one term in the expression of $d_{H}$. Let us now consider the second term. We would like to show that

[TABLE]

Note that the supremum is attained in (6.11). We will prove this by contradiction, for the rest of the proof, $M$ denotes a fixed positive number which value can change between expressions. Suppose that (6.11) is false. This means that for each $k\in\mathbb{N}$ (up to a subsequence), we can find $\mathbf{x}_{k}\in S$ which satisfies

[TABLE]

Since $\mathbf{x}_{k}\in S$ and $S$ is compact, the sequence $(\mathbf{x}_{k})_{k\in\mathbb{N}}$ has an accumulation point $\bar{\mathbf{x}}\in S$, i.e., (up to a subsequence) $\mathbf{x}_{k}\to\bar{\mathbf{x}}$ as $k\to\infty$. Since ${\rm dist}(\cdot,S_{k})$ is a Lipschitz function, combining with (6.12), for every $k\in\mathbb{N}$ (up to a subsequence),

[TABLE]

We next show that (6.13) contradicts the assumption $S={\rm cl}({\rm int}(S))$. From now on, we discard terms not in the subsequence and assume that (6.13) holds for all $k\in\mathbb{N}$. Combining Lemma 6.2 and Assumption 3.7, for every $k\in\mathbb{N}$

[TABLE]

Since $S={\rm cl}({\rm int}(S))$ and $\bar{\mathbf{x}}\in S$, consider a sequence $\{\mathbf{y}_{l}\}_{l\in\mathbb{N}}\subset{\rm int}(S)$ such that $\mathbf{y}_{l}\to\bar{\mathbf{x}}$ as $l\to\infty$. Since $\mathbf{y}_{l}\in{\rm int}(S)$, we have ${\rm dist}(\mathbf{y}_{l},\partial S)>0$ for all $l$. Up to a rearrangment of the terms, we may assume that ${\rm dist}(\mathbf{y}_{l},\partial S)$ is decreasing and ${\rm dist}(\mathbf{y}_{0},\partial S)\geq\delta_{0}$. For all $l$, denote by $k_{l}$ the smallest integer such that ${\rm dist}(\mathbf{y}_{l},\partial S)\geq\delta_{k_{l}}$. We must have $k_{l}\to\infty$ and we can discard terms so that $k_{l}$ is a valid subsequence. We have constructed a subsequence $k_{l}$ such that for every $l\in\mathbb{N}$, $\mathbf{y}_{l}\in S_{k_{l}}$ and $\mathbf{y}_{l}\to\bar{\mathbf{x}}$. This is in contradiction with (6.13) and hence (6.11) must be true. Combining (6.10) and (6.11) we have that $\lim_{k\to\infty}d_{H}(S,S_{k})=0$.

Let us now prove that $\lim_{k\to\infty}d_{H}(\partial S,\partial S_{k})=0$, we begin with the term $\sup_{\mathbf{x}\in\partial S_{k}}{\rm dist}(\mathbf{x},\partial S)$. Fix an arbitrary $k\in\mathbb{N}$ and $\bar{\mathbf{x}}\in\partial S_{k}$. We will distinguish the cases $\bar{\mathbf{x}}\in S$ and $\bar{\mathbf{x}}\not\in S$. Assume first that $\bar{\mathbf{x}}\not\in S$. We deduce from (6.10), that

[TABLE]

Assume now that $\bar{\mathbf{x}}\in S$. If $\bar{\mathbf{x}}\in\partial S$, we have ${\rm dist}(\bar{\mathbf{x}},\partial S)=0$. Assume that $\bar{\mathbf{x}}\in{\rm int}(S)$. From (6.14), we have that $S\setminus S_{k}\subseteq\left\{\mathbf{x}\in S:{\rm dist}(\mathbf{x},\partial S)<\delta_{k}\right\}$ and hence ${\rm cl}(S\setminus S_{k})\subseteq\left\{\mathbf{x}\in S:{\rm dist}(\mathbf{x},\partial S)\leq\delta_{k}\right\}$. Since $\bar{\mathbf{x}}\in\partial S_{k}\cap{\rm int}(S)$, we have $\bar{\mathbf{x}}\in{\rm cl}(S\setminus S_{k})$ and hence ${\rm dist}(\bar{\mathbf{x}},\partial S)\leq\delta_{k}$. Combining the two cases $\bar{x}\in S$ and $\bar{x}\not\in S$, we have in any case that ${\rm dist}(\bar{\mathbf{x}},\partial S)\leq\delta_{k}$ and hence

[TABLE]

Let us now prove that

[TABLE]

First since $S$ is closed by asumption, the supremum is attained for each $k\in\mathbb{N}$. Assume that (6.17) does not hold, this means there exists a constant $M>0$, such that we can find $\mathbf{x}_{k}\in\partial S$, $k\in\mathbb{N}$ with ${\rm dist}(\mathbf{x}_{k},\partial S_{k})\geq M$. If $\mathbf{x}_{k}\not\in S_{k}$ infinitely often, then, we would have up to a subsequence $\mathbf{x}_{k}\in S$ and ${\rm dist}(\mathbf{x}_{k},S_{k})\geq M$. This is exactly (6.12) and we alredy proved that it cannot hold true. Hence, $\mathbf{x}_{k}\not\in S_{k}$ only finitely many times and we may assume by discarding finitely many terms that $\mathbf{x}_{k}\in S_{k}$ for all $k\in\mathbb{N}$. Let $\bar{\mathbf{x}}\in\partial S$ be an accumulation point of $(\mathbf{x}_{k})_{k\in\mathbb{N}}$. Since $\bar{\mathbf{x}}\in\partial S$, there exists $\bar{\mathbf{y}}\not\in S$ such that $0<{\rm dist}(\bar{\mathbf{y}},S)\leq\lim\inf_{k\to\infty}{\rm dist}(\mathbf{x}_{k},\partial S_{k})/2$. Since $\mathbf{x}_{k}\in S_{k}$ for all $k$ sufficiently large, we have $\bar{\mathbf{y}}\in S_{k}$ for all $k$ sufficiently large but the fact that $0<{\rm dist}(\bar{\mathbf{y}},S)$ contradicts (6.10). Hence (6.17) must hold true and the proof is complete. ∎

Remark 6.7** (Refinements).**

The proof of Theorem 3.9 is based on the following fact

[TABLE]

Depending on the regularity of the boundary $\partial S$ of $S$, it should be possible to get sharper bounds on the distance as a function of $\delta_{k}$. This should involve the dependency on $\delta$ of the function

[TABLE]

For example, if the boundary $\partial S$ has bounded curvature, this function is equal to $\delta$ for sufficiently small $\delta$. Another example, if $S\subset\mathbb{R}^{2}$ is the interior region of a non self intersecting continuous polygonal loop, then the function is of the order of $\frac{\delta}{\sin\left(\frac{\theta}{2}\right)}$, where $\theta$ is the smallest angle between two consecutive segments of the loop.

6.3.4 Proof of Theorem 3.12

Proof.

Lemma 6.2 holds with $\mu$ in place of $\mu_{S}$ and $w_{-}$ in place of $\frac{1}{\lambda(S)}$. Indeed, we have

[TABLE]

and the rest of the proof remains the same with different constants. Similarly, Lemma 6.6 holds with $\mu$ in place of $\mu_{S}$, indeed, the proof only uses the fact that $\mu_{S}$ is a probability measure supported on $S$ which is also true for $\mu$. The proof then is identical to that of Theorem 3.9 by reflecting the corresponding change in the constants. ∎

6.4 Proof of Theorem 3.13

6.4.1 A preliminary Lemma

Lemma 6.8**.**

Let $\mu$ be a probability measure supported on a compact set $S$. Then for every $d\in\mathbb{N}$, $d>0$, and every $\mathbf{x}\in\mathbb{R}^{p}$,

[TABLE]

Proof.

Set $\mathbf{y}={\rm proj}(\mathbf{x},{\rm conv}(S))$, that is $\|\mathbf{y}-\mathbf{x}\|\,=\,{\rm dist}(\mathbf{x},{\rm conv}(S))$ and:

[TABLE]

Consider the affine function

[TABLE]

For any $\mathbf{z}\in S$, we have

[TABLE]

where we have used Cauchy-Schwartz and triangular inequalities. Furthermore, we have for any $\mathbf{z}\in S$,

[TABLE]

where we have used equation (6.18). Consider the affine function $q_{\mathbf{x}}\colon\mathbf{z}\to 1-f_{\mathbf{x}}(\mathbf{z})$. We have

[TABLE]

where the inequalities are obtained by combining (6.20) and (6.21). The result follows from (6.18), (6.22) and Theorem 3.1. ∎

6.4.2 Proof of Theorem 3.13

.

Proof.

First let us consider measurability issues. Fix $n$ and $d$ such that $M_{d}(\mu)$ is invertible. Let $\mathbf{X}$ be a matrix in $\mathbb{R}^{p\times n}$, we use the shorthand notation

[TABLE]

where for each $i$, $\mathbf{X}_{i}$ is the $i$-th column of the matrix $\mathbf{X}$. This corresponds to the empirical Christoffel function with input data given by the columns of $\mathbf{X}$. Consider the function $F\colon\mathbb{R}^{p\times n}\to[0,1]$ defined as follows:

[TABLE]

It turns out that $F$ is a semi-algebraic function (its graph is a semi-algebraic set). Roughly speaking a set is semi-algebraic if it can be defined by finitely many polynomial inequalities. We refer the reader to [Coste (2000)] for an introduction to semi-algebraic geometry, we mostly rely on content from Chapter 2. First, the function

[TABLE]

is semi-algebraic (by identifying the space of polynomials with the Euclidean space of their coefficients) and the set $\left\{(P,\mathbf{z}):P(\mathbf{z})=1\right\}$ is also semi-algebraic. Constrained partial minimization can be expressed by a first order formula, and, by Tarski-Seidenberg Theorem (see e.g. [Coste (2000), Theorem 2.6]), this operation preserves semi-algebraicity. Hence, the function $(\mathbf{X},\mathbf{z})\to\Lambda_{\mathbf{X},d}(\mathbf{z})$ is semi-algebraic. Furthermore, Theorem 3.1 ensures that $\Lambda_{\mu,d}(\mathbf{z})=1/\kappa(\mathbf{z},\mathbf{z})$ for any $\mathbf{z}$, where $\kappa(\mathbf{z},\mathbf{z})$ is a polynomial in $\mathbf{z}$ and hence $\mathbf{z}\to\Lambda_{\mu,d}(\mathbf{z})$ is semi-algebraic. Finally absolute value is semi-algebraic and using a partial minimization argument again, we have that $F$ is a semi-algebraic function.

As a semi-algebraic function, $F$ is Borel measurable. Indeed, using the good sets principle ([Ash (1972)] §1.5.1, p. 35) it is sufficient to prove that, for an arbitrary interval333Recall that the Borel $\sigma$-field $\mathcal{B}([0,1])$ is generated by the intervals $(a,b]$ of $[0,1]$; see [Ash (1972)] §1.4.6, p. 27. $(a,b]\subset[0,1]$, $F^{-1}((a,b])\in\mathcal{B}(\mathbb{R}^{p\times n})$. Any such set is the pre-image of a semi-algebraic set by a semi-algebraic map. As proved in [Coste (2000), Corollary 2.9], any such set must be semi-algebraic and hence measurable. Thus, with the notations of Theorem 3.13, $\|\Lambda_{\mu_{n},d}-\Lambda_{\mu,d}\|_{\infty}$ is indeed a random variable for each fixed $n,d$ such that $M_{d}(\mu)$ is invertible.

We now turn to the proof of the main result of the Theorem. For simplicity we adopt the following notation for the rest of the proof. For any continuous function $f:\mathbb{R}^{p}\to\mathbb{R}$, and any subset $V\subseteq\mathbb{R}^{p}$,

[TABLE]

which could be infinite. We prove that for any $\epsilon>0$,

[TABLE]

where the probability is taken with respect to the random choice of the sequence of independent samples from $\mu$ and the limit supremum is the set theoretic limit of the underlying events.

Fix $\epsilon>0$. Denote by $S$ the compact support of $\mu$. Note that $S$ contains also the support of $\mu_{n}$ with probability one. From Lemma 6.8, we have an upper bound on both $\Lambda_{\mu_{n},d}$ and $\Lambda_{\mu,d}$ of order $O\left({\rm dist}(\mathbf{x},{\rm conv}(S))^{-2}\right)$ which holds with probability one. Hence, it is possible to find a compact set $V_{\epsilon}$ containing $S$ (with complement $V_{\epsilon}^{c}=\mathbb{R}^{n}\setminus V_{\epsilon}$) such that, almost surely,

[TABLE]

Next, we have the following equivalence

[TABLE]

On the other hand, since both functions are non negative, from equation (6.27), almost surely,

[TABLE]

Hence the second event in the right hand side of (6.30) occurs with probability zero. As a consequence, except for a set of events of measure zero, we have

[TABLE]

which in turn implies

[TABLE]

By assumption the moment matrix $\mathbf{M}_{d}(\mu)$ is invertible and by the strong law of large numbers, almost surely, $\mathbf{M}_{d}(\mu_{n})$ must be invertible for sufficiently large $n$. Assume that $\mathbf{M}_{d}(\mu_{n})$ is invertible, we have

[TABLE]

Using the strong law of large numbers again, continuity of eigenvalues and the fact that for large enough $n$, $\mathbf{M}_{d}(\mu_{n})$ is invertible with probability one, the continuous mapping theorem ensures that almost surely, for $n$ sufficiently large, the smallest eigenvalue of $\mathbf{M}_{d}(\mu_{n})^{-1}$ is close to that of $\mathbf{M}_{d}(\mu)^{-1}$ and hence bounded away from zero. Since the first coordinate of $\mathbf{v}_{d}(\mathbf{x})$ is $1$, the denominator in (6.33) is bounded away from zero almost surely for sufficiently large $n$. In addition, since $V_{\epsilon}$ is compact, $\mathbf{v}_{d}(\mathbf{x})$ is bounded on $V_{\epsilon}$ and there exists a constant $K$ such that, almost surely, for sufficiently large $n$,

[TABLE]

where the matrix-norm in the right hand side is the operator norm induced by the Euclidean norm. Combining (6.32) and (6.34), we obtain

[TABLE]

The strong law of large numbers and the continuity of the matrix inverse $\mathbf{M}_{d}(\cdot)^{-1}$ at $\mu$ ensure that the right hand side of (6.4.2) is [math]. This concludes the proof. ∎