Mircea Moscu
email: mircea.moscu@oca.eu

Ricardo A Borsoi
email: raborsoi@gmail.com.

Cédric Richard
email: cedric.richard@unice.fr

José-Carlos M Bermudez
email: j.bermudez@ieee.org.

C M Bermudez

Graph topology inference with derivative-reproducing property in RKHS: algorithm and convergence analysis

Keywords: convergence analysis, kernel least mean squares, nonlinear topology inference, partial derivative sparsity

In many areas such as computational biology, finance or social sciences, knowledge of an underlying graph explaining the interactions between agents is of paramount importance but still challenging. Considering that these interactions may be based on nonlinear relationships adds further complexity to the topology inference problem. Among the latest methods that respond to this need is a topology inference one proposed by the authors, which estimates a possibly directed adjacency matrix in an online manner. Contrasting with previous approaches based on linear models, the considered model is able to explain nonlinear interactions between the agents in a network. The novelty in the considered method is the use of a derivative-reproducing property to enforce network sparsity, while reproducing kernels are used to model the nonlinear interactions. The aim of this paper is to present a thorough convergence analysis of this method. The analysis is proven to be sane both in the mean and mean square sense. In addition, stability conditions are devised to ensure the convergence of the analyzed method.

I. INTRODUCTION

Knowledge of the network topology in applications such as gene regulation systems [START_REF] Frazer | A second generation human haplotype map of over 3.1 million SNPs[END_REF], socio-economic interactions [START_REF] Heiberger | Predicting economic growth with stock networks[END_REF], or brain activity [START_REF] Kramer | Emergent network topology at seizure onset in humans[END_REF] is of utmost importance. The main reason is that processing these data via, e.g., filtering [START_REF] Sandryhaila | Discrete signal processing on graphs: Graph filters[END_REF], [START_REF] Nassif | A graph diffusion LMS strategy for adaptive graph signal processing[END_REF], [START_REF] Chen | Multitask diffusion adaptation over networks with common latent representations[END_REF], [START_REF] Coutino | Advances in distributed graph filtering[END_REF], [START_REF] Koppel | Decentralized online learning with kernels[END_REF] or spectral analysis [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF], [START_REF] Tremblay | Accelerated spectral clustering using graph filtering of random signals[END_REF], requires the network structure [START_REF] Djuric | Cooperative and Graph Signal Processing: Principles and Applications[END_REF]. Since most graph signal processing algorithms assume that the graph topology is known beforehand, significant progress has been made recently in the estimation of the graph topology from available data, be them links between regions of the brain, genes in a network, or sectors of a market economy. Moreover, the presence of nonlinear interactions in real-world applications led to the development of more general algorithms. One such proposition is the graph topology inference algorithm recently proposed by the authors [START_REF] Moscu | Online kernel-based graph topology identification with partial-derivative-imposed sparsity[END_REF], which is based on derivative-reproducing property to infer nonlinear interactions and promote sparsity. The current paper builds upon this preliminary work by proposing an in-depth and complex stability and performance analysis of this algorithm. The ability of reproducing kernels to model nonlinear relationships between nodal signals is employed in the graph inference process, alongside with kernel dictionaries to mitigate the increasing number of kernel functions due to the online setting, while distributing the computational burden over the agents in the network. Many real world examples, such as social graphs [START_REF] Speriosu | Twitter polarity classification with label propagation over lexical links and the follower graph[END_REF], show considerable edge sparsity, which needs a sparsity-inducing framework based on derivative-reproducing kernels. The present paper briefly recalls the considered algorithm and follows up with a complete convergence analysis, both in the mean and mean-square sense.

A. Prior works

Early advances in topology inference are put forward in [START_REF] Dempster | Covariance selection[END_REF], where log-likelihood-based method of covariance estimation is introduced. Similarly, in [START_REF] Friedman | Sparse inverse covariance estimation with the graphical Lasso[END_REF] the graphical Lasso is employed in order to estimate the precision matrix from the available data. In [START_REF] Segarra | Network topology inference from spectral templates[END_REF], the authors advocate that connectivity can be recovered from estimated spectral templates, while [START_REF] Sardellitti | Graph topology inference based on sparsifying transform learning[END_REF] adapts the method for band-limited signals, i.e., signals whose Graph Fourier Transform is sparse. The authors of [START_REF] Shafipour | Network topology inference from non-stationary graph signals[END_REF] introduce a method designed for topology inference for the case when the measured signals are non-stationary. An online adaptive algorithm is developed in [START_REF] Moscu | Learning causal networks topology from streaming graph signals[END_REF], where the authors assume a linear model governing the interactions. Linearity of interactions and signal stationarity are assumed in [START_REF] Shafipour | Online topology inference from streaming stationary graph signals[END_REF] to devise an ADMM algorithm.

The modeling of nonlinear interactions has proved to be a relevant research topic. Some existing methods of higherorder link inference yield satisfactory results when applied in optimizing power grids [START_REF] Coutino | Selfdriven graph Volterra models for higher-order link prediction[END_REF] or in predicting social behavior [START_REF] Yang | Learning connectivity and higher-order interactions in radial distribution grids[END_REF]. In the presence of nonlinear phenomena, works such as [START_REF] Harring | A comparison of methods for estimating quadratic effects in nonlinear structural equation models[END_REF], [START_REF] Finch | Modeling nonlinear structural equation models: A comparison of the two-stage generalized additive models and the finite mixture structural equation model[END_REF] focus on polynomial structural equation models, while the authors of [START_REF] Lim | Operatorvalued kernel-based vector autoregressive models for network inference[END_REF] use their nonlinear counterparts. They, however, have some limitations, such as assuming prior knowledge of certain connections or the form of the nonlinear basis functions. Reproducing kernels have seen extensive use in graph topology inference problems. Such a solution was developed in [START_REF] Lippert | A kernel method for unsupervised structured network inference[END_REF], where an unsupervised kernel-based method is implemented. One particularity of the algorithm is that it requires, as a parameter, the number of sought edges. It also offers the possibility of statistical significance testing when setting this parameter. Another work is [START_REF] Shen | Topology inference of directed graphs using nonlinear structural vector autoregressive models[END_REF] where kernels, chosen to best fit the available data, model nonlinear relationships between nodes based on measurements at successive time instants. The multi-kernel approach in [START_REF] Zhang | Going beyond linear dependencies to unveil connectivity of meshed grids[END_REF] uses partial correlations to encode graph topology and p -norm regression to enhance the performance. In [START_REF] Giannakis | Topology Identification and Learning over Graphs: Accounting for Nonlinearities and Dynamics[END_REF], a thorough analysis of the kernel-based topology inference problem is given. This work focuses on capturing nonlinear and dynamic links. A review of state-of-the art methods in topology inference is in [START_REF] Dong | Learning graphs from data: A signal representation perspective[END_REF].

None of these works, however, bases its rationale on kernel derivatives. Moreover, no analysis of this family of methods exists in the literature. For these reasons, we put forth the current paper, which first recalls the algorithm and then proposes an in-depth convergence analysis.

B. Notations and definitions

Normal font letters denote scalars, while boldface lowercase and uppercase letters stand for column vectors and matrices, respectively. Uppercase calligraphic letters denote sets. Their cardinality is denoted by | • |. The (u, v)-th entry of a matrix X is denoted by [X] u,v or by x uv . Notation α = col{{α p } p } denotes the column vector obtained by stacking all entries α p for all available indexes p. Finally, E{•} is the expectation operator.

A graph G consists of a set N of (N + 1) nodes, and a set E of edges such that if nodes m and n are linked, then we have (m, n) ∈ E. For undirected graphs, these node pairs are unordered. At node level, we collect a real-valued signal and we organize it as a column vector [y 1 (i), . . . , y N +1 (i)] , where y n (i) is the signal sample at node n and time instant i. The adjacency matrix A [START_REF] Biggs | Algebraic Graph Theory[END_REF], [START_REF] Sandryhaila | Big data analysis with signal processing on graphs: Representation and processing of massive data sets with irregular structure[END_REF], is defined as an (N + 1) × (N + 1) matrix whose entries a nm are zero if (m, n) / ∈ E and set to one otherwise. Throughout the paper we shall consider that there is no self-loop, i.e., a nn = 0, for all n.

C. Paper outline

The paper is organized as follows. Sections II and III recall the problem previously considered in [START_REF] Moscu | Online kernel-based graph topology identification with partial-derivative-imposed sparsity[END_REF] and the online algorithm, respectively. Section IV then proceeds with the main contribution of the current paper, consisting of the analysis of the algorithm, both in mean and mean square sense. In order to validate the pertinence of the analysis, Section V showcases a set of experimental simulations. Finally, a set of concluding remarks is given in Section VI.

II. LOCAL PROBLEM FORMULATION

Consider an (N + 1)-node graph with adjacency matrix A that models a system such as a brain network or a power grid. In this context, the electrical activity of different brain regions [START_REF] Rubinov | Complex network measures of brain connectivity: uses and interpretations[END_REF], [START_REF] Shen | Nonlinear structural vector autoregressive models with application to directed brain networks[END_REF], or the voltage angle per bus [START_REF] Zhang | Going beyond linear dependencies to unveil connectivity of meshed grids[END_REF], numbered from 1 to (N + 1), can be measured sequentially at different time instants i ∈ N, leading to a graph signal given by [y 1 (i), . . . , y N +1 (i)] ∈ R N +1 that changes over time. Signal y n (i) at each node n = 1, . . . , N + 1 of the graph is nonlinearly coupled to the signals at all nodes in its neighborhood, according to the topology described in A. Moreover, we assume that nonlinear relationships between signals at different nodes occur as have been reported in many applications; see, e.g., the case of brain connectivity [START_REF] Freeman | EEG analysis gives model of neuronal templatematching mechanism for sensory search with olfactory bulb[END_REF], [START_REF] De Zwart | Hemodynamic nonlinearities affect BOLD fMRI response timing and amplitude[END_REF].

Recent methods have considered additive nonlinear models [START_REF] Shen | Kernel-based structural equation models for topology identification of directed networks[END_REF], [START_REF] Moscu | Online graph topology inference with kernels for brain connectivity estimation[END_REF], which represent the local measurements at each node n as:

y n (i) = m∈N \{n} f nm (y m (i)) + n (i) , (1)
where

y m (i) = y m (i), . . . , y m (i -L m + 1) , for L m ≥ 1.
Parameter L m endows the algorithm with memory-like capacities. This characteristic is desirable in a number of applications such as brain topology inference where there is a 10-20 ms delay in signal propagation between nodes [START_REF] Petkoski | Transmission time delays organize the brain network synchronization[END_REF]. Functions f nm : R Lm → R represent the interactions between the different nodes in the network, and n (i) denotes innovation noise. By relating how each node m ∈ N \ {n} affects node n, functions f nm encode the connectivity in A directly as a nm = 0 if, and only if, f nm ≡ 0, where a nm is the (n, m)th entry of A. For ease of notation, we shall assume that n ≡ (N + 1), i.e., we identify n as the (N + 1)-th node of the graph, which allows us to denote N \ {n} = {1, . . . , N }.

The topology inference problem then consists of finding the set of functions f nm which best represent the available graph signal measurements y n (i) based on the previous model. This can be performed using both batch-based (see, e.g., [START_REF] Shen | Kernel-based structural equation models for topology identification of directed networks[END_REF]) and online strategies [START_REF] Moscu | Online graph topology inference with kernels for brain connectivity estimation[END_REF], [START_REF] Moscu | Convergence analysis of the graph-topology-inference kernel LMS algorithm[END_REF], [START_REF] Shen | Online identification of directional graph topologies capturing dynamic and nonlinear dependencies[END_REF]. Models such as (1) do not consider general nonlinear interactions between multiple nodes, as they assume an additive model for y n (i) [START_REF] Buja | Linear smoothers and additive models[END_REF]. This can be rather limiting since the nonlinearity introduced in the model through the functions f nm is local, i.e., it only acts on the signals at each node individually. To overcome this issue, we propose to consider the following general nonlinear model:

y n (i) = f n (y(i)) + n (i) , (2)
where y(i) = [y 1 (i) , . . . , y N (i) ] and function f n describes the interaction between the signals at nodes N \ {n} and the signal y n (i) at node n. Compared to (1), model (2) captures more complex relationships between nodes without relying on generalized linear model, rendering it more general. The method we propose relates f n to the underlying graph topology A by quantifying how a node m affects node n through the corresponding partial derivative, i.e.: node m does not affect node n (i.e., a mn = 0) ⇐⇒ ∂f n (y)

∂y m = 0 , (3)
under the assumption that f n is continuously differentiable, where the expectation is taken w.r.t. y.

Estimating the graph topology can then be performed by learning function f n that fits the data as per model (2), while additionally considering penalties to promote a sparsely connected graph, as in many real-world applications [START_REF] Danisch | Listing k-cliques in sparse real-world graphs[END_REF]. Afterwards, the topology can be recovered via relation [START_REF] Kramer | Emergent network topology at seizure onset in humans[END_REF]. By constraining f n to belong to a Reproducing Kernel Hilbert Space (RKHS) H associated with a positive definite repro-ducing kernel κ(•, •), this problem can be formulated in batch form as [START_REF] Moscu | Online kernel-based graph topology identification with partial-derivative-imposed sparsity[END_REF]:

min fn∈H 1 2i i =1 y n ( ) -f n (y( )) 2 + η   N m=1 1 i i p=1 Lm q=1 ∂f n (y(p)) ∂y m,q 2 + ψ f n H   , (4)
where y m,q represents the q-th entry of y m . The first term in (4) is a data fitting term that averages the reconstruction error over all available samples. The second one is a convex groupsparsity inducing penalty [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF], [START_REF] Rosasco | Nonparametric sparsity and regularization[END_REF], which promotes sparsity in the estimated derivatives averaged over all observations y(p) up to time instant i, and consequently on adjacency matrix A. The last term, with ψ : R → [0, ∞), is a monotonically increasing function of f n H with small magnitude in order to guarantee the uniqueness of the solution of (4). Parameter η makes a trade off between data the fitting term and the regularizers. Problem (4) allows us to introduce sparsity in f n without the inherent limitations of an additive model. However, the sparsity-promoting penalty term prevents the application of previously established Representation Theorems that guarantee the existence of a finite-dimensional representation of the solution. Assuming that kernel κ(•, •) is twice differentiable, the following holds [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF]:

H ∂f n (y) ∂y m,q = f n , κ ∂m,q (•, y) H , (5)
where:

κ ∂m,q (•, y(q)) = ∂κ(•, a) ∂a m,q a=y(q) .

Thus, for sufficiently smooth kernels, the derivative of functions in H belongs to the same H, in turn meaning that they can be evaluated in the form of simple inner products. The solution of (4) can now be written in a finite dimensional form, similarly to the proposition in [START_REF] Rosasco | Nonparametric sparsity and regularization[END_REF]. By generalizing the Representer Theorem in the aforementioned paper, we obtain a Representer Theorem for our approach. This result is stated in the following theorem:

Theorem 1. Let H be a RKHS whose associated reproducing kernel κ(•, •) is at least twice differentiable. Then, the solution of the optimization problem (4) can be written as:

f * n = i p=1 α p κ(•, y(p)) + N m=1 i q=1 Lm =1 β m,q, κ ∂ m, (•, y(q)) . (7)
Proof. The proof of this theorem is presented in Appendix A.

Expression [START_REF] Coutino | Advances in distributed graph filtering[END_REF] can then be substituted in problem (4) in order to obtain a finite-dimensional optimization problem.

In the sequel, for the sake of simplicity in the notation and ease of comprehension, we shall consider the case L m = 1 for all m, i.e., y m (i) = y m (i) and κ ∂m,q (•, y) = κ ∂m (•, y).

III. AN ONLINE ALGORITHM

The number of coefficients α p and β m,q in solution (7) can become prohibitive as i increases with each new measurement and associated kernel function. Kernel dictionaries are a solution to this issue. These dictionaries can be defined either a priori [START_REF] Chen | Convergence analysis of kernel LMS algorithm with pre-tuned dictionary[END_REF] or can admit a new candidate kernel function only if it satisfies a certain sparsification rule [START_REF] Gao | Kernel lms algorithm with forward-backward splitting for dictionary learning[END_REF], such as the approximate linear dependence criterion [START_REF] Engel | The kernel recursive least-squares algorithm[END_REF], [START_REF] Bueno | Gram-Schmidt-based sparsification for kernel dictionary[END_REF], the Nyström method [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF], dictionary-based random Fourier features [START_REF] Bouboulis | Efficient KLMS and KRLS algorithms: A random Fourier feature perspective[END_REF], orthogonal projection [START_REF] Takizawa | Adaptive nonlinear estimation based on parallel projection along affine subspaces in reproducing kernel hilbert space[END_REF] or the novelty criterion [START_REF] Platt | A resource-allocating network for function interpolation[END_REF].

We consider a dictionary-based framework, where each node n in the network is equipped with a dictionary of kernel functions (and their derivatives) of the following form:

D n = κ(•, y(ω j )) : ω j ∈ I i n , (8)
where I i n ⊂ {1, . . . , i -1} represents the set of time indices of elements selected for the dictionary, before instant i. To build D n iteratively, we consider the coherence criterion as the chosen method for dictionary sparsification [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. This entails the fact that, after a sufficient number of samples i, only a small number |I i n | i of coefficients will be needed. A candidate kernel function κ(•, y(i)) is added to D n if the following sparsification condition holds [START_REF] Richard | Online prediction of time series data with kernels[END_REF]:

max ωj ∈I i n |κ(y(i), y(ω j ))| ≤ ξ n , (9)
where ξ n ∈ [0, 1) determines the level of sparsity and the coherence of the dictionary. The number of entries in the dictionary satisfies |I i n | < ∞ when i → ∞ [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. Note that the dictionary D n can also be set a priori using, e.g., an uniform grid sampling [START_REF] Chen | Convergence analysis of kernel LMS algorithm with pre-tuned dictionary[END_REF]. In this case, with a slight abuse of notation, we represent the dictionary using a set of time indices I i n ⊂ Z \ N, attributing to the dictionary elements negative time indexes. We rewrite [START_REF] Coutino | Advances in distributed graph filtering[END_REF] as:

f n,Dn = |Dn| p=1 α p κ(•, y(ω p )) + N m=1 |Dn| q=1 β m,q κ ∂m (•, y(ω q )) . ( 10
)
For notation compactness, we shall write α = col{{α p } p }, and

β = col{{β m } m } with β m = col{{β m,q } q }.
The dictionary-based representation f n,Dn in (10) can now be plugged directly into the batch problem (4) to obtain a solution depending on a small amount of coefficients, for all i. However, this would still entail a growth in complexity with time i, since the kernel functions have to be evaluated at all available data points. To address this issue and devise a computationally scalable algorithm, we consider an alternative cost function to [START_REF] Sandryhaila | Discrete signal processing on graphs: Graph filters[END_REF]. Specifically, we substitute the averaging operations in the data fitting term and in the derivative penalty term in [START_REF] Sandryhaila | Discrete signal processing on graphs: Graph filters[END_REF], which use all data points up to time instant i, by the corresponding expected values. Setting ψ ≡ 0 without loss of generality, this leads to:

min f n,Dn 1 2 E y n (i) -f n,Dn (y(i)) 2 + η N m=1 E ∂f n,Dn (y(i)) ∂y m (i) 2 . ( 11
)
where the expectations are w.r.t. y(i) and conditioned on the dictionary (8).

Combining [START_REF] Djuric | Cooperative and Graph Signal Processing: Principles and Applications[END_REF] with [START_REF] Tremblay | Accelerated spectral clustering using graph filtering of random signals[END_REF] as shown in Appendix B, we get:

min γ J(γ) = 1 2 E y n (i) -γ s(i) 2 + η N m=1 γ R tt,m (i) γ , (12)
with

s(i) = z(i) k(i) , γ = β α , and t m (p) = m (p) ζ m (p) ,
where:

R tt,m (i) E{t m (i)t m (i)} (13)
k(i) = col κ(y(i), y(ω q )) |Dn| q=1 (14)
z(i) = z 1 (i), . . . , z N (i) (15)
[z m (i)] q = ∂κ(y(i), y(ω q )) ∂y m (ω q ) q=1,...,|Dn| (16)
ζ(i) = ζ 1 (i), . . . , ζ N (i) (17)
[ζ m (i)] q = ∂κ(y(i), y(ω q )) ∂y m (i) q=1,...,|Dn| (18)
m (i) = 1,m (i), . . . , N,m (i) (19)
[ m1,m2 (i)] q = ∂ 2 κ(y(i), y(ω q )) ∂y m1 (ω q )∂y m2 (i) q=1,...,|Dn| (20)
Note that γ is now the optimization variable in problem [START_REF] Moscu | Online kernel-based graph topology identification with partial-derivative-imposed sparsity[END_REF]. Quantities ( 14)-( 20) can be computed in closed form when an explicit expression for kernel κ(•, •) is provided. Expressions in the case of the Gaussian kernel are given in Appendix C. A closed form expression of R tt,m for Gaussian random variables is given in the supplementary material.

Calculating the gradient of ( 12) leads to:

∇J(γ) = -E s(i) y n (i) -s (i)γ + η N m=1 R tt,m (i)γ γ R tt,m (i) γ . ( 21
)
Approximating the first expectation in ( 21) by the instantaneous estimate s(i) y n (i) -s (i)γ(i) , and R tt,m (i) by a sample covariance matrix Rtt,m (i) available at time i, allows us to devise the following subgradient descent algorithm to iteratively minimize J(γ):

γ(i + 1) = γ(i) + µs(i) y n (i) -s (i)γ(i) -µη N m=1 Rtt,m (i)γ(i) γ (i) Rtt,m (i) γ(i) , (22)
with µ a small positive step size. We set x/0 0 to write the subgradient when the denominator in the last term of ( 22) equals 0. For the sake of brevity, we introduce the notation:

∆ m (i) = γ (i) Rtt,m (i) γ(i) . (23)
The covariance matrix Rtt,m (i) needs to be updated sequentially in environments where the data arrives sequentially. An extensive body of literature concerning sequential estimators exists. The following sample covariance matrix provides such an estimate:

Rtt,m (i) = α Rtt,m (i -1) + (1 -α) t m (i)t m (i) , (24)
with α ∈ [0, 1) a forgetting factor. In the absence of hypothesis about the covariance matrix, a successful approach so far has been shrinkage estimation [START_REF] Ledoit | Nonlinear shrinkage estimation of largedimensional covariance matrices[END_REF]. We shall not discuss these details here, which are out of the scope of this paper, but refer the reader to the literature. We shall however assume that the estimator is unbiased, that is:

E{ Rtt,m (i)} = R tt,m (i) . (25)
The time index i will be omitted for stationary variables. In ( 22)-( 23), observe that each ∆ m (i) represents an estimate of the squared norm of the partial derivative of f n with respect to y m . It thus allows us to infer the graph topology by comparing every quantity ∆ m (i) to a given threshold τ n such that, following (3), ânm (i) = 1 if ∆ m (i) ≥ τ n and 0 otherwise. These thresholds can be chosen to either obtain an estimated topology which realistically explains the studied process [START_REF] Shen | Nonlinear structural vector autoregressive models with application to directed brain networks[END_REF], or to obtain a connected graph, i.e., a graph in which there exists a path between all pairs of nodes.

IV. THEORETICAL ANALYSIS

Let v(i) = γ(i)-γ * be the weight-error vector, where γ(i) is the current estimate in [START_REF] Yang | Learning connectivity and higher-order interactions in radial distribution grids[END_REF] and γ * is the optimal coefficient vector, obtained as the solution to problem [START_REF] Moscu | Online kernel-based graph topology identification with partial-derivative-imposed sparsity[END_REF].

A. Relating the error in derivatives to the error in coefficients

One difficulty related to the iterative algorithm [START_REF] Yang | Learning connectivity and higher-order interactions in radial distribution grids[END_REF] is that it considers the evolution of the coefficients γ, instead of the partial derivatives ∂ f n,Dn (y)/∂y m . Nevertheless, we can study the convergence of the derivatives by means of filter coefficients as will be shown in the following.

Using ( 5) and the Cauchy-Schwarz inequality:

∂ f n,Dn (y) ∂y m - ∂f * n,Dn (y) ∂y m = (26) = f n,Dn , κ ∂m (•, y) H -f * n,Dn , κ ∂m (•, y) H = f n,Dn -f * n,Dn , κ ∂m (•, y) H ≤ κ ∂m (•, y) H f n,Dn -f * n,Dn H . (27)
Moreover, we also have:

f n,Dn -f * n,Dn H = s γ -s γ * = s (γ -γ * ) ≤ s γ -γ * . ( 28
)
By taking the expectation of both sides conditioned on γ , we have:

E ∂ f n,Dn (y) ∂y m - ∂f * n,Dn (y) ∂y m γ ≤ E κ ∂m (•, y) H s γ γ -γ * . ( 29
)
where the first expectation on the r.h.s. is not related to f n,Dn and f * n,Dn but depends on the statistics of the input data and on the reproducing kernel and its derivatives. Thus the error in the estimated derivatives is bounded by the error in the coefficients γ. This allows us to study the convergence behavior of the derivatives indirectly by means of the coefficients γ.

B. Simplifying assumptions

To make the analysis of algorithm ( 22) tractable, we need to introduce some simplifying assumptions. A1) We assume that vector y(i) is stationary and Gaussiandistributed, with zero-mean and covariance matrix R y , namely:

p y(i) = (2π) -N/2 det{R y } -1/2 exp - 1 2 y (i)R -1 y y(i) .

A2)

We assume that L m = 1 for simplicity, but the same analysis can be extended for L m > 1. A3) We consider the Gaussian kernel:

κ(a, b) = exp -a -b 2 /2σ 2 , ( 30
)
due to its capacities as an universal approximator [START_REF] Liu | Kernel adaptive filtering: a comprehensive introduction[END_REF]. A4) We make the Modified Independence Assumption, that is, s(i)s (i) is statistically independent of v(i). This assumption has been successfully used in the analysis of various adaptive filtering algorithms [START_REF] Parreira | Stochastic behavior analysis of the Gaussian kernel least-mean-square algorithm[END_REF]. It has also been shown in [START_REF] Minkoff | Comment on the "Unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF] to be less restrictive when compared to the classical independence assumption. A5) The elements {y(ω p )} in the dictionary are set a priori.

They are thus independent of y(i). A6) The entries of v(i) are jointly Gaussian distributed. This hypothesis has been used and validated in the analysis of both linear [START_REF] Chen | Sparse LMS for system identification[END_REF], [START_REF] Chen | Transient performance analysis of zero-attracting LMS[END_REF] and nonlinear (kernel-based) [START_REF] Gao | Transient performance analysis of zero-attracting Gaussian kernel LMS algorithm with pre-tuned dictionary[END_REF] adaptive algorithms. We also use the following approximation

Rtt,m (i) ≈ E Rtt,m (i) (25) = R tt,m (31)
to make the analysis tractable. As an example, note that the decay of the estimation error of R tt,m using (24) with i.i.d. samples is in the order of (1/i).

C. Mean weight error analysis

The estimation error e(i) is given by:

e(i) y n (i) -f n (y(i)) = y n (i) -s (i)γ(i) = y n (i) -s (i)γ * e0(i) -s (i)v(i) . (32)
Considering ( 22) and ( 32), we can write v(i + 1) as:

v(i + 1) = v(i) -µs(i)s (i)v(i) + µs(i)e 0 (i) (33) -µη N m=1 Rtt,m (i) [v(i) + γ * ] [v(i) + γ * ] Rtt,m (i) [v(i) + γ * ]
,

where e 0 (i) = y n (i) -s (i)γ * is the optimal estimation error.

Taking the expected value of both sides of (33), and using A4, we obtain the mean weight error model:

E{v(i + 1)} = (I -µR ss )E{v(i)} + µ r sy -R ss γ * -µη N m=1 E Rtt,m (i) [v(i) + γ * ] [v(i) + γ * ] Rtt,m (i) [v(i) + γ * ] , (34)
with R ss E{s(i)s (i)} and r sy E{s(i)y n (i)}.

The last term in ( 34) is a non-trivial expectation involving v(i) and Rtt,m (i). To proceed, we need to use the simplifying assumption Rtt,m (i) ≈ R tt,m discussed above. This yields:

E    Rtt,m (i)γ(i) γ (i) Rtt,m (i)γ(i)    ≈ E    R tt,m γ(i) γ (i)R tt,m γ(i)    . ( 35
)
By successively approximating the expectation of the ratio by the ratio of the expectations, and the expectation of the square root by the square root of the expectation [63, p. 70], we have:

E    R tt,m γ(i) γ (i)R tt,m γ(i)    ≈ R tt,m µ(i) Tr{R tt,m Σ(i)} + µ (i)R tt,m µ(i) , (36)
which represents the sought-after approximation of the regularization term, with µ(i) and Σ(i) given by:

µ(i) = E γ(i) = E{v(i)} + γ * , (37)
and:

Σ(i) = E γ(i)γ (i) -µ(i) µ (i) = E{v(i)v (i)} -E{v(i)}E{v(i)} , (38)
respectively. Using these results in [START_REF] Shen | Nonlinear structural vector autoregressive models with application to directed brain networks[END_REF], we obtain the following expression for the mean error recursion:

E{v(i + 1)} ≈ (I -µR ss )E{v(i)} + µ r sy -R ss γ * -µη N m=1 R tt,m µ(i) Tr{R tt,m Σ(i)} + µ (i) R tt,m µ(i) . ( 39
)
where µ(i) and Σ(i) have to be calculated at each time i in order to make recursion (39) functional. Expressions [START_REF] Shen | Kernel-based structural equation models for topology identification of directed networks[END_REF] and ( 38) are used, respectively. Specifically, a recursive model is provided in subsection IV-E for calculating E{v(i)v (i)}.

D. Mean stability analysis

Iterating [START_REF] Shen | Nonlinear structural vector autoregressive models with application to directed brain networks[END_REF] starting from i = 0, we arrive to the following expression:

E{v(i + 1)} = (I -µR ss ) i+1 E{v(0)} -µη i =0 (I -µR ss ) i- × N m=1 E    Rtt,m ( ) [v( ) + γ * ] [v( ) + γ * ] Rtt,m ( ) [v( ) + γ * ]    + µ i =0 (I -µR ss ) i-(r sy -R ss γ * ) , (40)
where E{v(0)} is the initial condition. E{v(i + 1)} converges when i → ∞ if, and only if, all terms on the r.h.s. of [START_REF] Moscu | Convergence analysis of the graph-topology-inference kernel LMS algorithm[END_REF] converge to finite values. The first term converges to zero as i → ∞ if the matrix (I -µR ss ) is stable. A sufficient condition to ensure the stability of (I -µR ss ) is to choose the step-size µ according to:

0 < µ < 2 λ max (R ss ) . ( 41
)
where λ max (•) denotes the maximum eigenvalue of its matrix argument. We shall now prove the convergence of the second and third series on the r.h.s. of [START_REF] Moscu | Convergence analysis of the graph-topology-inference kernel LMS algorithm[END_REF]. For compactness, we write the second series as i =0 (I -µR ss ) i- N m=1 q m ( ). To prove the convergence of this series, it is sufficient to prove that the series i =0 N m=1 [(I -µR ss ) i-q m ( )] k converges for all k, where [•] k is the k-th entry of its vector argument. A series is absolutely convergent if each term of the series can be bounded by a term of an absolutely convergent series. Using Jensen's inequality, we have:

[(I -µR ss ) i-q m ( )] k ≤ (I -µR ss ) i-q m ( ) ≈ (I -µR ss ) i-E    R tt,m [v( ) + γ * ] [v( ) + γ * ] R tt,m [v( ) + γ * ]    ≤ I -µR ss i- E R tt,m v( ) + γ * v( ) + γ * Rtt,m ≤ I -µR ss i-τ m R tt,m , (42)
where v Rtt,m denotes v R tt,m v. The last inequality follows because all norms are equivalent in finite dimensional spaces, that is, there exists a τ m > 0 such that:

v ≤ τ m v Rtt,m , (43)
for all v. The series i =0 I -µR ss i-N m=1 τ m R tt,m is absolutely convergent if the step-size µ is chosen according to [START_REF] Shen | Online identification of directional graph topologies capturing dynamic and nonlinear dependencies[END_REF]. The series

i =0 N m=1 [(I -µR ss ) i-q m ( )] k is therefore absolutely convergent.
To show that the third series in [START_REF] Moscu | Convergence analysis of the graph-topology-inference kernel LMS algorithm[END_REF] converges, it is sufficient to observe that:

[(I -µR ss ) i-r sy -R ss γ * ] k ≤ I -µR ss i- r sy -R ss γ * . ( 44
)
The series

i =0 I -µR ss i-
r sy -R ss γ * is absolutely convergent if the step-size µ is chosen according to [START_REF] Shen | Online identification of directional graph topologies capturing dynamic and nonlinear dependencies[END_REF].

To conclude, all series in (40) converge if ( 41) is satisfied. Therefore, for any initial condition, the algorithm (22) converges in the mean if the step-size µ is chosen according to condition [START_REF] Shen | Online identification of directional graph topologies capturing dynamic and nonlinear dependencies[END_REF].

E. Mean square error analysis

We examine now the mean-square-error behavior of the algorithm [START_REF] Yang | Learning connectivity and higher-order interactions in radial distribution grids[END_REF] by studying E{v(i)v (i)}. The mean-squareerror recursion is given by:

V (i + 1) = E{v(i + 1)v (i + 1)} = E sym 1 2 v(i)v (i) Q 1 -µ v(i)v (i)s(i)s (i) Q 2 + µ v(i)s (i)e 0 (i) Q 3 -µ 2 s(i)s (i)v(i)s (i)e 0 (i) Q 4 + µ 2 2 s(i)s (i)(e 0 (i)) 2 Q 5 + µ 2 2 s(i)s (i)v(i)v (i)s(i)s (i) Q 6 -µη v(i) N m=1 Rtt,m (i)γ(i) ∆ m (i) Q 7 + µ 2 η s(i)s (i)v(i) N m=1 Rtt,m (i)γ(i) ∆ m (i) Q 8 -µ 2 η s(i)e 0 (i) N m=1
Rtt,m (i)γ(i) ∆ m (i)

Q 9 + µ 2 η 2 2 N m=1 Rtt,m (i)γ(i) ∆ m (i) N p=1 Rtt,p (i)γ(i) ∆ p (i) Q 10 , (45)
where sym(X) = X + X . We shall now calculate the expectations separately. Helpful details are provided in [START_REF] Parreira | Stochastic behavior analysis of the Gaussian kernel least-mean-square algorithm[END_REF]:

Q 1 E{v(i)v (i)} = V (i)
appears on the r.h.s. of [START_REF] Rosasco | Nonparametric sparsity and regularization[END_REF]. Q 2 Using some algebraic manipulations and assumption A4:

E{v(i)v (i)s(i)s (i)} = E{v(i)v (i)}E{s(i)s (i)} = V (i)R ss , (46)
where R ss is calculated in closed-form in the supplementary material for this paper. Q 3 Using assumption A4, we obtain:

E{v(i)s (i)e 0 (i)} = E{v(i)}E{s (i)e 0 (i)} = E{v(i)} (r sy -R ss γ * ) , (47)
where all these quantities have already been calculated.

Q Considering hypothesis A4 and assuming the input signal to be i.i.d., we have:

E{s(i)s (i)v(i)s (i)e 0 (i)} u,v = N |Dn| a=1 E{s u (i)s a (i)v a (i)s v (i)e 0 (i)} = N |Dn| a=1 E{v a (i)}E{s u (i)s a (i)s v (i)[y n (i) -s (i)γ * ]} = - N |Dn| b=1 N |Dn| a=1 E{v a (i)}E{s u (i)s a (i)s v (i)s b (i)}γ * b + N |Dn| a=1 E{v a (i)}E{s u (i)s a (i)s v (i)y n (i)} , (48)
where both expectations E{s u (i)s a (i)s v (i)y n (i)} and E{s u (i)s a (i)s v (i)s b (i)} are provided in closed-form in the supplementary material for this paper. Q This term can be written as:

E{s(i)s (i)e 2 0 (i)} u,v = E{s u (i)s v (i)[y n (i) -s (i)γ * ] 2 } = E{s u (i)s v (i)y 2 n (i)} -2 p γ * p E{s u (i)s v (i)s p (i)y n (i)} + m γ * γ * m E{s u (i)s v (i)s (i)s m (i)} . (49)
Q We use the strategy in [START_REF] Parreira | Stochastic behavior analysis of the Gaussian kernel least-mean-square algorithm[END_REF] and assumption A4 in order to obtain the following approximation:

E{s(i)s (i)v(i)v (i)s(i)s (i)} u,v = m E{s u (i)s (i)v (i)v m (i)s m (i)s v (i)} ≈ m E{s u (i)s (i)s m (i)s v (i)}E{v (i)v m (i)} , (50)
where

E{v (i)v m (i)} = V (i) ,m . A closed-form ex- pression of the expectation E{s u (i)s (i)s m (i)s v (i)} is provided in the supplementary material for this paper. Q Recalling that v(i) = γ(i) -γ * , this term becomes: E v(i) N m=1 Rtt,m (i)γ(i) ∆ m (i) = N m=1 E γ(i)γ (i) Rtt,m (i) ∆ m (i) - N m=1 γ * E γ (i) Rtt,m (i) ∆ m (i) . ( 51
)
The second term on the r.h.s. of ( 51) can be calculated with [START_REF] De Zwart | Hemodynamic nonlinearities affect BOLD fMRI response timing and amplitude[END_REF]. To calculate the first term on the r.h.s. of ( 51), we use the following approximation:

E γ(i)γ (i) Rtt,m (i) ∆ m (i) ≈ E γ(i)γ (i)R tt,m ∆ m (i) ≈ Σ(i) + µ(i) µ (i) R tt,m Tr{R tt,m Σ(i)} + µ (i)R tt,m µ(i) , (52)
Q 8 By using some algebraic manipulations and A4, we get:

E s(i)s (i)v(i) N m=1 Rtt,m (i)γ(i) ∆ m (i) = E s(i)s (i) E v(i) N m=1 Rtt,m (i)γ(i) ∆ m (i) = R ss Q 7 . (53)
Q 9 By approximating Rtt,m (i) ≈ R tt,m and using assumption A4, we obtain:

E s(i)e 0 (i) N m=1 Rtt,m (i)γ(i) ∆ m (i) = E s(i)e 0 (i) E N m=1
Rtt,m (i)γ(i)

γ (i) Rtt,m (i)γ(i) = (r sy -R ss γ * ) (54)
× N m=1 R tt,m µ(i) Tr{R tt,m Σ(i)} + µ (i)R tt,m µ(i) .
Q 10 We successively have:

E N m=1 Rtt,m (i)γ(i) ∆ m (i) N p=1
Rtt,p (i)γ(i)

∆ p (i) ≈ N m=1 N p=1 E R tt,m γ(i)γ (i)R tt,p γ (i)R tt,m γ(i)γ (i)R tt,p γ(i) ≈ N m=1 N p=1 R tt,m Σ(i) + µ(i) µ (i) R tt,p E γ (i)R tt,m γ(i)γ (i)R tt,p γ(i) . (55)
Considering hypothesis A6, the denominator of this expression consists of the expectation of a linear combination of fourth order moments of Gaussian random variables. Using [START_REF] Kumar | Expectation of product of quadratic forms[END_REF], we have:

E γ (i)R tt,m γ(i)γ (i)R tt,p γ(i) = 4µ R tt,m ΣR tt,p µ + 2Tr R tt,m ΣR tt,p Σ + µ R tt,m µ + Tr R tt,m Σ × µ R tt,p µ + Tr R tt,p Σ , (56)
To conclude, using the results from relations ( 46)-( 55), we arrive at the following recursive equation for V (i):

V (i + 1) = V (i) -µ sym{V (i)R ss } + µ sym{Q 3 } -µ 2 sym{Q 4 } + µ 2 Q 6 -µη sym{Q 7 } + µ 2 η 2 Q 10 + µ 2 Q 5 + µ 2 η sym{R ss Q 7 } -µ 2 η sym{Q 9 } . (57)
The unknown moments are calculated in closed-form in the supplementary material for this paper.

Before going further, we define the Mean Square Deviation:

MSD(i) E γ(i) -γ * 2 = E v(i) 2 = Tr {V (i)} . ( 58
)
This quantity is useful in determining the performance of the algorithm.

F. Steady-state performance when η = 0

Assuming a small enough step-size µ that verifies [START_REF] Shen | Online identification of directional graph topologies capturing dynamic and nonlinear dependencies[END_REF], with relation [START_REF] Moscu | Convergence analysis of the graph-topology-inference kernel LMS algorithm[END_REF] we have:

v(∞) = lim i→∞ E {v(i)} = 0 . ( 59
)
We shall now evaluate V (∞) = lim i→∞ V (i) using [START_REF] Liu | Kernel adaptive filtering: a comprehensive introduction[END_REF]. With v(∞) = 0, we know that Q 3 (∞) = Q 4 (∞) = 0. Using vectorization of (57) together with the Kronecker product gives:

vec{V (∞)} = µ 2 (I 2 -F 0 ) -1 vec{Q 5 } , (60)
with:

F 0 = I 2 -µ(I ⊗ R ss + R ss ⊗ I) + µ 2 F 1 , (61)
where operator ⊗ is the Kronecker product. In these expressions, I 2 is the identity matrix of size k 2 s × k 2 s , and I is the identity matrix of size k s × k s , where k s = (N + 1)|D n | is the number of entries in block vector s. The entries of F 1 are:

[F 1 ] u+( -1)ks,m+(v-1)ks = E{s u (i)s (i)s m (i)s v (i)} . (62)
They are provided in closed-form in the supplementary material for this paper.

Finally, we conclude that:

MSD(∞) = Tr vec -1 µ 2 (I 2 -F 0 ) -1 vec{Q 5 } . (63)

V. EXPERIMENTAL VALIDATION

The performance of algorithm ( 22) has already been illustrated in the preliminary work [START_REF] Moscu | Online kernel-based graph topology identification with partial-derivative-imposed sparsity[END_REF], both with simulated and real-world data. Consequently, the experiments in this section will be restricted to model validation.

A. Signal model

We considered a simulation scenario with i.i.d. zero-mean Gaussian data y(i) with covariance matrix R yy satisfying the following model:

y(i) = Ay(i) + v(i) , ( 64
)
where A is the adjacency matrix defined as:

A =       0 1 0 1 1 1 0 1 0 1 1 0 0 1 0 0 1 1 0 1 1 0 1 1 0       . ( 65
)
and v(i) is a zero-mean white Gaussian noise with covariance matrix σ 2 v I 5 and σ v = 0.05. Note that R yy satisfies:

R yy = [I -A] -1 R vv ([I -A] -1 ) . (66)
This model, although not nonlinear, offers exact knowledge of the statistical properties necessary to validate the models derived in this paper. Validation on a purely nonlinear system follows in subsection V-E. The reader is also referred to [START_REF] Moscu | Online kernel-based graph topology identification with partial-derivative-imposed sparsity[END_REF] to assess the methods' behavior in real-world nonlinear settings. We used our algorithm with the Gaussian kernel [START_REF] Dong | Learning graphs from data: A signal representation perspective[END_REF], with kernel bandwidth σ = 1. Each node stored a dictionary D n with 6 elements chosen uniformly in [-1, 1] 4 . Simulations were averaged over 100 Monte-Carlo runs. Results presented hereafter focus on node n = 1.

B. Computing the optimal coefficients

In order to compute the optimal coefficients, we note that problem ( 12) can be written equivalently as:

min γ 1 2 γ R ss γ -γ r sy + η N m=1 C m γ 2 , ( 67
)
where C m C m = R tt,m . We solved this optimization problem using the CVX package to get γ * [START_REF] Grant | CVX: Matlab software for disciplined convex programming, version 2.1[END_REF], [START_REF]Graph implementations for nonsmooth convex programs[END_REF].

The MSD was evaluated experimentally using:

MSD(i) = E γ(i) -γ * 2 . (68)
C. Parameter η = 0

Figure 1(a) shows the theoretical and experimental learning curves of the entries of the parameter vector γ. Figure 1(b) shows both the theoretical and experimental MSD curves, as well as the steady-state MSD. The theoretical curves closely follow the theoretical ones, both in the analyses in the mean and mean square sense. They also validate the approximation (36) of the regularization term which is in the update [START_REF] Shen | Nonlinear structural vector autoregressive models with application to directed brain networks[END_REF].

D. Parameter η > 0

The conditions of this experiment were the same as for the previous case, with η = 1 • 10 -4 .

The learning curves are presented in Fig. 2. As for the previous case, we observe that the theoretical curves closely fit the experimental ones. This validates the assumptions and approximations employed in the analysis to make it tractable.

E. Validation with a nonlinear model

In order to evaluate how well the model can predict the the behavior of the algorithm in nonlinear settings, we consider an additional example with a 3-node graph with adjacency matrix given by:

A =   0 0 1 1 0 1 1 0 0   , (69)
y(i) =    y 1 -K 1 (y 3 + y 1 ) 3 (K 2 y 1 ) -1 y 2 + y2-K1(y3+y1) 3 (K2y1) -1 [0.5+exp(K1(y3+y1) 3 (K2y1) -1 )] 5 +1 y 3 + y 1 + K 1 (y 3 + y 1 ) 3 (K 2 y 1 ) -1    f (y(i))

+ (i)

with K 1 = 8000 and K 2 = 27, (i) being an i.i.d. zeromean Gaussian random vector with covariance matrix I 3 , and

f (y(i)) = f 1 (y(i)), f 2 (y(i)), f 3 (y(i))
the vector-valued function corresponding to the concatenation of f n in [START_REF] Heiberger | Predicting economic growth with stock networks[END_REF]. Notice that the sample index i was omitted inside f (y(i)) to simplify the notation. Graph signal y(i) can be sampled from this model by first sampling from (i) and using the relation y(i) = (id -f ) -1 ( (i)), with id(x) the identity function. We used our algorithm with the Gaussian kernel [START_REF] Dong | Learning graphs from data: A signal representation perspective[END_REF], with kernel bandwidth σ = 1. Each node stored a dictionary D n with 4 elements chosen uniformly in [-1, 1] × [-1, 1]. Simulations were averaged over 100 Monte-Carlo runs. Results presented hereafter focus on node n = 1.

Figure 3 depicts the behavior of the algorithm in mean and mean-square sense without regularization (i.e., for η = 0). Figure 4 shows the mean and mean-square performance in the case when η = 0.3. These two figures together show that the conducted analysis is indeed able to predict the behavior of the algorithm without the need of running multiple Monte-Carlo runs. Moreover, even with the highly nonlinear nature of f and y(i) not being Gaussian distributed, the theoretical curves still follow the experimental ones relatively closely, although a small amount of bias can be observed in steady state.

VI. CONCLUSION

In this paper, an online graph-topology inference method recently proposed by the authors was recalled and analyzed. While previous works mainly focused on models based on additive interactions between the signals at each node, we considered arbitrary nonlinear interactions between the nodes, which render our model much more general. By encoding the links as partial derivatives of some nonlinear functions, our method benefits from the kernel machinery framework to estimate a possibly directed, sparse adjacency matrix. However, it also demonstrates significantly increased complexity that makes its analysis non-trivial. This work proposed a thorough analysis of the algorithm, as well as performance bounds and stability conditions. The derived models characterize the performance of the algorithm as a function of important parameters such as the kernel bandwidth or the step-size. In turn, this allows for precise tuning of such parameters in order to obtain a desired performance in transient or in steady state. Moreover, the behavior of the intuitive partial-derivative-based sparsity regularizer was analyzed in-depth, thus completely characterizing the method. The experimental results, as well as the algorithm analysis, show that the proposed method can lead to more accurate estimates for more general nonlinear systems.

APPENDIX A PROOF OF THEOREM 1

Proof. Using orthogonal projection, we can decompose any function f n ∈ H κ as the sum of two functions:

f n = f n + f ⊥ n , (70)
where

f ⊥ n is orthonormal to f n , that is, f n , f ⊥ n H = 0)
, and f n can be written as:

f n = i =1 α κ (•, y( )) + N m=1 i p=1 Lm q=1 β m,p,q κ ∂m,q (•, y(p)) , (71)
where α and β m,p,q are coefficients. This means that f n lies in the span of κ (•, y( )) and κ ∂m,q (•, y(p)), from which the orthogonality condition implies that f ⊥ n , κ (•, y( )) H = 0 and f ⊥ n , κ ∂m,q (•, y(p)) H = 0, for all , m, p, q. Assume that the kernel κ(•, •) is at least twice differentiable. Then, the following relation holds [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF]:

H ∂f n (y) ∂y m,q = f n , κ ∂m,q (•, y) H , (72)
where y m,q is the q-th entry of y m . Let us plug then decomposition (70) of f n in problem (4). For the first term, using the reproducing property, we have:

f n (y( )) = f n , κ(•, y( )) H (73) = f n + f ⊥ n , κ(•, y( )) H (74) = f n , κ(•, y( )) H , (75)
for all = 1, . . . , i, where

f ⊥ n , κ(•, y( )) H = 0 because f ⊥ n
is orthogonal to each term in (71). For the second term:

∂f n (y(p)) ∂y m,q = ∂f n (y) ∂y m,q y=y(p) (76)
= f n + f ⊥ n , κ ∂m,q (•, y(p)) H (77) = f n , κ ∂m,q (•, y(p)) H , (78)
for all p = 1, . . . , i, where f ⊥ n , κ ∂m,q (•, y(p)) H = 0 since f ⊥ n is perpendicular to each term in (71). For the last term, we have:

ψ f n H = ψ f n + f ⊥ n H (79) = ψ f n H + f ⊥ n H . (80)
Since f ⊥ n H does not influence the first two terms in problem (4), if ψ is monotonically increasing, then the solution is such that f ⊥ n H = 0.

APPENDIX B COMPUTING THE TERMS IN THE OPTIMIZATION PROBLEM

THAT DEPEND ON f n,Dn Let us recall the optimization problem [START_REF] Djuric | Cooperative and Graph Signal Processing: Principles and Applications[END_REF]:

min f n,Dn 1 2 E y n (i) -f n,Dn (y(i)) 2 + η N m=1 E ∂f n,Dn (y(i)) ∂y m (i) 2 . ( 81
)
We now evaluate the terms involving f n,Dn , using the reproducing property and (5). For the first term, we have:

f n,Dn (y(i)) = f n,Dn , κ(•, y(i)) H (82) = |Dn| p=1 α p κ(•, y(ω p )), κ(•, y(i)) H + N m=1 |Dn| q=1 β m,q κ ∂m (•, y(ω q )), κ(•, y(i)) H ,
from which:

κ(•, y(ω p )), κ(•, y(i)) H = κ(y(ω p ), y(i)) , κ ∂m (•, y(ω q )), κ(•, y(i)) H = ∂κ(y(ω q ), y(i)) ∂y m (ω q ) . (83a) (83b)
For the second term, we have:

∂f n,Dn (y(i)) ∂y m2 (i) = f n,Dn , κ ∂m 2 (•, y(i)) H = |Dn| p=1 α p κ(•, y(ω p )), κ ∂m 2 (•, y(i)) H + N m2=1 |Dn| q=1 β m2,q κ ∂m 2 (•, y(ω q )) , κ ∂m 2 (•, y(i)) H , (84)
from which: 16)-( 20) are given by:

κ(•, y(ω p )), κ ∂m 2 (•, y(i)) H = ∂κ(y(ω p ), y(i)) ∂y m2 (i) , κ ∂m 1 (•, y(ω q )), κ ∂m 2 (•, y(i)) H = ∂κ(y(ω q ), y(i)) ∂y m1 (ω q )y m2 (i) . ( 85a
) (85b) APPENDIX C CALCULATION OF z m (i), ζ m (i), AND m1,m2 (i) With the Gaussian kernel κ(a, b) = exp -a -b 2 /2σ 2 , quantities (
[z m (i)] q = κ(y(i), y(ω q )) y m (i) -y m (ω q ) σ 2 [ζ m (i)] q = -[z m (i)] q , [ m1,m2 (i)] q =        -κ(y(i), y(ω q )) y m1 (i) -y m1 (ω q ) σ 2 y m2 (i) -y m2 (ω q ) σ 2 , m 1 = m 2 -κ(y(i), y(ω q )) y m1 (i) -y m1 (ω q ) 2 σ 4 - 1 σ 2 , m 1 = m 2
.

Note that we can write:

- 1 2 y c 4 B 0 + R -1 y y -x 4 B 4 y = - 1 2 y c 4 B 0 + R -1 y y + x 4 B 4 y + y B 4 x 4 = - 1 2 y + c 4 B 0 + R -1 y -1 B 4 x 4 c 4 B 0 + R -1 y y + c 4 B 0 + R -1 y -1 B 4 x 4 + 1 2 x 4 B 4 c 4 B 0 + R -1 y -1 B 4 x 4 . (9)
Thus, relation [START_REF] Koppel | Decentralized online learning with kernels[END_REF] becomes:

= (2π) - N + 1 2 det{R y } - 1 2 exp 1 2 x 4 Q 4 x 4 R • • • R ( y h1 -x h2 ) ι1 ( y h3 -x h4 ) ι2 ( y h5 -x h6 ) ι3 ( y h7 -x h8 ) ι4 × exp - 1 2 y + c 4 B 0 + R -1 y -1 B 4 x 4 c 4 B 0 + R -1 y y + c 4 B 0 + R -1 y -1 B 4 x 4 × exp 1 2 x 4 B 4 c 4 B 0 + R -1 y -1 B 4 x 4 d y 1 • • • d y N +1 = (2π) - N + 1 2 det{R y } - 1 2 exp 1 2 x 4 Q 4 x 4 det c 4 B 0 + R -1 y -1 1 2 det c 4 B 0 + R -1 y -1 - 1 2 × exp 1 2 x 4 B 4 c 4 B 0 + R -1 y -1 B 4 x 4 R • • • R ( y h1 -x h2 ) ι1 ( y h3 -x h4 ) ι2 ( y h5 -x h6 ) ι3 ( y h7 -x h8 ) ι4 × exp - 1 2 y + c 4 B 0 + R -1 y -1 B 4 x 4 c 4 B 0 + R -1 y y + c 4 B 0 + R -1 y -1 B 4 x 4 d y 1 • • • d y N +1 = det{R y } - 1 2 det c 4 B 0 + R -1 y -1 1 2 exp 1 2 x 4 Q 4 + B 4 c 4 B 0 + R -1 y -1 B 4 x 4 × E p ( y) ( y h1 -x h2 ) ι1 ( y h3 -x h4 ) ι2 ( y h5 -x h6 ) ι3 ( y h7 -x h8 ) ι4 (10) = ν({h i } 8 i=1 ) , (11)
with the normal distribution:

p ( y) = N µ = -c 4 B 0 + R -1 y -1 B 4 x 4 , Σ = c 4 B 0 + R -1 y -1 . (12)
Let us focus on the expectation in [START_REF] Tremblay | Accelerated spectral clustering using graph filtering of random signals[END_REF]. We have, for ι i = 1, i = {1, 2, 3, 4}:

E p (y) ( y h1 -x h2 ) ι1 ( y h3 -x h4 ) ι2 ( y h5 -x h6 ) ι3 ( y h7 -x h8 ) ι4 = E y h1 y h3 y h5 y h7 -x h8 E y h1 y h3 y h5 -x h6 E y h1 y h3 y h7 + x h6 x h8 E y h1 y h3 -x h4 E y h1 y h5 y h7 + x h4 x h8 E y h1 y h5 + x h4 x h6 E y h1 y h7 -x h4 x h6 x h8 E y h1 -x h2 E y h3 y h5 y h7 + x h2 x h8 E y h3 y h5 + x h2 x h6 E y h3 y h7 -x h2 x h6 x h8 E y h3 + x h2 x h4 E y h5 y h7 -x h2 x h4 x h8 E y h5 -x h2 x h4 x h6 E y h7 + x h2 x h4 x h6 x h8 . (13)
Since h i , i ∈ {1, 3, 5, 7} represent solely generic placeholders for any actual index, it suffices to compute only a few expectations from the previous relation, which can then be used for any other combination. We have:

E p ( y) y h1 y h3 y h5 y h7 = µ h1 µ h3 µ h5 µ h7 + Σ h1,h3 Σ h5,h7 + Σ h1,h5 Σ h3,h7 + Σ h1,h7 Σ h3,h5

+ µ h1 µ h3 Σ h5,h7 + µ h1 µ h5 Σ h3,h7 + µ h1 µ h7 Σ h3,h5 + µ h3 µ h5 Σ h1,h7 + µ h3 µ h7 Σ h1,h5 + µ h5 µ h7 Σ h1,h3 , (14)
E p (y) y h1 y h3 y h5 = µ h1 µ h3 µ h5 + µ h5 Σ h1,h3 + µ h3 Σ h1,h5 + µ h1 Σ h3,h5 , (15)
E p (y) y h1 y h3 = µ h1 µ h3 + Σ h1,h3 , (16)
E p (y) y h1 = µ h1 . ( 17
) APPENDIX D CASES CORRESPONDING TO R tt,m (c 1 = c 2 = 0, c 3 = 1) Since t m (p) = m (p) ζ m (p)
, the matrix R tt,m can be written as:

R tt,m = E{t m (i)t m (i)} = E m (i) ζ m (i) m (i) ζ m (i) = R ,m R ζ ,m R ζ ,m R ζζ,m , (18)
where in the last step we used the fact that the signals are identically distributed (for different i).

For ease of comprehension, we recall the definitions of z(i), ζ(i) and m (i):

z(i) = z 1 (i), . . . , z N (i) , [z m (i)] q = ∂κ(y(i), y(ω q )) ∂y m (ω q ) q=1,...,|Dn| ,
ζ(i) = ζ 1 (i), . . . , ζ N (i) , [ζ m (i)] q = ∂κ(y(i), y(ω q )) ∂y m (i) q=1,...,|Dn|
, m (i) = 1,m (i), . . . , N,m (i) , [ m1,m2 (i)] q = ∂ 2 κ(y(i), y(ω q )) ∂y m1 (ω q )∂y m2 (ω q ) q=1,...,|Dn| .

Their explicit forms for the particular case of the Gaussian kernel are:

[z m (i)] q = -[ζ m (i)] q = exp - 1 2σ 2 y(i) -y(ω q ) 2 [y(i)] m -[y(ω q )] m σ 2 , [ m1,m2 ] q =          -exp - 1 2σ 2 y(i) -y(ω q ) 2 [y(i)] m1 -[y(ω q )] m1 σ 2 [y(i)] m2 -[y(ω q )] m2 σ 2 m 1 = m 2 -exp - 1 2σ 2 y(i) -y(ω q ) 2 [y(i)] m1 -[y(ω q )] m1 2 σ 4 - 1 σ 2 m 1 = m 2 . A. Computing R ,m R ,m (a-1)|Dn|+p, (b-1)|Dn|+q = E{ m (i) m (i)} (a-1)|Dn|+p, (b-1)|Dn|+q = E [ a,m ] p [ b,m ] q =                                                              E (y a (i) -y a (ω p )) σ 2 (y m (i) -y m (ω p )) σ 2 (y b (i) -y b (ω q )) σ 2 (y m (i) -y m (ω q )) σ 2 × exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 , a = m, b = m E (y a (i) -y a (ω p )) σ 2 (y m (i) -y m (ω p )) σ 2 y b (i) -y b (ω q ) 2 σ 4 - 1 σ 2 × exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 , a = m, b = m E (y b (i) -y b (ω q )) σ 2 (y m (i) -y m (ω q )) σ 2 y a (i) -y a (ω p ) 2 σ 4 - 1 σ 2 × exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 , a = m, b = m E y a (i) -y a (ω p ) 2 σ 4 - 1 σ 2 y b (i) -y b (ω q ) 2 σ 4 - 1 σ 2 × exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 , a = m, b = m , (19)
h 1 = a h 2 = a h 3 = m h 4 = m h 5 = b h 6 = N + b h 7 = m h 8 = N + m , when a = m, b = m 1 σ 8 ν({h i } 8 i=1 ) - 1 σ 2 R zz (a-1)|Dn|+p,(m-1)|Dn|+p , with                              h 1 = a h 2 = a h 3 = b h 4 = b h 5 = b h 6 = N + b h 7 = b h 8 = N + b , when a = m, b = m 1 σ 8 ν({h i } 8 i=1 ) - 1 σ 2 R zz (b-1)|Dn|+q,(m-1)|Dn|+q , with                              h 1 = b h 2 = b h 3 = a h 4 = a h 5 = a h 6 = N + a h 7 = a h 8 = N + a , when a = m, b = m 1 σ 8 ν({h i } 8 i=1 ) - 1 σ 2 R zz (a-1)|Dn|+p,(a-1)|Dn|+p - 1 σ 2 R zz (b-1)|Dn|+q,(b-1)|Dn|+q + 1 σ 4 R kk p,q , with                              h 1 = a h 2 = a h 3 = a h 4 = a h 5 = b h 6 = N + b h 7 = b h 8 = N + b , when a = m, b = m . ( 20
)
h 1 = • h 2 = • ι 1 = 0 h 3 = m h 4 = m h 5 = b h 6 = N + b h 7 = m h 8 = N + m , when b = m - 1 σ 6 ν({h i } 8 i=1 ) + 1 σ 2 R kz q,(b-1)|Dn|+p , with                            h 1 = • h 2 = • ι 1 = 0 h 3 = b h 4 = b h 5 = b h 6 = N + b h 7 = b h 8 = N + b , when b = m , (22)
where the bullet • means that the index is irrelevant, due to its corresponding term not being present in the generic form [START_REF] Frazer | A second generation human haplotype map of over 3.1 million SNPs[END_REF].

C. Computing R ss (including R zz,m , R zz,m1,m2 , R kz,m and R kk,m ) Since s(i) = z(i) k(i)
, all the matrices detailed in this section can be used in computing R ss .

R ζζ,m p, q = E{ζ m (i)ζ m (i)} p, q = E [ζ m ] p [ζ m ] q = E (y m (i) -y m (ω p )) σ 2 (y m (i) -y m (ω q )) σ 2 exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 . (23) = 1 σ 4 ν({h i } 8 i=1 ), with                          h 1 = • h 2 = • ι 1 = 0 h 3 = m h 4 = m h 5 = • h 6 = • ι 3 = 0 h 7 = m h 8 = N + m . (24)
Generalizing to the case m 1 = m 2 , we straightforwardly have:

R ζζ,m1,m2 p, q = E{ζ m1 (i)ζ m2 (i)} p, q = E [ζ m1 ] p [ζ m2 ] q = E (y m1 (i) -y m1 (ω p )) σ 2 (y m2 (i) -y m2 (ω q )) σ 2 exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 . ( 25
) = 1 σ 4 ν({h i } 8 i=1 ), with                          h 1 = • h 2 = • ι 1 = 0 h 3 = m 1 h 4 = m 1 h 5 = • h 6 = • ι 3 = 0 h 7 = m 2 h 8 = N + m 2 . ( 26
)
Particularizing to R kz,m and using the fact that for the Gaussian kernel [z m ] q = -[ζ m ] q , we have:

R kz,m p, q = -E{k(i)ζ m (i)} p, q = E [k] p [z m ] q = E (y m2 (i) -y m2 (ω q )) σ 2 exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 . ( 27
) = 1 σ 2 ν({h i } 8 i=1 ), with                        h 1 = • h 2 = • ι 1 = 0 h 3 = • h 4 = • ι 2 = 0 h 5 = • h 6 = • ι 3 = 0 h 7 = m h 8 = N + m . (28)
Particularizing to R kk , we have:

R kk p, q = E{k(i)k (i)} p, q = E [k] p [k] q = E exp - 1 2σ 2 y(i) -y(ω p ) 2 - 1 2σ 2 y(i) -y(ω q ) 2 . ( 29
) = ν({h i } 8 i=1 ), with                        h 1 = • h 2 = • ι 1 = 0 h 3 = • h 4 = • ι 2 = 0 h 5 = • h 6 = • ι 3 = 0 h 7 = • h 8 = • ι 4 = 0 . ( 30
) APPENDIX E CASES CORRESPONDING TO E{s u (i)s v (i)s a (i)s b (i)} (c 1 = c 2 = c 3 = 1) Recall that s(i) = z(i) k(i)
. We have the following terms: + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 + y(i) -y(ω s ) 2 (35) + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 + y(i) -y(ω s ) 2 = exp -1 2σ 2 y(i) -y(ω p ) 2 + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 + y(i) -y(ω s ) 2

a) Term E{z u (i)z a (i)z b (i)z v (i)}: u, v = 1, . . . , N |D n |, a, b = 1, . . . , N |D n | E{z u (i)z a (i)z a (i)z v (i)} = E (y m1 (i) -y m1 (ω p )) σ 2 (y m2 (i) -y m2 (ω q )) σ 2 (y m3 (i) -y m3 (ω r )) σ 2 (y m4 (i) -y m4 (ω s )) σ 2 × exp - 1 2σ 2 y(i) -y(ω p ) 2 + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 + y(i) -y(ω s ) 2 (31) = 1 σ 8 ν({h i } 8 i=1 ), with                h 1 = h 2 = m 1 = u |Dn| h 3 = h 4 -N = m 2 = a |Dn| h 5 = h 6 -2N = m 3 = b |Dn| h 7 = h 8 -3N = m 4 = v |Dn| ,          ω p = mod(u -1, |D n |) + 1 ω q = mod(a -1, |D n |) + 1 ω r = mod(b -1, |D n |) + 1 ω s = mod(v -1, |D n |) + 1 . ( 32
) b) Term E{k u (i)z a (i)z b (i)z v (i)}: u = 1, . . . , |D n |, v = 1, . . . , N |D n |, a, b = 1, . . . , N |D n | E{k u (i)z a (i)z b (i)z v (i)} = E (y m2 (i) -y m2 (ω q )) σ 2 (y m3 (i) -y m3 (ω r )) σ 2 (y m4 (i) -y m4 (ω s )) σ 2 exp - 1 2σ 2 y(i) -y(ω p ) 2 + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 + y(i) -y(ω s ) 2 (33) = 1 σ 6 ν({h i } 8 i=1 ), with              h 1 = h 2 = m 1 = • h 3 = h 4 -N = m 2 = a |Dn| h 5 = h 6 -2N = m 3 = b |Dn| h 7 = h 8 -3N = m 4 = v |Dn| ,          ω p = u ω q = mod
= 1 σ 4 ν({h i } 8 i=1 ), with            h 1 = h 2 = m 1 = • h 3 = h 4 -N = m 2 = • h 5 = h 6 -2N = m 3 = b |Dn| h 7 = h 8 -3N = m 4 = v |Dn| ,          ω p = u ω q = a ω r = mod(b -1, |D n |) + 1 ω s = mod(v -1, |D n |) + 1 . ( 36
= ν({h i } 8 i=1 ), with

         h 1 = h 2 = m 1 = • h 3 = h 4 -N = m 2 = • h 5 = h 6 -2N = m 3 = • h 7 = h 8 -3N = m 4 = • ,          ω p = u ω q = a ω r = b ω s = v . ( 40
)
APPENDIX F CASES CORRESPONDING TO E{s u (i)s a (i)s b (i)y n (i)} (c 1 = 0, c 2 = c 3 = 1)

These terms can easily be computed using the results in section E. (y m2 (i) -y m2 (ω q )) σ 2

(y m3 (i) -y m3 (ω r )) σ 2 y n (i)

× exp -1 2σ 2 y(i) -y(ω p ) 2 + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 (41) (y m3 (i) -y m3 (ω r )) σ 2 y n (i)

= 1 σ 6 ν({h i } 8 i=1 ), with              h 1 = h 2 = m 1 = u |Dn| h 3 = h 4 -N = m 2 =
× exp -1 2σ 2 y(i) -y(ω p ) 2 + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2

= ν({h i } 8 i=1 ), with

         h 1 = h 2 = m 1 = • h 3 = h 4 -N = m 2 = • h 5 = h 6 -2N = m 3 = • h 7 = h 8 -3N = m 4 = n ,          ω p = • ω q = • ω r = b ω s = • . (58)

Fig. 1 .

 1 Fig. 1. Analysis validation in mean and mean square sense, for η = 0

Fig. 2 .

 2 Fig. 2. Analysis validation in mean and mean square sense, for η = 1 • 10 -4

Fig. 3 .

 3 Fig. 3. Analysis validation in mean and mean-square sense, for η = 0

Fig. 4 .

 4 Fig. 4. Analysis validation in mean and mean-square sense, for η = 0.3

(a - 1 , 1 . 2 (

 112 |D n |) + 1 ω r = mod(b -1, |D n |) + 1 ω s = mod(v -1, |D n |) + E{k u (i)k a (i)z b (i)z v (i)}: u = 1, . . . , |D n |, v = 1, . . . , N |D n |, a = 1, . . . , |D n |, b = 1, . . . , N |D n | E{k u (i)k a (i)z b (i)z v (i)} = E (y m3 (i) -y m3 (ω r )) σ y m4 (i) -y m4 (ω s )) σ 2 exp -1 2σ 2 y(i) -y(ω p ) 2

  ) d) Term E{k u (i)k a (i)k b (i)z v (i)}: u = 1, . . . , |D n |, v = 1, . . . , N |D n |, a, b = 1, . . . , |D n | E{k u (i)k a (i)k b (i)z v (i)} = E (y m4 (i) -y m4 (ω s )) σ 2 exp -1 2σ 2 y(i) -y(ω p ) 2

h 1 = 1 .

 11 h 2 = m 1 = • h 3 = h 4 -N = m 2 = • h 5 = h 6 -2N = m 3 = • h 7 = h 8 -3N = m 4 = p = u ω q = a ω r = b ω s = mod(v -1, |D n |) + E{k u (i)k a (i)k a (i)k v (i)}: u, v = 1, . . . , |D n |, a, b = 1, . . . , |D n | E{k u (i)k a (i)k b (i)k v (i)}

  a) Term E{z u (i)z a (i)z b (i)y n (i)}: u = 1, . . . , N |D n |, , a, b = 1, . . . , N |D n |We adapt relations (31)-[START_REF] Sandryhaila | Big data analysis with signal processing on graphs: Representation and processing of massive data sets with irregular structure[END_REF], thus obtaining:E{z u (i)z a (i)z b (i)y n (i)} = E (y m1 (i) -y m1 (ω p )) σ 2

a |Dn| h 5 = h 6 -

 56 2N = m 3 = b |Dn| h 7 = h 8 -3N = m 4 = n , p = mod(u -1, |D n |) + ω q = mod(a -1, |D n |) + ω r = mod(b -1, |D n |) + 1 ω s = • . (42) b) Term E{z u (i)k a (i)z b (i)y n (i)}: u = 1, . . . , N |D n |, a = 1, . . . , |D n |, b = 1, . . . , N |D n | We adapt relations (33)-(34), thus obtaining:E{z u (i)k a (i)z b (i)y n (i)} = E(y m1 (i) -y m1 (ω p )) σ 2

h 1 = h 2 = m 1 = u |Dn| h 3 = h 4 -N = m 2 = • h 5 = h 6 -2N = m 3 = b |Dn| h 7 = h 8 -h 1 =h 1 = h 2 = m 1 = • h 3 = h 4 -N = m 2 = • h 5 = h 6 -

 123456781123456 3N = m 4 = n , p = mod(u -1, |D n |) + ω q = a ω r = mod(b -1, |D n |) + 1 ω s = • . (44) c) Term E{z u (i)k a (i)k b (i)y n (i)}: u = 1, . . . , N |D n |, a, b = 1, . . . , |D n | We adapt relations (35)-[START_REF] De Zwart | Hemodynamic nonlinearities affect BOLD fMRI response timing and amplitude[END_REF], thus obtaining:E{z u (i)k a (i)k b (i)y n (i)} = E (y m1 (i) -y m1 (ω p )) σ 2 y n (i) × exp -1 2σ 2 y(i) -y(ω p ) 2 + y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 (45) = 1 σ 2 ν({h i } 8 i=1 ), with h 2 = m 1 = u |Dn| h 3 = h 4 -N = m 2 = • h 5 = h 6 -2N = m 3 = • h 7 = h 8 -3N = m 4 = n , p = mod(u -1, |D n |) + 1 ω q = a ω r = b ω s = • .(46) a) Term E{z b (i)y n (i)}: , b = 1, . . . , N |D n | We adapt relations (49)-(50), thus obtaining:E{z b (i)y n (i)} = E (y m3 (i) -y m3 (ω r )) σ 2 y n (i) exp -1 2σ 2 y(i) -y(ω r ) 2 2N = m 3 = y |Dn| h 7 = h 8 -3N = m 4 = n , • ω q = • ω r = mod(b -1, |D n |) + 1 ω s = • . (56) b) Term E{k y (i)y n (i)}: , b = 1, . . . , |D n | We adapt relations (53)-(54), thus obtaining: E{k b (i)y n (i)} = E y n (i) exp -1 2σ 2 y(i) -y(ω r ) 2

Note that the i index was omitted in (1) to simplify the notation.

This work was funded in part by ANR under grant ANR-19-CE48-0002. M. Moscu is with the Université Côte d'Azur, Nice, France

Supplementary material for the paper

Graph topology inference with derivative-reproducing property in RKHS: algorithm and convergence analysis

The goal of the supplementary material is to provide detailed computations of quantities R ss , R tt,m , E{s u (i)s v (i)s a (i)s b (i)}, E{s u (i)s v (i)s a (i)y n (i)}, E{s u (i)s v (i)y 2 n (i)} and E{s u (i)y n (i)}. We introduce vector y(i) = [y (i), y n (i)] , whose k-th entry is denoted by y k (i). The expectations involved in computing R tt,m (with R ss being one of its particular cases),

n (i)} and E{s u (i)y n (i)} can be expressed generically as 1 :

where x 4 = [y (ω p ), y (ω q ), y (ω r ), y (ω s )] represents the fixed dictionary elements. Also,

(3)

In the previous quantities, c 1 , c 2 , c 3 represent binary selection variables, accounting for the following possible cases:

Note that h 1 , h 2 , . . . , h 8 are different indexes of y(i) and x 4 (not necessarily distinct), and the binary variables ι i ∈ {0, 1}, for i = 1, . . . , 4, allow us to accommodate lower-order cases and be more flexible. We now have, for y(i) jointly Gaussian distributed:

= (2π)

which is equivalent to:

which becomes: 39)-( 40), thus obtaining:

These terms can easily be computed using the results in section F. 41)-( 42), thus obtaining:

We adapt relations ( 45)- [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF], thus obtaining:

E{y n (i)k a (i)k a (i)y n (i)} = E y n (i)y n (i) exp -1 2σ 2 y(i) -y(ω q ) 2 + y(i) -y(ω r ) 2 (53)

= ν({h i } 8 i=1 ), with

APPENDIX H CASES CORRESPONDING TO E{s b (i)y n (i)} (c 1 = c 2 = c 3 = 0)

These terms can easily be computed using the results in section G.