Gilles Mordant

Stephen Zhang

Quadratically Regularized Optimal Transport: nearly optimal potentials and convergence of discrete Laplace operators

Keywords: Regularised Optimal Transport, rates of dual potentials, manifold learning, porous medium equation

 suggesting that optimal transport with quadratic regularisation can be used to construct a graph whose discrete Laplace operator converges to the Laplace-Beltrami operator. We derive first order optimal potentials for the problem under consideration and find that the resulting solutions exhibit a surprising resemblance to the well-known Barenblatt-Prattle solution of the porous medium equation. Then, relying on these first order optimal potentials, we derive the pointwise L 2 -limit of such discrete operators built from an i.i.d. random sample on a smooth compact manifold. Simulation results complementing the limiting distribution results are also presented.

Introduction and main results

1.1 Discrete optimal transport with quadratic regularisation.

In this entire paper, we will consider a d-dimensional compact smooth Riemannian manifold (M, g) isometrically embedded in R p via the embedding ι : M → R p . In the sequel, we denote by ι * , the differential of this embedding.

Consider {ι(x i )} N i=1 a set of points of the manifold embedded in the ambient Euclidean space. These can be random or deterministic. Then, let

C ij = C(x i , x j ) = 1 2 ∥ι(x i ) -ι(x j )∥ 2 2
be the matrix of pairwise distances. Let µ = N -1 i δ xi be the uniformly weighted empirical measure. The discrete quadratically regularised optimal transport (QOT) problem reads min π∈Π(µ,µ)

⟨C, π⟩ + ε 2 ∥π∥ 2 2 ⇔ min π∈Π(µ,µ) ∥π + ε -1 C∥ 2 2 , (1)
where Π(µ, µ) = {π : π1 = π ⊤ 1 = 1/N } denotes the set of bistochastic couplings. Following Lorenz, Manns, and Meyer (2021, Equation D), we have the following duality result in terms of the dual potential u.

sup u ⟨u, µ⟩ - 1 4ε ∥[u ⊕ u -C] + ∥ 2 2 , (2)
where [x] + is the positive part of x. Denoting by u ⋆ the optimal solution of (2) the relationship at optimality between the primal and dual variables is given by

π ⋆ = [u ⋆ ⊕ u ⋆ -C] + ε . ( 3
)
As it is an optimal transport plan belonging to Π(µ, µ), we recall the constraints

N i=1 π ⋆ i,j = 1 N , ∀j ∈ {1, . . . , N }, (4)
which are crucial to understand the optimal dual potentials.

1.2 Discrete operators based on QOT.

A consequence of the formula (3) is that the optimal transport plan is sparse -entries π ⋆ ij are identically zero whenever C ij becomes too large. Because of this, [START_REF] Matsumoto | Beyond knn: Adaptive, sparse neighborhood graphs via optimal transport[END_REF] proposed to use the (rescaled) optimal transport plan as the adjacency matrix of a graph between pairs of points. The resulting weighted, undirected graph can in turn be used for downstream applications including semi-supervised learning, manifold learning, or dimensionality reduction in single cell RNA sequencing applications. The authors found that remarkable performance was achieved in these examples. More precisely, the weight matrix W that they consider is given by

W i,j := π ⋆ i,j N j=1 π ⋆ i,j
.

(5) [START_REF] Matsumoto | Beyond knn: Adaptive, sparse neighborhood graphs via optimal transport[END_REF] further raise the question whether a discrete Laplace operator constructed from this matrix converges to a Laplace-Beltrami operator in the limit of infinitely many samples. This motivates the present paper.

Main contributions and outline

Our main contribution is twofold. First, we establish the correct asymptotic order of the potentials as a function of the regularisation parameter ε and the sample size N in the discrete setting. This constitutes the content of Sections 2 and 3. Interestingly we find that the rates match with the solution of the porous medium equation, which is believed to be linked to optimal transport with a quadratic regularisation. Then, in Section 4, we then prove that, under suitable conditions, the discrete operator indeed can converge to the Laplace-Beltrami operator for random samples, as stated in Theorem 1.

Manifold setting and notation

Let us now describe a bit more the manifold setting that we consider for random samples. Let X be a p-dimensional random variable whose range is supported on M. Let us further assume for simplicity that X has a uniform distribution on the manifold, i.e., the density dP (x) = vol -1 (M)dx, ∀x ∈ ι(M).

Whenever we write E or Var, unless otherwise denoted we mean to be with respect to P . Let us denote by s(x), the scalar curvature of the manifold at x and by I x the second fundamental form of the isometric embedding ι at x. In the sequel, ∇ will denote the covariant derivative while ∆ will be the Laplace-Beltrami operator. Further, set

ω(x) = 1 |S d-1 | S d-1 ∥I x (θ, θ)∥ 2 dθ,
as well as

N(x) = 1 |S d-1 | S d-1 I x (θ, θ)dθ,
where S d-1 denotes the (d -1)-dimensional unit sphere in T x M.

For the sake of simplicity, let us make the following assumptions.

Assumption 1. The manifold M is correctly shifted and rotated so that ι * T x0 M is spanned by e 1 , . . . , e d .

Assumption 2. The manifold is properly rotated and translated so that e d+1 , . . . , e p diagonalise the second fundamental form I x0 .

These two conditions are not particularly important, they just help simplify both notation and result statements. Under these assumptions, we use the notation v 1 , v 2 for the vector v whose d first components are the vector v 1 and its p -d last components are the vector v 2 . The p × r matrix Jp,r is then defined as Jp,r := 0 p-r×r

I r×r

In what follows, for some η sufficiently small, define the fattened manifold in the embedded space by N , i.e. the set N := {x ∈ R p : inf y∈M ∥x -ι(y)∥ ≤ η}. We recall that the Laplace operator for a function g on a certain embedded smooth manifold ι(M) is defined for p ∈ ι(M) as

∆g(p) := (∆ R p g ext )(p)
with ∆ R p the usual Laplace operator in the Euclidean space, g ext (x) := g(π M (x)), x ∈ N , and π M projects N onto the manifold.

2 Ansatz for the potentials in the continuous case

An important question is understanding the relationship between the optimal potential u ⋆ and the chosen regularisation parameter ε. For x 0 ∈ M, it holds for r sufficiently small 1 that

E [f (X; r)1{∥X -ι(x 0 )∥ ≤ r}] = |S d-1 | d vol(M) f (ι(x 0 ); r)r d + |S d-1 | d(d + 2) 1 2 vol(M) ∆f (x 0 ; r) + s(x 0 )f (x 0 ; r) 6 vol(M) + d(d + 2)ω(x 0 )f (x 0 ; r) 24 vol(M) r d+2 (6)
+ O(f (x 0 ; r)r d+3 ).

This a slight variation of Lemma B.5 of [START_REF] Wu | Think globally, fit locally under the manifold setup: Asymptotic analysis of locally linear embedding[END_REF] for a uniform density and in which the functions are allowed to depend on the parameter r. Because of this modification the asymptotic expansion has been slightly refined.

Consider the continuous setting of the problem (1): following Lorenz, Manns, and Meyer (2021) we write π ∈ L 2 (M) to be the density of a candidate transport plan w.r.t. product measure on M × M, i.e. π(x, y)dxdy = 1. Then, the optimal transport plan π ⋆ in the quadratically regularised problem must satisfy

1 vol(M) M π ⋆ (ι(x 0 ), y)dy = E [π ⋆ (ι(x 0 ), X)] = 1 vol(M) 2 .

Taking the relation π

⋆ = ε -1 [u ⋆ ⊕ u ⋆ -C] + where u ⋆ ∈ L 2 (M)
is the corresponding optimal dual potential, making the ansatz that u ⋆ ∼ ε α and invoking (6) we have, for ε ≪ 1,

E u ⋆ (ι(x 0 )) + u ⋆ (X) - 1 2 ∥ι(x 0 ) -X∥ 2 2 + ∼ E ε α -∥X -ι(x 0 )∥ 2 1{∥X -ι(x 0 )∥ ≤ ε α/2 } ∼ |S d-1 | d ε α(d+1) 2 + O ε α(d+2) 2
, (in the above multiplicative constants were dropped). The above quantity must behave asymptotically like ε at leading order, and so matching exponents gives us

α + dα 2 = 1 ⇔ α = 2 2 + d .
In the continuous case, the optimal potentials must thus behave in the first order like ε 2 2+d .

1 For the entire paper, "r sufficiently small" must be understood as r being smaller than the injectivity radius of the manifold.

There is a belief in the community that there should be some link between quadratically regularised optimal transport and a class of nonlinear partial differential equations known as the porous medium equation on R d for index m = 2 (see e.g. [START_REF] Lavenant | Dynamical optimal transport on discrete surfaces[END_REF]), i.e., the equation

∂u ∂t = ∆(u m ),
where u = u(x, t) and with an initial condition on u at time t = 0. Starting from a Dirac mass of integral m at the origin, the solution of the porous medium equation for m = 2 is given by the Barenblatt-Prattle formula [START_REF] Vázquez | The porous medium equation: mathematical theory[END_REF]:

u(x, t) = max 0, t -d 2+d m - 1 4(d + 2) ∥x∥ 2 t 2 2+d
. or, rewriting terms,

u(x, t) = max 0, t -1 m t 2 2+d - 1 4(d + 2) ∥x∥ 2 .
As the porous medium equation conserves mass, the integral of u(x, t) over R d is m. A key property of the porous medium equation which distinguishes it from the standard diffusion equation is that the solution remains compactly supported. This is a property that also applies to the transport plans derived from quadratically regularized optimal transport [START_REF] Lorenz | Quadratically regularized optimal transport[END_REF]. Perhaps closer to the theory of optimal transport, the porous medium equation of index m can also be understood as the 2-Wasserstein gradient flow of the Tsallis entropy of order m (see for example, the discussion in [START_REF] Peyré | Entropic approximation of Wasserstein gradient flows[END_REF]). The Tsallis entropy generalizes the Gibbs entropy: for m = 1, it coincides with the Gibbs entropy, while for m = 2 it is corresponds to the squared L 2 norm of the density. It is remarkable that the squared L 2 norm is the functional that generates the porous medium equation as Wasserstein gradient flow, which is also the regularizing functional used in quadratically regularized optimal transport exhibiting analogous sparsity and scaling behaviour. Furthermore in the entropy regularized setting where m = 1, optimal transport enjoys the celebrated connection to a theory of large deviations for Brownian motions and the Schrödinger problem [START_REF] Léonard | A survey of the Schrödinger problem and some of its connections with optimal transport[END_REF]. One interesting theoretical question would be whether similar connections could hold in more general cases, e.g. m > 2.

Although our work does not formally establish the existence of such a connection, it is interesting that the same types of exponents appear and that the solution of the porous medium equation is so close in form to the solution of the quadratically regularised optimal transport problem.

Optimal potentials in the discrete case

We now turn the to study of the optimal rates for the potentials in the discrete case. Let us slightly change the setting and consider a sample of size N + 1 where one point, x 0 ∈ M, is fixed and the remaining ones are an i.i.d. random sample on the manifold. Set X 0 = ι(x 0 ). Relabel the sample points so that ∥X

0 -X 1 ∥ 2 ≤ ∥X 0 -X 2 ∥ 2 ≤ . . . ≤ ∥X 0 -X N ∥ 2 .
Remark 1 (No loss of generality in choosing X 0 ). In the results below, the same analysis has to be carried out for each point X i . One can thus view our (convenient) choice of working with one distinguished, deterministic point as a conditioning on an arbitrary X i . Still, as the expectations of the quantities for X 0 fixed are constants with uniformly decaying terms, the reasoning would apply for each X i using the tower property of conditional expectation.

Using Lemma B.5. with f = 1, we get that the quantile function of the local distribution of squared distances at x 0 is approximately

p → p |S d-1 |d -1 vol(M) -1 2/d
, so that the duality constraint in the discrete problem is approximately

N + 1 ε N j=0 u ⋆ (X 0 ) + u ⋆ (X j ) - U (j:N ) |S d-1 |d -1 vol(M) -1 2/d + = 1, (7)
where U (j:N ) is the j-th sorted element of an i.i.d. sample of size N of random variables uniformly distributed on [0, 1] and we take U (0:N ) = 0. Note that although we write down all N order statistics U (j:N ) and the expression for the quantile function is only a good approximation for p ≪ 1, as long as u(x 0 ) + u(X j ) is small, only the first few terms will be nonzero. We get the equivalent problem

N + 1 ε(N + 1) 2/d κ d N j=0 ũ(x 0 ) + ũ(X j ) -U 2/d (j:N ) + = 1, (8)
where we have set

ũ(•) = (N + 1) 2/d κ d u ⋆ (•),
and

κ d := vol(M)d |S d-1 | 2/d . Choose k 2/d < 2ũ ≤ (k + 1) 2/d
and plug ũ (choosing a constant approximation to the potential) as a choice for the potential in (8). It yields,

N + 1 ε(N + 1) 2/d κ d   2(k + 1)ũ - k j=1 U 2/d (j:N )   = 1. ( 9
)
Let us turn to the size of the sum in (9). First, basic calculations show that

E U (j:N ) 2/d = Γ( 2 d + j)Γ(N + 1) Γ( 2 d + N + 1)Γ(j)
, So that understanding the problem (7), even in expectation and for constant potentials is not so easy for d > 2.

For d = 1, we get

k j=1 E U (j:N ) 2/d = 1 (N + 2)(N + 1) 1 3 k(k + 1)(k + 1) , while for d = 2, it holds that k j=1 E U (j:N ) 2/d = 1 (N + 1) 1 2 k(k + 1) .
In general,

k j=1 E U (j:N ) 2/d ≈ (N + 1) -2/d k j=1 j 2/d 1 + 2 -d d 2 j + O 1 d 2 ,
so that the leading order is

(N + 1) -2/d d d + 2 k d+2 d .
Equation ( 9) then becomes,

εN 2/d-1 κ d ≈ 1 - d d + 2 k d+2 d = 2 d + 2 k d+2 d so that ũ ≈ 1 2 2 + d 2κ d 2 d+2 ε 2 d+2 N (2-d)2 d(d+2)
and then

u ≈ κ d d+2 d d + 2 2 2 d+2 ε 2 d+2 N (2-d)2 d(d+2) -2/d = vol(M)d |S d-1 | (d + 2) 2 2 d+2 ε 2 d+2 N -4 (d+2)
Finally, using this first order to approximately solve the equation yields that the optimal potential must behave like

K ε,N := 2κ d d + 2 2 d+2 ε 2 d+2 N -4 (d+2) =: C d ε 2 d+2 N -4 (d+2)
These results provide a reasonable ansatz, still these are only approximations. We thus now assess the quality of this first order approximation of the solution by evaluating how the dual constraints are fulfilled when plugging-in the first-order approximation of the solution.

Validity of the derived finite sample rate

We will use the function

f (y) = C d ε 2/(d+2) N -4/(d+2) -∥ι(x 0 ) -y∥ 2 = K ε,N -∥ι(x 0 ) -y∥ 2
and apply to the result of Lemma B.5 from [START_REF] Wu | Think globally, fit locally under the manifold setup: Asymptotic analysis of locally linear embedding[END_REF] to evaluate the constraints arising from the dual formulation of the problem, recall (3) and (4). Doing so, one gets

E   N + 1 ε N j=0 K ε,N -∥X j -ι(x 0 )∥ 2 +   = N + 1 ε K ε,N + N E f (X)1 ∥ι(x 0 ) -X∥ ≤ K 1/2 ε,N = N + 1 ε K ε,N + N (N + 1) ε |S d-1 | d vol(M) K 1+d/2 ε,N + N (N + 1) ε |S d-1 | d(d + 2) vol(M) K 1+d/2 ε,N -d + s(x 0 )K ε,N 6 + d(d + 2)ω(x 0 )K ε,N 24 + N (N + 1) ε O(K 5+d 2 ε,N ) = N + 1 ε K ε,N + N (N + 1) ε |S d-1 | vol(M) K 1+d/2 ε,N 1 d - 1 d + 2 + O N (N + 1) ε K 2+d/2 ε,N
In the display above, terms have orders

O(N ε -1 K ε,N ) = O(ε -d/(d+2) N (d-2)/(d+2) ), O(N 2 ε -1 K 2+d/2 ε,N
) and we remark that

O(N 2 ε -1 K 1+d/2 ε,N ) = O(1)
. This latter term is the leading order. Then, to fulfill the constraint, we need that

C 2+d 2 d |S d-1 | vol(M) 1 d - 1 d + 2 = 1 so that C d = vol(M) |S d-1 | d(d + 2) 2 2 d+2
, which matches with the expression above. One gets that the chosen rate gives the correct constraint in expectation at the first order.

One can also rewrite the conditions that O(ε

-d/(d+2) N (d-2)/(d+2) ) = o(1) as ε -d N (d-2) → 0 and the condition O(N 2 ε -1 K 2+d/2 ε,N ) = O(K ε,N ) = o(1) as ε 2 N -4 → 0.
Note that the latter condition was already somewhat required to apply Lemma B.5 from [START_REF] Wu | Think globally, fit locally under the manifold setup: Asymptotic analysis of locally linear embedding[END_REF]. Together, these results indicate that the asymptotic scaling on ε is N 1-2/d ≪ ε ≪ N 2 for the constraints to be asymptotically fulfilled in expectation. Note that in this analysis that ε need not go to zero asymptotically. Rather, the need is for ε to be asymptotically sufficiently small relative to N . See the remark below.

Remark 2. In the developments above, the result holds uniformly in x 0 under quite mild assumptions as, for a closed 2 and smooth manifold, the different kinds of curvatures appearing in the expansions are bounded, recall Remark 1.

Because of Remark 2, one can derive that all the constraints will asymptotically be fulfilled in expectation when replacing the sum of optimal potentials by K ε,N .

Remark 3. There is a difference in scaling between the discrete and continuous settings in our analysisto get empirical input distributions that are consistent with the continuous setting in the limit of large N , for samples X 1 , . . . , X N we take μ = N -1 i δ Xi as the corresponding empirical distribution. Suppose π ij is an admissible coupling for such a discrete problem. Then π is concentrated on the support of μ ⊗ μ and admits a density, (dπ/dμ ⊗ dμ)(X i , X j ) = π ij /N 2 . Then, note that the corresponding empirical entropy term would behave like

H(π|μ ⊗ μ) = dπ log dπ dμ ⊗ dμ = ij π ij log π ij N 2
Up to a constant, this is equal to the discrete entropy of π, i.e. ij π ij log π ij . Thus, we expect no scaling behaviour between N and the entropic regularizer.

On the other hand, for the quadratic regularizer, one would have

∥π∥ 2 μ⊗μ = dπ dμ ⊗ dμ dπ = ij π ij π ij N 2 = N -2 ∥π∥ 2 .
Thus, there is the presence of a factor N -2 . This can be understood in that N appears in the density of π w.r.t. empirical product measure, which is lost as an additive term in the case of a log, but remains in the quadratic case. Thus, noting that ε∥π∥ 2 μ⊗μ = (εN -2 )∥π∥ 2 2 , it is apparent that the requirement that ε → 0 in the continuous setting corresponds to εN -2 → 0 in the discrete setting. This is in agreement with the scaling we derived earlier.

Replacing the optimum by a uniform approximation

Similarly to Lorenz, Manns, and Meyer (2021, Lemma 3.1), the Newton Hessian of the optimisation problem (2) is diag(σ1 N +1 )

where

σ ij = 1 if u i + u j -C ij ≥ 0, 0 otherwise.
Because of the constraints, the potential must be chosen such that u j > 0, ∀j ≤ N + 1. The function to optimise is thus strictly concave for the set of such potentials and thus admits a unique optimum. The update step in the semismooth Newton algorithm used in Matsumoto, Zhang, and Schiebinger (2022, Appendix A)3 and originally developed in [START_REF] Lorenz | Quadratically regularized optimal transport[END_REF] takes the form

G -1 i,i   j (K ε,N -c i,j ) + - ε N + 1   , 1 ≤ i ≤ N,
when the regularisation parameter of the algorithm is set to zero. From (6) again, it holds that

EG ii = O(N K d/2 ε,N ).
From this and the computations of Section 3.1, the expectation of the update step is of order

1 N K d/2 ε,N O(K ε,N ) + O N K 2+d/2 ε,N = O ε 2-d d+2 N -2(2-d)-(d+2) d+2 + O(K 2 ε,N ) = O ε 2-d d+2 N d-6 d+2 + O(K 2 ε,N ) = O ε -d d+2 N d-2 d+2 ε 2 d+2 N -4 d+2 + O(K 2 ε,N ),
which goes to zero in view of the conditions on ε and N mentioned above. Remark 4 (QOT, nonparametric statistics and optimality). The form of the optimal transport plan is very much alike an Epanechnikov kernel, which is very often used in nonparametric statistics. The latter kernel is

u → Γ(2 + d/2) π d 2 (1 -u ⊤ u)1 {u ⊤ u≤1} .
Even though this statement is debated [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF], Section 1.2.4), the Epanechnikov kernel is often claimed to be the optimal nonnegative kernel in terms of asymptotic MISE for the estimation of a twice differentiable density. Thus, the compact support and the fact that the optimal dual potential is a function-which is likely more adaptive to the data than a uniform bandwidth, might explain the outstanding performances observed in the examples of [START_REF] Matsumoto | Beyond knn: Adaptive, sparse neighborhood graphs via optimal transport[END_REF].

4 Graphs Laplacians based on Quadratically Regularized OT

Limiting operators

Before proving the main result, we state two useful lemmas.

Lemma 1. For fixed x 0 ∈ M, r sufficiently small and X uniformly distributed on M under Assumptions 1 and 2, it holds that

E f (X; r)(X -ι(x 0 ))(X -ι(x 0 )) ⊤ 1{∥X -ι(x 0 )∥ ≤ r} = |S d-1 | d(d+2) vol(M) f (ι(x 0 ); r)r d+2 I d×d 0 0 0 + O(r 2 )
Proof. The proof follows along the same lines as Proposition 3.1 in [START_REF] Wu | Think globally, fit locally under the manifold setup: Asymptotic analysis of locally linear embedding[END_REF].

Lemma 2. For fixed x 0 ∈ M, r sufficiently small and X uniformly distributed on M under Assumptions 1 and 2, it holds that

E f (X; r)e ⊤ k (X -ι(x 0 ))(X -ι(x 0 )) ⊤ e l e ⊤ m (X -ι(x 0 ))(X -ι(x 0 )) ⊤ e n 1{∥X -ι(x 0 )∥ ≤ r} = f (ι(x 0 ); r) (d + 4) vol(M) r d+4 C k,l,m,n + O(r d+5 ), where C k,l,m,n = S d-1 ⟨ι * θ, e k ⟩⟨ι * θ, e l ⟩⟨ι * θ, e m ⟩⟨ι * θ, e n ⟩dθ Proof. First set Br (x 0 ) := ι -1 (B R p r (ι(x 0 )) ∩ ι(M)).
Then, the quantity of interest can be written

I := 1 vol(M) Br(x0)
⟨ι(y) -ι(x 0 ), e k ⟩⟨ι(y) -ι(x 0 ), e l ⟩⟨ι(y) -ι(x 0 ), e m ⟩⟨ι(y) -ι(x 0 ), e n ⟩f (y)dV (y).

Recalling that for (t, θ)

∈ [0, ∞) × S d-1 ι • exp x0 (θt) -ι(x 0 ) = ι * θt + O(t 2 ) r = r + O(r 3 ) dV (exp x0 (θt)) = t d-1 + O(t d+1 ) f (exp x0 (θt)) = f (x 0 ) + O(t),
it holds that

I = 1 vol(M) S d-1 r 0 f (x 0 )t d+3 ⟨ι * θ, e k ⟩⟨ι * θ, e l ⟩⟨ι * θ, e m ⟩⟨ι * θ, e n ⟩ + O(t d+4 )dtdθ = f (x 0 ) vol(M)(d + 4) r d+4 S d-1
⟨ι * θ, e k ⟩⟨ι * θ, e l ⟩⟨ι * θ, e m ⟩⟨ι * θ, e n ⟩dθ + O(r d+5 ); as claimed.

We can now state our main theorem. Note that we consider functions defined on the ambient space, R p as opposed to only on M, since in manifold learning M is unknown and the operator will thus be applied to function on R p . Theorem 1. Consider g ∈ C 2 , g : R p → R. Denote by Q 0 the Hessian of g at x 0 . To simplify notation, set X 0 = x 0 . Take a sample {ι(X j )} N j=1 from the uniform distribution on M which is embedded in R p . Then, under Assumptions 1 and 2, defining

∆ OT g(X 0 ) := N j=0 W ε 0,j g(X 0 ) -g(X j ) ,
with W ε the approximate solution of the quadratically regularised OT problem, it holds that,

-2K -1 ε,N ∆ OT g(X 0 ) L 2 --→ dL 0 0, J⊤ p,p-d N(x 0 ) 2 + 1 2 tr Q 0 I d×d 0 0 0 , provided that K ε,N → 0 and N/ε → 0 when N → ∞.
Proof. As g is twice differentiable, we can write

g(X j ) = g(x 0 ) + L 0 (x 0 -X j ) + 1 2 (X j -x 0 ) ⊤ Q 0 (X j -x 0 ) + o ∥x 0 -X j ∥ 2 ,
where L 0 is the gradient of g at x 0 and Q 0 is the Hessian of g evaluated at x 0 . Plugging this result in the definition of (∆ OT g)(x 0 ), we derive

(∆ OT g)(x 0 ) = N j=1 W ε 0,j -L 0 (X 0 -X j ) - 1 2 (X j -X 0 ) ⊤ Q 0 (X j -X 0 ) + o ∥X 0 -X j ∥ 2 .
We will split this sum into three terms and control each one separately. We will first consider the expectation and then the variance.

Step 1: Expectation.

Let us start with the second term

- 1 2 N j=0 W ε 0,j (X j -X 0 ) ⊤ Q 0 (X j -X 0 ) =: B
The quantity in the above display is a scalar so that it is equal to its trace. Further, using the linearity and the cyclical property of the trace, it holds that

B = - 1 2 tr   Q 0 N j=0 W ε 0,j (X j -X 0 )(X j -X 0 ) ⊤   = - 1 2 (N + 1) tr   Q 0 N j=0 (K ε,N -c 0,j ) + ε (X j -X 0 )(X j -X 0 ) ⊤   .
Using Lemma 1, it holds that

E N j=1 (K ε,N -c 0,j ) + ε (X j -X 0 )(X j -X 0 ) ⊤ = N ε |S d-1 | d(d + 2) vol(M) K d+2 2 ε,N K ε,N I d×d 0 0 0 + O(K ε,N ) .
Let us now address the first term, i.e.,

-L 0 N j=0 W ε i,j (X 0 -X j ).
The second part of Lemma B.5 in [START_REF] Wu | Think globally, fit locally under the manifold setup: Asymptotic analysis of locally linear embedding[END_REF] reads, in our case,

E [(X -ι(x 0 ))f (X; r)1{∥X -ι(x 0 )∥ ≤ r}] = |S d-1 | (d + 2) vol(M) J ⊤ p,d ι * ∇f (x 0 ; r) d , f (x 0 ; r) J⊤ p,p-d N(x 0 ) 2 r d+2 + O(r d+4 ).
It follows that

E N j=1 W ε 0,j (X 0 -X j ) = N (N + 1) ε |S d-1 | (d + 2) vol(M) 0, K ε,N J⊤ p,p-d N(x 0 ) 2 K d+2 2 ε,N + O N (N + 1) ε K (d+4)/2 ε,N
.

In view of the developments above, the expectation of the Taylor residual is negligible.

Step 2: Variance.

Let us deal with

-L 0 N j=1 W ε 0,j (x 0 -X j ).
We have that

Var N j=1 W ε 0,j (x 0 -X j ) = N (N + 1) 2 ε 2 Var((K ε,N -c 0,j ) + (X -x 0 )) Further, Var[(K ε,N -c(X, x 0 )) + (X -x i )] = E[(K ε,N -c(X, x 0 )) 2 + (X -x 0 )(X -x 0 ) ⊤ ] -E[(K ε,N -c(X, x 0 )) + (X -x i )]E[(K ε,N -c(X, x 0 )) + (X -x 0 ) ⊤ ] = |S d-1 | d(d + 2) vol(M) κ d+2 2 d K 2 ε,N εN -2 I d×d 0 0 0 + O(K ε,N ) - |S d-1 | 2 (d + 2) 2 vol 2 (M) K 2 ε,N vv ⊤ κ d+2 d ε 2 N -4 + O K (d+4)/2 ε,N εN -2 K ε,N ,
where v are vectors that depend on the curvature as above. Because of the rescaling by K -1 ε,N , we finally get

Var   K -1 ε,N N j=1 W ε i,j (x i -X j )   = O K -2 ε,N N 2 N ε 2 K 2 ε,N N -2 ε = O N ε .
Let us now turn to the covariance matrix of

V := vec   j W ε i,j (x i -X j )(x i -X j ) ⊤  
which, using Equations (1.3.14), (1.3.16) and (1.3.31) in [START_REF] Kollo | Advanced multivariate statistics with matrices[END_REF], is equal to

N E[(W ε j ) 2 (x i -X j ) ⊗ (x i -X j ) ⊤ ⊗ (x i -X j ) ⊗ (x i -X j ) ⊤ ] -N EV E ⊤ V.
Relying on Lemma 2, we get that the leading order of the variance of

K -1 ε,N V is O K -2 ε,N N 2 N ε 2 K 2 ε,N K d+4 2 ε,N = O N K ε,N ε .
The claim follows.

Infinitesimal generator limit and spectral convergence

A relatively general analysis of convergence of graph Laplacians was carried out by Ting, Huang, and Jordan (2010), wherein consistency results are established for a general class of constructions leveraging connections to diffusion processes. We remark that when M is endowed with a uniform measure, a constant approximation of the potential is valid and so the operator resulting from quadratically regularized optimal transport falls under their framework (Ting, Huang, and Jordan, 2010, Theorem 3).

The assumptions are compatible with the ones that we make here, namely that M is a smooth, compact manifold, and the authors consider a general kernel of the form

K N (x, y) = w (N ) x (y)K 0 ∥y-x∥ h N r (N ) x (y)
.

In our setting where iid samples are drawn uniformly on M, we invoke a constant potential approximation u ∼ ε 2 2+d N -4 2+d , we have (up to a multiplicative constant)

K N (x, y) = ε 2 2+d N -4 d+2 -∥y -x∥ 2 + = 1 - ∥y -x∥ ε 1 2+d N -2 d+2 2 + = φ ∥y -x∥ h (N )
Where the choice of kernel is the Epanechnikov kernel φ(r) = (1 -r 2 ) + . The condition under which their theorem holds is that N h m+2 / log N → ∞. In our case, this simplifies to ε/(N log N ) → ∞, and this is compatible with the range of scalings N 1-2/d ≪ ε ≪ N 2 from our previous analysis.

It is further possible understand how the eigenvalues and eigenvectors of the discrete operator relate to the continuous one, relying on the recent results by [START_REF] Trillos | Error estimates for spectral convergence of the graph Laplacian on random geometric graphs toward the Laplace-Beltrami operator[END_REF]. Their results apply in the setting that we consider for an intrinsic dimension d ≥ 2. Upon choosing

ε = N 3d+2 2d(d+2) (log N ) p d (d+2

Equispaced points on the circle

We finally consider an example for which the computations can be explicitly carried out: the case of equidistant points on the circle.

First-order optimal potentials

We finish this section about the rates in the discrete case in a one-dimensional deterministic example.

Consider N points that are equispaced on the circle each with mass 1/N . Set D j to be the j-th squared Euclidean distance in the sorted list of all distances from one point to the others. We thus have

D j = 2 sin πj N 2 = 4π 2 j 2 N 2 + O j 4 N 4
We aim at solving

N j=1 (y -D j ) + = ε N .
We get, for y small, that there exists k such that

4π 2 N 2 (k) 2 ≤ y ≤ 4π 2 N 2 (k + 1) 2 , that (2k + 1)y -2 4π 2 N 2 k j=1 j 2 = ε N .
Thus, there exists α such that

(2k + 1)(k + α) 2 -2 k j=1 j 2 = εN 4π 2 (2k + 1)(k + α) 2 - 1 3 k(k + 1)(2k + 1) = εN 4π 2 .
It follows that by matching the largest order for k

4k 3 3 ≈ εN 4π 2
and thus

y ≈ 4π 2 N 2 3εN 16π 2 2/3
Alternatively, to asses the quality of the approximation above, consider 4 sin 2 πk N ≤ y < 4 sin 2 π(k + 1) N , to derive

(2k + 1)4 sin 2 π(k + α) N -2 k j=1 4 sin 2 πj N = ε N (10)
It further holds that

k j=1 4 sin 2 πj N = 1 + 2k - sin π(2k+1) N sin(π/N ) = 1 + 2k - (2k + 1)π/N -(2k+1) 3 π 3 6N 3 + O (2k + 1) 5 /N 5 π/N + O(1/N 3 ) = 1 + 2k - (2k + 1) -(2k+1) 3 π 2 6N 2 + O (2k + 1) 5 /N 4 1 + O(1/N 2 ) = (2k + 1) 3 π 2 6N 2 + O 2k + 1 N 2 + O (2k + 1) 5 N 4 = 4k 3 2 3N 2 + O k 2 N 2 + O (2k + 1) 5 N 4 .
Plugging this result in (10), one gets

(2k + 1) 4π 2 (k + α) 2 + O k 4 N 4 -2 4k 3 π 2 3 + O k 2 + O (2k + 1) 5 N 2 = εN which gives 4k 3 3 + O k 2 + O k 5 N 4 + O (2k + 1) 5 N 2 = εN 4π 2 .
This matches with the other approximation.

Limiting operator

Let us place ourselves in the same setting as Section 5.1 again. We have seen that the optimal potential must behave as κε 2/3 N -4/3 at the first order.

Theorem 2. Consider g ∈ C 2 , g : R 2 → R. Consider again a set {x i } N i=1 of N equispaced points on the unit circle. For simplicity, choose i ≤ N such that x i = (0, 1). Denote by Q i , the Hessian of g at x i . Then, defining

(∆ OT g)(x i ) := N j=1 W i,j g(x i ) -g(x j ) , with W as above, there exist constants C 1 , C 2 , such that C 1 ε -2 3 N 4 3 (∆ OT g)(x i ) → C 2 ∂g(z, y) ∂y (z,y)=(0,1) + tr Q i 0 0 0 1 = C 2 ∂g(z, y) ∂y (z,y)=(0,1) + ∂ 2 g(z, y) ∂z 2 (z,y)=(0,1) , for lim N →∞ εN = ∞.
Proof. As g is twice differentiable, we can write

g(x j ) = g(x i ) + L i (x i -x j ) + 1 2 (x j -x i ) ⊤ Q i (x j -x i ) + o ∥x i -x j ∥ 2 ,
where L i is linear and Q i is the Hessian of g evaluated at x i . Plugging this result in the definition of (∆ OT g)(x i ), we derive

(∆ OT g)(x i ) = N j=1 W i,j -L i (x i -x j ) - 1 2 (x j -x i ) ⊤ Q i (x j -x i ) + o ∥x i -x j ∥ 2 .
We will split this sum into three terms and control each one separately. Let us start with the the second term

- 1 2 N j=1 W i,j (x j -x i ) ⊤ Q i (x j -x i ) =: B.
The quantity in the above display is a scalar so that it is equal to its trace. Further, using the linearity and the cyclical property of the trace, it holds that

B = - 1 2 tr   Q i N j=1 W i,j (x j -x i )(x j -x i ) ⊤   .
Relabelling the points from closest to i to furthest, computing explicitly (

x j -x i )(x j -x i ) ⊤ , we get B = -N tr   Q i ⌊N/2⌋ j=1 (κε 2/3 N -4/3 -4π 2 j 2 N 2 + O(j 4 /N 4 )) + ε M j   ,
where

M j := sin 2 ( 2πj N ) -2 sin 2 ( πj N ) sin( 2πj N ) -2 sin 2 ( πj N ) sin( 2πj N ) 4 sin 4 ( πj N ) sin 2 ( 2πj N )
.

A first order development gives

M j ≈ 4π 2 N 2 j 2 -π 2N j 3 -π 2N j 3 π 2 j 4 N 2 . Thus, B ≈ - 4π 2 N tr Q i E 1 0 0 0 , where E := 4π 2 εN 2 κ 4π 2 (εN ) 2/3 1 3 κ 4π 2 3/2 εN + O (εN ) 2/3 - 1 5 κ 4π 2 5/2 (εN ) 5/3 + O (εN ) 4/3 . So that B ≈ ε 2 3 N -4 3 C 1 tr Q i 1 0 0 0 + O ε 1/3 N -2/3 N ,
for a constant C 1 . Let us now consider the Taylor residual. It holds that

N j=1 W i,j o ∥x i -x j ∥ 2 ≈ N N j=1 (κε 2/3 N -4/3 -∥x i -x j ∥ 2 ) + ε o ∥x i -x j ∥ 2 ≈ 2N ⌊N/2⌋ j=1 (κε 2/3 N -4/3 -4π 2 j 2 N 2 + O(j 4 /N 4 )) + ε o j 2 N 2 .
Comparing with the developments above, we see that the sum is a weighted sum of o(j 2 /N 2 ) whereas, for B it was the sum of terms behaving like j 2 /N 2 with the same weights. It follows that

ε -2 3 N 4 3 N j=1 W i,j o ∥x i -x j ∥ 2 = o(1).
Let us finally address the first term, i.e.,

-L i N j=1 W i,j (x i -x j ),
Owing to the symmetry of the problem, this is constant times the normal vector at the point. It is thus non zero in the y direction. Remark that the vector x j -x i has a component in the y direction equal to 2 sin 2 (πj/N ) and it thus has the same leading order as M j in terms of ε, N . The claim follows.

Simulations

In this section we exhibit the size of the optimal potentials obtained from the semismoothed Newton algorithm proposed in [START_REF] Lorenz | Quadratically regularized optimal transport[END_REF] and adapted to our setting in [START_REF] Matsumoto | Beyond knn: Adaptive, sparse neighborhood graphs via optimal transport[END_REF].

d-Sphere

We now exhibit the behaviour of the optimal potentials for N random points on the d-Sphere and various parameters ε. In d = 1, 2, 3, N = 10 3 points were sampled uniformly by sampling from S d first a standard Gaussian and normalizing. We numerically solved the corresponding discrete optimal transport problem with ε in the range [10 -3 , 10 5 ] and plotted log(u) against log(ε). For ε sufficiently large, we estimated the exponent α for the relationship u ∼ ε α . Our empirical findings agree with the exponent 2 2+d .

Figure 1: Scaling of dual potential u on the d-Sphere.

Torus

Next, we investigate the behaviour of the operator ∆ OT in the discrete setting where points are sampled from the uniform distribution on the 2-dimensional torus with major and minor radii R = 1, r = 1/2. We fix a point x 0 = (0, 1/2, 0) at which the tangent space T x0 M is spanned by e 1 , e 3 . We then consider a function f (x, y, z) = 3x 2 +5y 2 +7z 2 . For N = 100, 250, 500, 1000, 2500 points sampled from the torus, we calculated the (N + 1) × (N + 1) coupling π by solving (1), normalized following (5), and then computed the quantity K -1 ε,N (∆ OT f )(x 0 ). Motivated by the asymptotic scalings we derived, we tried setting ε ∝ N α for varying exponents: α = 2 which should correspond to a fixed regularization level in the continuous case (and we do not expect convergence to the Laplacian in this case), and α = 1.05, 1.125, 1.25, 1.5, 1.75 which all fall within the regime where Theorem 1 applies. We show in Figure 2 the values of K -1 ε,N (∆ OT f )(x 0 ) over 10 repeats at each value of N .

We see that when α = 2, the quantity K -1 ε,N (∆ OT f )(x 0 ) stabilizes around a fixed value as N increases. This agrees with our understanding that ε ∝ N 2 in the discrete setting corresponds to the continuous case of empirical distributions with a fixed value of ε. On the other hand, when 1 < α < 2 we observe a pattern of values appears to converge around a different value. Importantly, for various 1 < α < 2, these values are similar -this supports the scaling relation of Theorem 1 and suggests that the quantity is converging to the value (up to a constant independent of ε, N ) of the Laplace-Beltrami operator at x 0 . ε,N (∆ OT f )(x 0 ) for varying sample sizes N , and ε ∝ N α , for various choices α = 2, 1. 75, 1.5, 1.25, 1.125, 1.05.

) 2 ,

 2 with p d = 3/4 if d = 2 and p d = 1/d if d ≥ 3, the rate of convergence of the eigenvalues and eigenfunctions 4 is O (log N ) p d N 1/d almost surely (García Trillos et al., 2020, Theorems 1 and 5).

Figure 2 :

 2 Figure 2: Estimate of Laplacian K -1ε,N (∆ OT f )(x 0 ) for varying sample sizes N , and ε ∝ N α , for various choices α = 2, 1.75, 1.5, 1.25, 1.125, 1.05.

Recall that a manifold is closed if it is compact and without boundary.

Note that the measures in their paper are not probability measures, which explains the slight difference.

We refer to the paper for an explicit description of how the eigenvector is interpolated to compute the norm between that interpolation and the eigenfunction on the manifold.

* G. Mordant gratefully acknowledges the funding by the DFG for CRC1456. † S. Zhang gratefully acknowledges funding from the Australian Government RTP Program.


Côme Le Breton

Théodore Papadopoulo

Maureen Clerc

Investigating the origin of phase delays in the alpha band

Results

B] There is a gradual phase shift in the alpha band, in the anteroposterior direction, that is visible on single recordings, and which is subject specific. C] Phase difference peaks are stable across recordings. Non zero phase difference peaks cannot be observed in control conditions such as during calculus blocks, in the beta band or along the coronal direction. These phase peaks would, at 10Hz, amount to propagation speed of 7.3 m.s -1 . This is incompatible (too fast) with transmission delays in the cortex or through white matter. D] The geometry of the model with 2 sources explains the phase difference observed in the real data. With a phase shift of 0.93 π , one can observe similar phase patterns as those of the real data. E] The PLI shows spurious functionnal connectivity between electrodes along the anteroposterior axis that are likely due to only two functionally connected sources. In scalp EEG, volume conduction considerably complicates the problem, and allows to consider a fourth "mechanism", in which only two sources can generate phase delays at many electrodes.

Phase differences in the α band appear along the anteroposterior axis during the ECR condition, and this pattern is stable across recordings taken years apart.

A simplistic model of two sources is enough to model phase differences incompatible with transmission delays.

Care must therefore be brought to the apriori misconception that PLI index is not sensitive to volume conduction.

Code

  Côme Le Breton, Théodore Papadopoulo, Maureen Clerc Athena research team, Inria Sophia AntipolisTraveling (or propagating) wave are oscillating electrical voltages observed sequentially or simultaneously at different brain locations, relating them functionally. They are observed at various scales, all over the brain, in various frequency bands and following varied stimuli (or in their absence), making them an ubiquitous phenomenon.Ermentrout and Kleinfeld [2001] proposed three possible mechanisms: 1) A unique oscillator activating cortical areas with different time delays (resulting in a spurious wave), 2) A chain of successively activated oscillators, 3) Coupled local oscillators generating wave propagation due to phase delay.

MethodologyA]

  EEG data acquisition split over several recordings (up to 4 per session) and over several sessions (up to 8 per subject). Each recording lasts 275 seconds and alternates between 20 second blocks of Eyes Open Fixation (EOF) and Eyes Closed Relaxation (ECR). EEG channels of interest are re-referenced to M2. Instantaneous phase is extracted with the Morlet wavelet transform in the alpha band (8-12Hz) and split according to condition. B] Histograms of phase differences are then constructed, for each considered pair and each condition. C] To show the stability of the distributions across recordings, the parameters of the fitted von Mises distributions are extracted from each block, and the locations are plotted on a circular histogram. D] To explain the phase difference observed on real data, a simulation based on a model of 2 sources is devised. One source is located in the frontal area and another in the occipital area. E] The Phase Lag Index is computed on short windows and averaged over 5 blocks of ECR.