Online Distributed Learning Over Networks in RKH Spaces Using Random Fourier Features

TL;DR

This paper introduces a novel distributed online learning method in RKHS using Random Fourier Features, enabling fixed-size parameter updates and reducing communication complexity in networked environments.

Contribution

It presents the first complete protocol for distributed online learning in RKHS with fixed-size approximations using Random Fourier Features.

Findings

Achieves convergence and bounded network regret.

Demonstrates effective performance in simulations.

Reduces communication overhead in distributed settings.

Abstract

We present a novel diffusion scheme for online kernel-based learning over networks. So far, a major drawback of any online learning algorithm, operating in a reproducing kernel Hilbert space (RKHS), is the need for updating a growing number of parameters as time iterations evolve. Besides complexity, this leads to an increased need of communication resources, in a distributed setting. In contrast, the proposed method approximates the solution as a fixed-size vector (of larger dimension than the input space) using Random Fourier Features. This paves the way to use standard linear combine-then-adapt techniques. To the best of our knowledge, this is the first time that a complete protocol for distributed online learning in RKHS is presented. Conditions for asymptotic convergence and boundness of the networkwise regret are also provided. The simulated tests illustrate the performance of the proposed scheme.

Tables

Table 1. TABLE I: Dataset Information.

Method	Adult	Banana	Waveform	MNIST
Training size	32562	4000	4000	60000
Testing size	16282	1300	1000	10000
dimensions	123	2	21	784

Table 2. TABLE II: Comparing the performances of the Distributed Pegasos versus the non-cooperative Pegasos for graphs with M = 5 𝑀 5 M=5 nodes.

Method	Adult	Banana	Waveform	MNIST
Distributed-Pegasos (1)	19%	11.80%	11.82%	0.79%
Distributed-Pegasos (2)	17.43%	10.84%	10.49%	0.68%
Distributed-Pegasos (5)	15.87%	10.34%	9.56%	0.59%
Non-cooperative-Pegasos (1)	19.11%	14.52%	13.75%	1.42%
Non-cooperative-Pegasos (2)	18.31%	12.52%	12.59%	1.19%
Non-cooperative-Pegasos (5)	17.29%	11.32%	11.86%	1.01%

Table 3. TABLE III: Comparing the performances of the Distributed Pegasos versus the non-cooperative Pegasos for graphs with M = 20 𝑀 20 M=20 nodes.

Method	Adult	Banana	Waveform	MNIST
Distributed-Pegasos (1)	24.04%	16.38%	16.26%	1.03%
Distributed-Pegasos (2)	22.34%	13.23%	13.93%	0.77%
Distributed-Pegasos (5)	18.94%	10.83%	11.20%	0.57%
Non-cooperative-Pegasos (1)	20.81%	21.74%	18.40%	2.93%
Non-cooperative-Pegasos (2)	20.52%	18.64%	16.54%	2.19%
Non-cooperative-Pegasos (5)	19.88%	15.96%	14.86%	1.87%

Table 4. TABLE IV: Parameters for each method.

Method	Adult	Banana	Waveform	MNIST
Kernel-Pegasos	$\begin{array}{c}\sigma =\sqrt{10}\\ \lambda =0.0000307\end{array}$	$\begin{array}{c}\sigma =0.7\\ \lambda =\frac{1}{316}\end{array}$	$\begin{array}{c}\sigma =\sqrt{10}\\ \lambda =0.001\end{array}$	$\begin{array}{c}\sigma =4\\ \lambda ={10}^{-7}\end{array}$
RFF-Pegasos	$\begin{array}{c}\sigma =\sqrt{10}\\ \lambda =0.0000307\\ D=2000\end{array}$	$\begin{array}{c}\sigma =0.7\\ \lambda =\frac{1}{316}\\ D=200\end{array}$	$\begin{array}{c}\sigma =\sqrt{10}\\ \lambda =0.001\\ D=2000\end{array}$	$\begin{array}{c}\sigma =4\\ \lambda ={10}^{-7}\\ D=100000\end{array}$

Table 5. TABLE V: Mean training times for QKLMS and RFFKLMS.

Experiment	QKLMS time	RFFKLMS time	QKLMS dictionary size
Example 5	0.55 sec	0.35 sec	$M=1088$
Example 6	0.47 sec	0.15 sec	$M=104$
Example 7	0.02 sec	0.0057 sec	$M=7$
Example 8	0.03 sec	0.008 sec	$M=32$

Equations

$$f_{n-1}=\sum _{i=1}^{n-1}{\alpha }_i\kappa (\cdot ,{\mathbf{x}}_i),$$

$$\kappa (\mathbf{x}-\mathbf{y})=E_{\mathbf{\omega },b}[z_{\mathbf{\omega },b}(\mathbf{x})z_{\mathbf{\omega },b}(\mathbf{y})],$$

$$\kappa ({\mathbf{x}}_n-{\mathbf{x}}_m)\approx \frac{1}{D}\sum _{i=1}^D{z}_{{\mathbf{\omega }}_i,b_i}({\mathbf{x}}_m)z_{{\mathbf{\omega }}_i,b_i}({\mathbf{x}}_n).$$

$$p(\mathbf{\omega })={\left(\sigma /\sqrt{2\pi }\right)}^D{e}^{-\frac{{\sigma }^2\Vert \mathbf{\omega }{\Vert }^2}{2}},$$

$${\mathbf{z}}_{\Omega }(\mathbf{u})=\sqrt{\frac{2}{D}}\left(\begin{array}{c}\cos ({\mathbf{\omega }}_1^T\mathbf{u}+b_1)\\ ⋮\\ \cos ({\mathbf{\omega }}_D^T\mathbf{u}+b_D)\end{array}\right),$$

$$\Omega =\left(\begin{array}{cccc}{\mathbf{\omega }}_1& {\mathbf{\omega }}_2& ...& {\mathbf{\omega }}_D\\ b_1& b_2& ...& b_D\end{array}\right),$$

$$f_{n-1}({\mathbf{x}}_n)$$

$$a_{k,l}=\{\begin{array}{ll}\frac{1}{\max \{|{\mathcal{N}}_k|,|{\mathcal{N}}_l|\}},& \mathrm{if}l\in {\mathcal{N}}_k,\mathrm{and}l\ne k\\ 1-{\sum }_{i\in {\mathcal{N}}_k\setminus k}{a}_{k,i},& \mathrm{if}l=k\\ 0,& \mathrm{otherwise}.\end{array}$$

$$f_n={\psi }_{k,n-1}-{\mu }_n{\nabla }_f\mathcal{L}({\mathbf{x}}_n,y_n,{\psi }_{k,n-1}).$$

$${\mathbf{\psi }}_{k,n}$$

$${\mathbf{\theta }}_{k,n}$$

$${\underline{\mathbf{\theta }}}_n=\mathbf{A}{\underline{\mathbf{\theta }}}_{n-1}-{\mathbf{M}}_n{\mathbf{G}}_n,$$

$$\Vert \nabla \mathcal{L}({\mathbf{z}}_{\Omega }(\mathbf{x}),y,\mathbf{\theta })\Vert$$

$$\le V{U}_1+U_1^2{U}_2.$$

$$\mathbf{A}=\mathbf{X}+\mathbf{B}{\mathbf{B}}^T,$$

$${\underline{\mathbf{r}}}_n:={\underline{\mathbf{\theta }}}_{n+1}-\mathbf{A}{\underline{\mathbf{\theta }}}_n=-{\mathbf{M}}_{n+1}{\mathbf{G}}_{n+1}.$$

$${\underline{\mathbf{\theta }}}_{n+1}$$

$$={\mathbf{A}}^{n+1}{\underline{\mathbf{\theta }}}_0+\sum _{j=0}^n{\mathbf{A}}^{n-j}{\underline{\mathbf{r}}}_j.$$

$$\sum _{i=1}^N\sum _{k\in \mathcal{N}}(\mathcal{L}({\mathbf{x}}_{k,i},y_{k,i},{\mathbf{\psi }}_{k,i})-\mathcal{L}({\mathbf{x}}_{k,i},y_{k,i},\mathbf{g}))\le \gamma \sqrt{N}+\delta ,$$

$$\begin{array}{cc}[l]{\mathbf{\theta }}_{k,n}=& (1-\frac{1}{n}){\mathbf{\psi }}_{k,n-1}\\ & +{\mathbf{I}}_{+}(1-y_n{\mathbf{\psi }}_{k,n-1}^T{\mathbf{z}}_{\Omega }({\mathbf{x}}_{k,n}))\frac{y_{k,n}}{\lambda n}{\mathbf{z}}_{\Omega }({\mathbf{x}}_{k,n}),\end{array}$$

$${\mathbf{\theta }}_{k,n}={\mathbf{\psi }}_{k,n-1}+\mu {\epsilon }_{k,n}{\mathbf{z}}_{\Omega }({\mathbf{x}}_{k,n}),$$

$$y_{k,n}=\sum _{m=1}^M{a}_m\kappa ({\mathbf{c}}_m,{\mathbf{x}}_{k,n})+{\eta }_{k,n},$$

$$y_{k,n}$$

$$={\mathbf{a}}^T{Z}_{\Omega }^T{\mathbf{z}}_{\Omega }({\mathbf{x}}_{k,n})+{ϵ}_{k,n}+{\eta }_{k,n},$$

$$={\mathbf{\theta }}_o^T{\mathbf{z}}_{\Omega }({\mathbf{x}}_{k,n})+{ϵ}_{k,n}+{\eta }_{k,n},$$

$${\underline{\mathbf{y}}}_n={\mathbf{V}}_n^T{\underline{\mathbf{\theta }}}_o+{\underline{\mathbf{ϵ}}}_n+{\underline{\mathbf{\eta }}}_n,$$

$${\underline{\mathbf{\theta }}}_{\ast }=E[\mathbf{V}{\mathbf{V}}^T{]}^{-1}E[\mathbf{V}\underline{\mathbf{y}}],$$

$$E[\mathbf{V}\underline{\mathbf{y}}]$$

$$=E[\mathbf{V}{\mathbf{V}}^T]{\underline{\mathbf{\theta }}}_o+E[\mathbf{V}\underline{\mathbf{ϵ}}],$$

$${\underline{\mathbf{\theta }}}_{\ast }$$

$$={\underline{\mathbf{\theta }}}_o+E[\mathbf{V}{\mathbf{V}}^T{]}^{-1}E[\mathbf{V}\underline{\mathbf{ϵ}}]\approx {\underline{\mathbf{\theta }}}_o.$$

$$\underset{n}{\lim }E[{\mathbf{\theta }}_{k,n}-{\mathbf{\theta }}_o]=0_D,\text{ for all }k=1,2,\dots ,K.$$

$$r_{i,j}=$$

$$+\frac{1}{2}\exp \left(\frac{-\Vert {\mathbf{\omega }}_i+{\mathbf{\omega }}_j{\Vert }^2{\sigma }_X^2}{2}\right)\cos (b_i+b_j).$$

Full text

I Introduction

THE topic of distributed learning, has grown rapidly over the last years. This is mainly due to the exponentially increasing volume of data, that leads, in turn, to increased requirements for memory and computational resources. Typical applications include sensor networks, social networks, imaging, databases, medical platforms, e.t.c., [1]. In most of those, the data cannot be processed on a single processing unit (due to memory and/or computational power constrains) and the respective learning/inference problem has to be split into subproblems. Hence, one has to resort to distributed algorithms, which operate on data that are not available on a single location but are instead spread out over multiple locations, e.g., [2, 3, 4].

In this paper, we focus on the topic of distributed online learning and in particular to non linear parameter estimation and classification tasks. More specifically, we consider a decentralized network which comprises of nodes, that observe data generated by a non linear model in a sequential fashion. Each node communicates its own estimates of the unknown parameters to its neighbors and exploits simultaneously a) the information that it receives and b) the observed datum, at each time instant, in order to update the associated with it estimates. Furthermore, no assumptions are made regarding the presence of a central node, which could perform all the necessary operations. Thus, the nodes act as independent learners and perform the computations by themselves. Finally, the task of interest is considered to be common across the nodes and, thus, cooperation among each other is meaningful and beneficial, [5, 6].

The problem of linear online estimation has been considered in several works. These include diffusion-based algorithms, e.g., [7, 8, 9], ADMM-based schemes, e.g., [10, 11], as well as consencus-based ones, e.g., [12, 13]. The multitask learning problem, in which there are more than one parameter vectors to be estimated, has also been treated, e.g., [14, 15]. The literature on online distributed classification is more limited; in [16], a batch distributed SVM algorithm is presented, whereas in [17], a diffusion based scheme suitable for classification is proposed. In the latter, the authors study the problem of distributed online learning focusing on strongly-convex risk functions, such as the logistic regression loss, which is suitable to tackle classification tasks. The nodes of the network cooperate via the diffusion rationale. In contrast to the vast majority of works on the topic of distributed online learning, which assume a linear relationship between input and output measurements, in this paper we tackle the more general problem, i.e., the distributed online non–linear learning task. To be more specific, we assume that the data are generated by a model $y=f(\bm{x})$, where $f$ is a non-linear function that lies in a Reproducing Kernel Hilbert Space (RKHS). These are inner-product function spaces, generated by a specific kernel function, that have become popular models for non-linear tasks, since the introduction of the celebrated Support Vectors Machines (SVM) [18, 19, 20, 6].

Although there have been methods that attempt to generalize linear online distributed strategies to the non-linear domain using RKHS, mainly in the context of the Kernel LMS e.g., [21, 22, 23], these have major drawbacks. In [21] and [23], the estimation of $f$, at each node, is given as an increasingly growing sum of kernel functions centered at the observed data. Thus, a) each node has to transmit the entire sum at each time instant to its neighbors and b) the node has to fuse together all sums received by its neighbors to compute the new estimation. Hence, both the communications load of the entire network as well as the computational burden at each node grow linearly with time. Clearly, this is impractical for real life applications. In contrast, the method of [22] assumes that these growing sums are limited by a sparsification strategy; how this can be achieved is left for the future. Moreover, the aforementioned methods offer no theoretical results regarding the consensus of the network. In this work, we present a complete protocol for distributed online non-linear learning for both regression and classification tasks, overcoming the aforementioned problems. Moreover, we present theoretical results regarding network-wise consensus and regret bounds. The proposed framework offers fixed-size communication and computational load as time evolves. This is achieved through an efficient approximation of the growing sum using the random Fourier features rationale [24]. To the best of our knowledge, this is the first time that such a method appears in the literature.

Section II presents a brief background on kernel online methods and summarizes the main tools and notions used in this manuscript. The main contributions of the paper are presented in section III. The proposed method, the related theoretical results and extensive experiments can be found there. Section IV presents a special case of the proposed framework for the case of a single node. In this case, we demonstrate how the proposed scheme can be seen as a fixed-budget alternative for online kernel based learning (solving the problem of the growing sum). Finally, section V offers some concluding remarks. In the rest of the paper, boldface symbols denote vectors, while capital letters are reserved for matrices. The symbol $\otimes$ denotes the Kronecker product of matrices and the symbol $\cdot^{T}$ the transpose of the respective matrix or vector. Finally, the symbol $\|\cdot\|$ refers to the respective $\ell_{2}$ matrix or vector norm.

II Preliminaries

II-A RKHS

Reproducing Kernel Hilbert Spaces (RKHS) are inner product spaces of functions defined on $X$, whose respective point evaluation functional, i.e., $T_{x}:\mathcal{H}\rightarrow X:T_{x}(f)=f(x)$, is linear and continuous for every $x\in X$. This is usually portrayed by the reproducing property [18, 6, 25], which links inner products in $\mathcal{H}$ with a specific (semi-)positive definite kernel function $\kappa$ defined on $X\times X$ (associated with the space $\mathcal{H}$). As $\kappa(\cdot,x)$ lies in $\mathcal{H}$ for all $x\in X$, the reproducing property declares that $\langle\kappa(\cdot,y),\kappa(\cdot,x)\rangle_{\mathcal{H}}=\kappa(x,y)$, for all $x,y\in X$. Hence, linear tasks, defined on the high dimensional space, $\mathcal{H}$, (whose dimensionality can also be infinite) can be equivalently viewed as non-linear ones on the, usually, much lower dimensional space, $X$, and vice versa. This is the essence of the so called kernel trick: Any kernel-based learning method can be seen as a two step procedure, where firstly the original data are transformed from $X$ to $\mathcal{H}$, via an implicit map, $\Phi(x)=\kappa(\cdot,x)$, and then linear algorithms are applied to the transformed data. There exist a plethora of different kernels to choose from in the respective literature. In this paper, we mostly focus on the popular Gaussian kernel, i.e., $\kappa(\bm{x},\bm{y})=e^{\|\bm{x}-\bm{y}\|^{2}/(2\sigma^{2})}$, although any other shift invariant kernel can be adopted too.

Another important feature of RKHS is that any regularized ridge regression task, defined on $\mathcal{H}$, has a unique solution, which can written in terms of a finite expansion of kernel functions centered at the training points. Specifically, given the set of training points $\{(x_{n},y_{n}),n=1,\dots,N,x_{n}\in X,y_{n}\in\mathbb{R}\}$, the representer theorem [26, 18], states that the unique minimizer, $f_{*}\in\mathcal{H}$, of $\sum_{n=1}^{N}l(f(x_{n}),y_{n})+\lambda\|f\|^{2}_{\mathcal{H}}$, admits a representation of the form $f_{*}=\sum_{n=1}^{N}a_{n}\kappa(\cdot,x_{n})$, where $l$ is any convex loss function that measures the error between the actual system’s outputs, $y_{n}$, and the estimated ones, $f(x_{n})$, and $\|\cdot\|_{\mathcal{H}}$ is the norm induced by the inner product.

II-B Kernel Online Learning

The aforementioned properties have rendered RKHS a popular tool for addressing non linear tasks both in batch and online settings. Besides the widely adopted application on SVMs, in recent years there has been an increased interest on non linear online tasks around the squared error loss function. Hence, there have been kernel-based implementations of LMS [27, 28], RLS [29, 30], APSM [31, 32] and other related methods [33], as well as online implementations of SVMs [34], focusing on the primal formulation of the task. Henceforth, we will consider online learning tasks based on the training sequences of the form $\mathcal{D}=\{(\bm{x}_{n},y_{n}),\;n=1,2,\dots\}$, where $\bm{x}_{n}\in\mathbb{R}^{d}$ and $y_{n}\in\mathbb{R}$. The goal of the assumed learning tasks is to learn a non-linear input-output dependence, $y=f(\bm{x})$, $f\in\mathcal{H}$, so that to minimize a preselected cost. Note that these types of tasks include both classification (where $y_{n}=\pm 1$) and regression problems (where $y_{n}\in\mathbb{R}$). Moreover, in the online setting, the data are assumed to arrive sequentially.

As a typical example of these tasks, we consider the KLMS, which is one of the simplest and most representative methods of this kind. Its goal is to learn $f$, so that to minimize the MSE, i.e., $\mathcal{L}(f)=E[(y-f(\bm{x}))^{2}]$. Computing the gradient of $\mathcal{L}$ and estimating it via the current set of observations (in line with the stochastic approximation rationale, e.g., [6]), the estimate at the next iteration, employing the gradient descent method, becomes $f_{n}=f_{n-1}+\mu\epsilon_{n}\kappa(\bm{x}_{n},\cdot)$, where $\epsilon_{n}=y_{n}-f_{n-1}(\bm{x}_{n})$ and $\mu$ is the step-size (see, e.g., [6, 35, 36] for more). Assuming that the initial estimate is zero, the solution after $n-1$ steps turns out to be

[TABLE]

where $\alpha_{i}=\mu\epsilon_{i}$. Observe that this is in line with the representer theorem. Similarly, the system’s output can be estimated as $f_{n-1}(\bm{x}_{n})=\sum_{i=1}^{n-1}\alpha_{i}\kappa(\bm{x}_{n},\bm{x}_{i})$. Clearly, this linear expansion grows indefinitely as $n$ increases; hence the original form of KLMS is impractical. Typically, a sparsification strategy is adopted to bound the size of the expansion [37, 38, 39]. In these methods, a specific criterion is employed to decide whether a particular point, $\bm{x}_{n}$, is to be included to the expansion, or (if that point is discarded) how its respective output $y_{n}$ can be exploited to update the remaining weights of the expansion. There are also methods that can remove specific points from the expansion, if their information becomes obsolete, in order to increase the tracking ability of the algorithm [40].

II-C Approximating the Kernel with random Fourier Features

Usually, kernel-based learning methods involve a large number of kernel evaluations between training samples. In the batch mode of operation, for example, this means that a large kernel matrix has to be computed, increasing the computational cost of the method significantly. Hence, to alleviate the computational burden, one common approach is to use some sort of approximation of the kernel evaluation. The most popular techniques of this category are the Nyström method [41, 42] and the random Fourier features approach [24, 43]; the latter fits naturally to the online setting. Instead of relying on the implicit lifting, $\Phi$, provided by the kernel trick, Rahimi and Recht in [24] proposed to map the input data to a finite-dimensional Euclidean space (with dimension lower than $\mathcal{H}$ but larger than the input space) using a randomized feature map $\bm{z}_{\Omega}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D}$, so that the kernel evaluations can be approximated as $\kappa(\bm{x}_{n},\bm{x}_{m})\approx\bm{z}_{\Omega}(\bm{x}_{n})^{T}\bm{z}_{\Omega}(\bm{x}_{m})$. The following theorem plays a key role in this procedure.

Theorem 1**.**

Consider a shift-invariant positive definite kernel $\kappa(\bm{x}-\bm{y})$ defined on $\mathbb{R}^{d}$ and its Fourier transform $p(\bm{\omega})=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\kappa(\bm{\delta})e^{-i\bm{\omega}^{T}\bm{\delta}}d\bm{\delta}$, which (according to Bochner’s theorem) it can be regarded as a probability density function. Then, defining $z_{\bm{\omega},b}(\bm{x})=\sqrt{2}\cos(\bm{\omega}^{T}\bm{x}+b)$, it turns out that

[TABLE]

where $\bm{\omega}$ is drawn from $p$ and $b$ from the uniform distribution on $[0,2\pi]$.

Following Theorem 1, we choose to approximate $\kappa(\bm{x}_{n}-\bm{x}_{m})$ using $D$ random Fourier features, $\bm{\omega}_{1},\bm{\omega}_{2},\dots,\bm{\omega}_{D}$, (drawn from $p$) and $D$ random numbers, $b_{1},b_{2},\dots,b_{D}$ (drawn uniformly from $[0,2\pi]$) that define a sample average:

[TABLE]

Evidently, the larger $D$ is, the better this approximation becomes (up to a certain point). Details on the quality of this approximation can be found in [24, 43, 44, 45]. We note that for the Gaussian kernel, which is employed throughout the paper, the respective Fourier transform is

[TABLE]

which is actually the multivariate Gaussian distribution with mean $\bm{0}_{D}$ and covariance matrix $\frac{1}{\sigma^{2}}\bm{I}_{D}$.

We will demonstrate how this method can be applied using the KLMS paradigm. To this end, we define the map $\bm{z}_{\Omega}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D}$ as follows:

[TABLE]

where $\Omega$ is the $(d+1)\times D$ matrix defining the random Fourier features of the respective kernel, i.e.,

[TABLE]

provided that $\bm{\omega}$’s and $b$’s are drawn as described in theorem 1. Employing this notation, it is straightforward to see that (3) can be recast as $\kappa(\bm{x}_{n}-\bm{x}_{m})\approx\bm{z}_{\Omega}(\bm{x}_{m})^{T}\bm{z}_{\Omega}(\bm{x}_{n})$. Hence, the output associated with observation $\bm{x}_{n}$ can be approximated as

[TABLE]

It is a matter of elementary algebra to see that (7) can be equivalently derived by approximating the system’s output as $f(\bm{x})\approx\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x})$, initializing $\bm{\theta}_{0}$ to $\bm{0}_{D}$ and iteratively applying the following gradient descent type update: $\bm{\theta}_{n}=\bm{\theta}_{n-1}+\mu e_{n}\bm{z}_{\Omega}(\bm{x}_{n})$.

Clearly, the procedure described here, for the case of the KLMS, can be applied to any other gradient-type kernel based method. It has the advantage of modeling the solution as a fixed size vector, instead of a growing sum, a property that is quite helpful in distributed environments, as it will be discussed in section III.

III Distributed kernel-based Learning

In this section, we discuss the problem of online learning in RKHS over distributed networks. Specifically, we consider $K$ connected nodes, labeled $k\in{\mathcal{N}}=\{1,2,\dots K\}$, which operate in cooperation with their neighbors to solve a specific task. Let ${\mathcal{N}}_{k}\subseteq{\mathcal{N}}$ denote the neighbors of node $k$. The network topology is represented as an undirected connected graph, consisting of $K$ vertices (representing the nodes) and a set of edges connecting the nodes to each other (i.e., each node is connected to its neighbors). We assign a nonnegative weight $a_{k,l}$ to the edge connecting node $k$ to $l$. This weight is used by $k$ to scale the data transmitted from $l$ and vice versa. This can be interpreted as a measure of the confidence level that node $k$ assigns to its interaction with node $l$. We collect all coefficients into a $K\times K$ symmetric matrix $A=(a_{k,l})$, such that the entries of the $k$-th row of $A$ contain the coefficients used by node $k$ to scale the data arriving from its neighbors. We make the additional assumption that $A$ is doubly stochastic, so that the weights of all incoming and outgoing “transmissions” sum to 1. A common choice, among others, for choosing these coefficients, is the Metropolis rule, in which the weights equal to:

[TABLE]

Finally, we assume that each node, $k$, receives streaming data $\{(\bm{x}_{k,n},y_{k,n}),n=1,2,\dots\}$, that are generated from an input-output relationship of the form $y_{k,n}=f(\bm{x}_{k,n})+\eta_{k,n}$, where $\bm{x}_{k,n}\in\mathbb{R}^{d}$, $y_{k,n}$ belongs to $\mathbb{R}$ and $\eta_{k,n}$ represents the respective noise, for the regression task. The goal is to obtain an estimate of $f$. For classification, $y_{n,k}=\phi(f(\bm{x}_{k,n}))$, where, $\phi$ is a thresholding function; here we assume that $y_{n,k}\in\{-1,1\}$. Once more, the goal is to optimally estimate the classifier function $f$.

Each one of the nodes aims to estimate $f\in\mathcal{H}$ by minimizing a specific convex cost function, $\mathcal{L}(\bm{x},y,f)$, using a (sub)gradient descent approach. We employ a simple Combine-Then-Adapt (CTA) rationale, where at each time instant, $n$, each node, $k$, a) receives the current estimates, $f_{l,n-1}$, from all neighbors (i.e., from all nodes $l\in{\mathcal{N}}_{k}$), b) combines them to a single solution, $\psi_{k,n-1}=\sum_{l\in{\mathcal{N}}_{k}}a_{k,l}f_{l,n-1}$ and c) apply a step update procedure:

[TABLE]

The implementation of such an approach in the context of RKHS presents significant challenges. Keep in mind that, the estimation of the solution at each node is not a simple vector, but instead it is a function, which is expressed as a growing sum of kernel evaluations centered at the points observed by the specific node, as in (1). Hence, the implementation of a straightforward CTA strategy would require from each node to transmit its entire growing sum (i.e., the coefficients $a_{i}$ as well as the respective centers $\bm{x}_{i}$) to all neighbors. This would significantly increase both the communication load among the nodes, as well as the computational cost at each node, since the size of each one of the expansions would become increasingly larger as time evolves (as for every time instant, they gather the centers transmitted by all neighbors). This is the rationale adopted in [21, 22, 23] for the case of KLMS. Clearly, this is far from a practical approach. Alternatively, one could devise an efficient method to sparsify the solution at each node and then merge the sums transmitted by its neighbors. This would require (for example) to search all the dictionaries, transmitted by the neighboring nodes, for similar centers and treat them as a single one, or adopting a single pre-arranged dictionary (i.e., a specific set of centers) for all nodes and then fuse each observed point with the best-suited center. However, no such strategy has appeared in the respective literature, perhaps due to its increased complexity and lack of a theoretical elegance.

In this paper, inspired by the random Fourier features approximation technique, we approximate the desired input-output relationship as $y=\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x})$ and propose a two step procedure: a) we map each observed point $(\bm{x}_{k,n},y_{k,n})$ to $(\bm{z}_{\Omega}(\bm{x}_{k,n}),y_{k,n})$ and then b) we adopt a simple linear CTA diffusion strategy on the transformed points. Note that in the proposed scheme, each node aims to estimate a vector $\bm{\theta}\in\mathbb{R}^{D}$ by minimizing a specific (convex) cost function, $\mathcal{L}(\bm{x},y,\bm{\theta})$. Here, we imply that the model can be closely approximated by $y_{k,n}\approx\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x}_{k,n})+\eta_{k,n}$, for regression, and $y_{k,n}\sim\phi(\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x}_{k,n}))$ for classification, for all $k,n$, for some $\bm{\theta}$. We emphasize that $\mathcal{L}$ need not be differentiable. Hence, a large family of loss functions can be adopted. For example:

• •

Squared error loss: $\mathcal{L}(\bm{x},y,\bm{\theta})=(y-\bm{\theta}^{T}\bm{x})^{2}$.

• •

Hinge loss: $\mathcal{L}(\bm{x},y,\bm{\theta})=\max(0,1-y\bm{\theta}^{T}\bm{x})$.

We end up with the following generic update rule:

[TABLE]

where $\nabla_{\bm{\theta}}\mathcal{L}(\bm{z}_{\Omega}(\bm{x}_{k,n}),y_{k,n},\bm{\psi}_{k,n})$ is the gradient, or any subgradient of $\mathcal{L}(\bm{x},y,\bm{\theta})$ (with respect to $\bm{\theta}$), if the loss function is not differentiable. Algorithm 1 summarizes the aforementioned procedure. The advantage of the proposed scheme is that each node transmits a single vector (i.e., its current estimate, $\bm{\theta}_{k,n}$) to its neighbors, while the merging of the solutions requires only a straightforward summation.

III-A Consensus and regret bound

In the sequel, we will show that, under certain assumptions, the proposed scheme achieves asymptotic consensus and that the corresponding regret bound grows sublinearly with the time. It can readily be seen that (8)-(9) can be written more compactly (for the whole network) as follows:

[TABLE]

where $\underline{\bm{\theta}}_{n}:=(\bm{\theta}_{1,n}^{T},\ldots,\bm{\theta}_{K,n}^{T})^{T}\in\mathbb{R}^{KD}$, $\bm{M}_{n}:=\mathrm{diag}\{\mu_{1,n},\ldots,\mu_{K,n}\}\otimes I_{D}$, $\bm{G}_{n}:=[(\bm{u}_{1,n}^{T},\ldots,\bm{u}_{K,n}^{T}]^{T}\in\mathbb{R}^{KD}$, where $\bm{u}_{k,n}=\nabla\mathcal{L}(\bm{z}_{\Omega}(\bm{x}_{k,n}),y_{k,n},\bm{\psi}_{k,n})$, and $\bm{A}:=A\otimes I_{D}$. The necessary assumptions are the following:

Assumption 1**.**

The step size is time decaying and is bounded by the inverse square root of time, i.e., $\mu_{k,n}=\mu_{n}\leq\mu n^{-1/2}$.

Assumption 2**.**

The norm of the transformed input is bounded, i.e., $\exists U_{1}$ such that $\|\bm{z}_{\Omega}(\bm{x}_{k,n})\|\leq U_{1},\ \forall k\in{\mathcal{N}},\forall n\in\mathbb{N}$. Furthermore, $y_{k,n}$ is bounded, i.e., $|y_{k,n}|\leq V\ \forall k\in{\mathcal{N}},\forall n\in\mathbb{N}$ for some $V>0$.

Assumption 3**.**

The estimates are bounded, i.e., $\exists U_{2}$ s.t. $\|\bm{\theta}_{k,n}\|\leq U_{2},\ \forall k\in{\mathcal{N}},\forall n\in\mathbb{N}$.

Assumption 4**.**

The matrix comprising the combination weights, i.e., $A$, is doubly stochastic (if the weights are chosen with respect to the Metropolis rule, this condition is met).

Note that assumptions 2 and 3 are valid for most of the popular cost functions. As an example, we can study the squared error loss, i.e., $\mathcal{L}(\bm{x},y,\bm{\theta})=1/2(y-\bm{\theta}^{T}\bm{x})^{2}$, where:

[TABLE]

Following similar arguments, we can also prove that many other popular cost functions (e.g., the hinge loss, the logistic loss, e.t.c.) have bounded gradients too.

Proposition 1** (Asymptotic Consensus).**

All nodes converge to the same solution.

Proof.

Consider a $KD\times KD$ consensus matrix $\bm{A}$ as in (10). As $\bm{A}$ is doubly stochastic, we have the following [9]:

• •

$\|\bm{A}\|=1$.

• •

Any consensus matrix $\bm{A}$ can be decomposed as

[TABLE]

where $\bm{B}=[\bm{b}_{1},\ldots,\bm{b}_{D}]$ is an $KD\times D$ matrix, and $\bm{b}_{k}=1/\sqrt{K}(\bm{1}\otimes\bm{e}_{k})$, where $\bm{e}_{k}$, $k=1,\dots,D$ represent the standard basis of $\mathbb{R}^{D}$ and $\bm{X}$ is a $KD\times KD$ matrix for which it holds that $\left\|\bm{X}\right\|<1$.

• •

$\bm{A}\underline{\breve{\bm{\theta}}}=\underline{\breve{\bm{\theta}}}$, for all $\underline{\breve{\bm{\theta}}}\in\mathcal{O}:=\{\underline{\bm{\theta}}\in\mathbb{R}^{KD}:\ \underline{{\bm{\theta}}}=[\bm{\theta}^{T},\ldots,\bm{\theta}^{T}]^{T},\ \bm{\theta}\in\mathbb{R}^{D}\}$. The subspace $\mathcal{O}$ is the so called consensus subspace of dimension $D$, and $\bm{b}_{k},\ k=1,\ldots,D$, constitute a basis for this space. Hence, the orthogonal projection of a vector, $\underline{\bm{\theta}}$, onto this linear subspace is given by $P_{\mathcal{O}}(\underline{\bm{\theta}}):=\bm{B}\bm{B}^{T}\underline{\bm{\theta}}$, for all $\underline{\bm{\theta}}\in\mathbb{R}^{KD}$.

In [9], it has been proved that, the algorithmic scheme achieves asymptotic consensus, i.e., $\|\bm{\theta}_{k,n}-\bm{\theta}_{l,n}\|\rightarrow 0$, as $n\rightarrow\infty,$ for all $k,l\in{\mathcal{N}}$, if and only if $\lim_{n\rightarrow\infty}\|\bm{\underline{\bm{\theta}}}_{n}-P_{\mathcal{O}}(\bm{\underline{\bm{\theta}}}_{n})\|={0}$. We can easily check that the quantity

[TABLE]

approaches [math], as $n\rightarrow\infty$, since $\lim_{n\to\infty}\bm{M}_{n}=\bm{O}_{KD}$ (assumption 1) and the matrix $\bm{G}_{n}$ is bounded for all $n$. Rearranging the terms of (12) and iterating over $n$, we have:

[TABLE]

If we left-multiply the previous equation by $(\bm{I}_{KD}-\bm{B}\bm{B}^{T})$ and follow similar steps as in [9, Lemma 2], it can be verified that $\underset{n\rightarrow\infty}{\lim}\|\left(\bm{I}_{Km}-\bm{B}\bm{B}^{T}\right)\bm{\underline{\theta}}_{n+1}\|=0$, which completes our proof. ∎

Proposition 2**.**

Under assumptions 1-4 (and a cost function with bounded gradients) the networkwise regret is bounded by

[TABLE]

for all $\bm{g}\in\mathcal{B}_{[\bm{0}_{D},U_{2}]}$, where $\gamma,\ \delta$ are positive constants and $\mathcal{B}_{[\bm{0}_{D},U_{2}]}$ is the closed ball with center $\bm{0}_{D}$ and radius $U_{2}$.

Proof.

See appendix A. ∎

Remark 1**.**

It is worth pointing out that the theoretical properties, which were stated before, are complementary. In particular, the consensus property (Proposition 1) indicates that the nodes converge to the same solution and the sublinearity of the regret implies that on average the algorithm performs as well as the best fixed strategy. In fact, without the regret related proof we cannot characterize the solution in which the nodes converge.

III-B Diffusion SVM (Pegasos) Algorithm

The case of the regularized hinge loss function, i.e., $\mathcal{L}(\bm{x},y,\bm{\theta})=\frac{\lambda}{2}\|\bm{\theta}\|^{2}+\max\{0,1-y\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x})\}$, for a specific value of the regularization parameter $\lambda>0$, generates the Distributed Pegasos (see [34]). Note that the Pegasos solves the SVM task in the primal domain. In this case, the gradient becomes $\nabla_{\bm{\theta}}\mathcal{L}(\bm{x},y,\bm{\theta})=\lambda\bm{\theta}-\bm{I}_{+}(1-y\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x}))y\bm{z}_{\Omega}(\bm{x})$, where $\bm{I}_{+}$ is the indicator function of $(0,+\infty)$, which takes a value of 1, if its argument belongs in $(0,+\infty)$, and zero otherwise. Hence the step-update equation of algorithm 1 becomes:

[TABLE]

where, following [34], we have used a decreasing step size, $\mu_{n}=\frac{1}{\lambda n}$. This scheme satisfies the required assumptions, hence consensus is guaranteed.

We have tested the performance of Distributed-Pegasos versus the non-cooperative Pegasus on four datasets downloaded from Leon Bottou’s LASVM web page [46]. The chosen datasets are: a) the Adult dataset, b) the Banana dataset (where we have used the first 4000 points as training data and the remaining 1300 as testing data), c) the Waveform dataset (where we have used the first 4000 points as training data and the remaining 1000 as testing data) and d) the MNIST dataset (for the task of classifying the digit 8 versus the rest). The sizes of the datasets are given in Table I. In all experiments, we generate random graphs (using MIT’s random_graph routine, see [47]) and compare the proposed diffusion method versus a noncooperative strategy (where each node works independent of the rest). For each realization of the experiments, a different random connected graph with $M=5$ or $M=20$ nodes was generated, with probability of attachment per node equal to 0.2 (i.e, there is a 20% probability that a specific node $k$ is connected to any other node $l$). The adjacency matrix, $A$, of each graph was generated using the Metropolis rule. For the non-cooperative strategies, we used a graph that connects each node to itself, i.e., $A=I_{5}$ or $A=I_{20}$ respectively. The latter, implies that no information is exchanged between the nodes, thus each node is working alone. Moreover, for each realization, the corresponding dataset was randomly split into $M$ subsets of equal size (one for every node).

We note that the value of $D$ affects significantly the quality of the approximation via the Fourier features rationale and thus it also affects the performance of the experiments. The value of $D$ must be large enough so that the approximation is good, but not too large so that to the communicational and computational load become affordable. In practice, we can find a value for $D$ so that any further increase results to almost negligible performance variation (see also section IV). All other parameters were optimized (after trials) to give the lowest number of test errors. Their values are reported on Table IV. The algorithms were implemented in MatLab and the experiments were performed on a i7-3770 machine running at $3.4$GHz with 32 Mb of RAM. Tables II and III report the mean test errors obtained by both procedures. For $M=5$, the mean algebraic complexity of the generated graphs lies between $0.61$ and $0.76$ (different for each experiment), while the corresponding mean algebraic degree lies around 1.8. For $M=20$, the mean algebraic complexity of the generated graphs lies around $0.70$, while the corresponding mean algebraic degree lies around 3.9. The number inside the parentheses indicates the times of data reuse (i.e., running the algorithm again over the same data, albeit with a continuously decreasing step-size $\mu_{n}$), which has been suggested that improves the classification accuracy of Pegasos (see [34]). For example, the number 2 indicates that the algorithm runs over a dataset of double size, that contains the same data pairs twice. For the three first datasets (Adult, Banana, Waveform) we have run 100 realizations of the experiment, while for the fourth (MNIST) we have run only 10 (to save time). Besides the ADULT dataset, all other simulations show that the distributed implementation significantly outperforms the non-cooperative one. For that particular dataset, we observe that for a single run the non-cooperative strategy behaves better (for $M=20$), but as data reuse increases the distributed implementation reaches lower error floors.

III-C Diffusion KLMS

Adopting the squared error in place of $\mathcal{L}$, i.e., $\mathcal{L}(\bm{x},y,\bm{\theta})=(y-\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x}))^{2}$, and estimating the gradient by its current measurement, we take the Random Fourier Features Diffusion KLMS (RFF-DKLMS) and the step update becomes:

[TABLE]

where $\varepsilon_{k,n}=y_{n}-\bm{\psi}_{k,n-1}^{T}\bm{z}_{\Omega}(\bm{x}_{k,n})$. Although proposition 1 cannot be applied here (as it requires a decreasing step-size), we can derive sufficient conditions for consensus following the results of the standard Diffusion LMS [8]. Henceforth, we will assume that the data pairs are generated by

[TABLE]

where $\bm{c}_{1},\dots,\bm{c}_{M}$ are fixed centers, $\bm{x}_{k,n}$ are zero-mean i.i.d, samples drawn from the Gaussian distribution with covariance matrix $\sigma_{x}^{2}\bm{I}_{d}$ and $\eta_{k,n}$ are i.i.d. noise samples drawn from ${\mathcal{N}}(0,\sigma_{\eta}^{2})$. Following the RFF approximation rationale (for shift invariant kernels), we can write that

[TABLE]

where $Z_{\Omega}=\left(\bm{z}_{\Omega}(\bm{c}_{1}),\dots,\bm{z}_{\Omega}(\bm{c}_{M})\right)$, $\bm{a}=(a_{1},\dots,a_{M})^{T}$, $\bm{\theta}_{o}=Z_{\Omega}\bm{a}$ and $\epsilon_{k,n}$ is the approximation error between the noise-free component of $y_{k,n}$ (evaluated only by the linear kernel expansion of (15)) and the approximation of this component using random Fourier features, i.e., $\epsilon_{k,n}=\sum_{m=1}^{M}a_{m}\kappa(\bm{c}_{m},\bm{x}_{k,n})-\bm{\theta}_{o}^{T}\bm{z}_{\Omega}(\bm{x}_{k,n})$. For the whole network we have the following

[TABLE]

where

• •

$\underline{\bm{y}}_{n}:=(y_{1,n},y_{2,n},\dots,y_{K,n})^{T}$,

• •

$\bm{V}_{n}:=\mathop{\operator@font diag}(\bm{z}_{\Omega}(\bm{x}_{1,n}),\bm{z}_{\Omega}(\bm{x}_{2,n}),\dots,\bm{z}_{\Omega}(\bm{x}_{K,n}))$, is a $DK\times K$ matrix,

• •

$\underline{\bm{\theta}}_{o}=\left(\bm{\theta}_{o}^{T},\bm{\theta}_{o}^{T},\dots,\bm{\theta}_{o}^{T}\right)^{T}\in\mathbb{R}^{DK}$,

• •

$\underline{\bm{\epsilon}}_{n}=\left(\epsilon_{1,n},\epsilon_{2,n},\dots,\epsilon_{K,n}\right)^{T}\in\mathbb{R}^{K}$,

• •

$\underline{\bm{\eta}}_{n}=\left(\eta_{1,n},\eta_{2,n},\dots,\eta_{K,n}\right)^{T}\in\mathbb{R}^{K}$.

Let $\bm{x}_{1},\dots,\bm{x}_{K}\in\mathbb{R}^{d}$, $\underline{\bm{y}}\in\mathbb{R}^{K}$, be the random variables that generate the measurements of the nodes; it is straightforward to prove that the corresponding Wiener solution, i.e., $\underline{\bm{\theta}}_{*}=\textrm{argmin}_{\underline{\bm{\theta}}}E[\|\underline{\bm{y}}-\bm{V}^{T}\underline{\bm{\theta}}\|^{2}]$, becomes

[TABLE]

provided that the autocorrelation matrix $\bm{R}=E[\bm{V}\bm{V}^{T}]$ is invertible, where $\bm{V}=\mathop{\operator@font diag}(\bm{z}_{\Omega}(\bm{x}_{1}),\bm{z}_{\Omega}(\bm{x}_{2}),\dots,\bm{z}_{\Omega}(\bm{x}_{K}))$ is a $DK\times K$ matrix that collects the transformed random variables for the whole network. Assuming that the input-output relationship of the measurements at each node follows (16), the cross-correlation vector takes the form

[TABLE]

where for the last relation we have used that $\bm{\eta}$ is a zero mean vector representing noise and that $\bm{V}$ and $\underline{\bm{\eta}}$ are independent. For large enough $D$, the approximation error vector $\underline{\bm{\epsilon}}$ approaches $\bm{0}_{K}$, hence the optimal solution becomes:

[TABLE]

Here we actually imply that (16) can be closely approximated by $\bm{y}_{n}\approx\bm{V}_{n}\underline{\bm{\theta}}_{o}+\underline{\bm{\eta}}_{n}$; hence, the RFF-DKLMS is actually the standard diffusion LMS applied on the data pairs $\{(\bm{z}_{\Omega}(\bm{x}_{k,n}),y_{k,n}),\;k=1,\dots,K,\ n=1,2\dots\}$. The difference is that the input vectors $\bm{z}_{\Omega}(\bm{x}_{k,n})$ may have non zero mean and do not follow, necessarily, the Gaussian distribution. Hence, the available results regarding convergence and stability of diffusion LMS (e.g., [48, 49]) cannot be applied directly (in these works the inputs are assumed to be zero mean Gaussian to simplify the formulas related to stability). To this end, we will follow a slightly different approach. Regarding the autocorrelation matrix, we have the following result:

Lemma 1**.**

Consider a selection of samples $\bm{\omega}_{1},\bm{\omega}_{2},\dots,\bm{\omega}_{D}$, drawn from (4) such that $\bm{\omega}_{i}\not=\bm{\omega}_{j}$, for any $i\not=j$. Then, the matrix $\bm{R}=E[\bm{V}\bm{V}^{T}]$ is strictly positive definite (hence invertible).

Proof.

Observe that the $DK\times DK$ autocorrelation matrix is given by $\bm{R}=E[\bm{V}\bm{V}^{T}]=\mathop{\operator@font diag}(R_{zz},R_{zz}\dots,R_{zz})$, where $R_{zz}=E[\bm{z}_{\Omega}(\bm{x}_{k})\bm{z}_{\Omega}(\bm{x}_{k})^{T}]$, for all $k=1,2,\dots,K$. It suffices to prove that the $D\times D$ matrix $R_{zz}$ is strictly positive definite. Evidently, $\bm{c}^{T}R_{zz}\bm{c}=\bm{c}^{T}E\left[\bm{z}_{\Omega}(\bm{x}_{k})\bm{z}_{\Omega}(\bm{x}_{k})^{T}\right]\bm{c}=E\left[\left(\bm{z}_{\Omega}(\bm{x}_{k})^{T}\bm{c}\right)^{2}\right]\geq 0$, for all $\bm{c}\in\mathbb{R}^{D}$. Now, assume that there is a $\bm{c}\in\mathbb{R}^{D}$ such that $E\left[\left(\bm{z}_{\Omega}(\bm{x}_{k})^{T}\bm{c}\right)^{2}\right]=0$. Then $\bm{z}_{\Omega}(\bm{x})^{T}\bm{c}=0$ for all $\bm{x}\in\mathbb{R}^{D}$, or equivalently, $\sum_{i=1}^{D}c_{i}\cos(\bm{\omega}_{i}^{T}\bm{x}+b_{i})=0$, for all $\bm{x}\in\mathbb{R}^{D}$. Thus, $\bm{c}=\bm{0}$. ∎

As expected, the eigenvalues of $R_{zz}$ play a pivotal role in the convergence’s study of the algorithm. As $R_{zz}$ is a strictly positive definite matrix, its eigenvalues satisfy $0<\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{D}$.

Proposition 3**.**

If the the step update $\mu$ satisfies: $0<\mu<\frac{2}{\lambda_{D}}$, where $\lambda_{D}$ is the maximum eigenvalue of $R_{zz}$, then the RFF-DKLMS achieves asymptotic consensus in the mean, i.e.,

[TABLE]

Proof.

See Appendix B. ∎

Remark 2**.**

If $\bm{x}_{k,n}\sim{\mathcal{N}}(\bm{0},\sigma_{X}\bm{I}_{d})$, it is possible to evaluate explicitly the entries of $R_{zz}$, i.e.,

[TABLE]

Proposition 4**.**

For stability in the mean-square sense, we must ensure that both $\mu$ and $A$ satisfy:

[TABLE]

where $\boxtimes$ denotes the unbalanced block Kronecker product.

Proof.

See Appendix C. ∎

In the following, we present some experiments to illustrate the performance of the proposed scheme. We demonstrate that the estimation provided by the cooperative strategy is better than having each node working alone (i.e., lower MSE). Similar to section III-C, each realization of the experiments uses a different random connected graph with $M=20$ nodes and probability of attachment per node equal to 0.2. The adjacency matrix, $A$, of each graph was generated using the Metropolis rule (resulting to graphs with mean algebraic connectivity around $0.69$), while for the non-cooperative strategies, we used a graph that connects each node to itself, i.e., $A=I_{20}$. All parameters were optimized (after trials) to give the lowest MSE. The algorithms were implemented in MatLab and the experiments were performed on a i7-3770 machine running at $3.4$GHz with 32 Mb of RAM.

III-C1 Example 1. A Linear Expansion in terms of kernels

In this set-up, we generate $5000$ data pairs for each node using the following model: $y_{k,n}=\sum_{m=1}^{M}a_{m}\kappa(\bm{c}_{m},\bm{x}_{k,n})+\eta_{k,n}$, where $\bm{x}_{k,n}\in\mathbb{R}^{5}$ are drawn from ${\mathcal{N}}(\bm{0},I_{5})$ and the noise are i.i.d. Gaussian samples with $\sigma_{\eta}=0.1$. The parameters of the expansion (i.e., $a_{1},\dots,a_{M}$) are drawn from ${\mathcal{N}}(0,25)$, the kernel parameter $\sigma$ is set to $5$, the step update to $\mu=1$ and the number of random Fourier features to $D=2500$. Figure 1(a) shows the evolution of the MSE over all network nodes for 100 realizations of the experiment. We note that the selected value of step size satisfies the conditions of proposition 3.

III-C2 Example 2

Next, we generate the data pairs for each node using the following simple non-linear model: $y_{k,n}=\bm{w}_{0}^{T}\bm{x}_{k,n}+0.1\cdot(\bm{w}_{1}^{T}\bm{x}_{k,n})^{2}+\eta_{k,n}$, where $\eta_{k,n}$ represent zero-mean i.i.d. Gaussian noise with $\sigma_{\eta}=0.05$ and the coefficients of the vectors $\bm{w}_{0},\bm{w}_{1}\in\mathbb{R}^{5}$ are i.i.d. samples drawn from ${\mathcal{N}}(0,1)$. Similarly to Example 1, the kernel parameter $\sigma$ is set to $5$ and the step update to $\mu=1$. The number of random Fourier coefficients for RFFKLMS was set to $D=300$. Figure 3(b) shows the evolution of the MSE over all network nodes for 1000 realizations of the experiment over $15000$ samples.

III-C3 Example 3

Here we adopt the following chaotic series model [50]: $d_{k,n}=\frac{d_{k,n-1}}{1+d^{2}_{k,n-1}}+u^{3}_{k,n-1}$, $y_{k,n}=d_{k,n}+\eta_{k,n}$, where $\eta_{n}$ is zero-mean i.i.d. Gaussian noise with $\sigma_{\eta}=0.01$ and $u_{n}$ is also zero-mean i.i.d. Gaussian with $\sigma_{u}=0.15$. The kernel parameter $\sigma$ is set to $0.05$, the number of Fourier features to $D=100$ and the step update to $\mu=1$. We have also initialized $d_{1}$ to $1$. Figure 1(c) shows the evolution of the MSE over all network nodes for 1000 realizations of the experiment over 500 samples.

III-C4 Example 4

For the final example, we use another chaotic series model [50]: $d_{k,n}=u_{k,n}+0.5v_{k,n}-0.2d_{k,n-1}+0.35d_{k,n-2}$, $y_{k,n}=\phi(d_{k,n})+\eta_{k,n}$,

[TABLE]

where $\eta_{k,n}$, $v_{k,n}$ are zero-mean i.i.d. Gaussian noise with $\sigma_{\eta}=0.001$ and $\sigma_{v}^{2}=0.0156$ respectively, and $u_{k,n}=0.5v_{k,n}+\hat{\eta}_{k,n}$, where $\hat{\eta}_{n}$ is also i.i.d. Gaussian with $\sigma^{2}=0.0156$. The kernel parameter $\sigma$ is set to $0.05$ and the step update to $\mu=1$. We have also initialized $d_{1},d_{2}$ to $1$. Figure 3(d) shows the evolution of the MSE over all network nodes for $1000$ realizations of the experiment over $1000$ samples. The number of random Fourier features was set to $D=200$.

IV Revisiting Online Kernel based learning

In this section, we investigate the use of random Fourier features as a general framework for online kernel-based learning. The framework presented here can be seen as a special case of the general distributed method presented in section III for a network with a single node. Similar to the case of the standard KLMS, the learning algorithms considered here adopt a gradient descent rationale to minimize a specific loss function, $\mathcal{L}(\bm{x},y,f)$ for $f\in\mathcal{H}$, so that $f$ approximates the relationship between $\bm{x}$ and $y$, where $\mathcal{H}$ is the RKHS induced by a specific choice of a shift invariant (semi)positive definite kernel, $\kappa$. Hence, in general, these algorithms can be summarized by the following step update equation: $f_{n}=f_{n-1}+\mu_{n}\nabla_{f}\mathcal{L}(\bm{x}_{n},y_{n},f_{n-1})$. Algorithm 2 summarizes the proposed procedure for online kernel-based learning. The performance of the algorithm depends on the quality of the adopted approximation. Hence, a sufficiently large $D$ has to be selected.

Although algorithm 2 is given in a general setting, in the following we focus on the fixed-budget KLMS. As it has been discussed in section II, KLMS adopts the MSE cost function, which in the proposed framework takes the form: $\mathcal{L}(\bm{x},y,\bm{\theta})=E[(y_{n}-\bm{\theta}^{T}\bm{z}_{\Omega}(\bm{x}_{n}))^{2}]$. Hence, the respective step update equation of algorithm 2 becomes

[TABLE]

where $\varepsilon_{n}=y_{n}-\bm{\theta}_{n-1}^{T}\bm{z}_{\Omega}(\bm{x}_{n})$. Observe that, contrary to the typical implementations of KLMS, where the system’s output is a growing expansion of kernel functions and hence special care has to be carried out to prune the so called dictionary, the proposed approach employs a fixed-budget rationale, which doesn’t require any further treatment. We call this scheme the Random Fourier Features KLMS (RFF-KLMS) [51, 52].

The study of the convergence properties of RFFKLMS is based on those of the standard LMS. Henceforth, we will assume that the data pairs are generated by

[TABLE]

where $\bm{c}_{1},\dots,\bm{c}_{M}$ are fixed centers, $\bm{x}_{n}$ are zero-mean i.i.d, samples drawn from the Gaussian distribution with covariance matrix $\sigma_{x}^{2}\bm{I}_{d}$ and $\eta_{n}$ are i.i.d. noise samples drawn from ${\mathcal{N}}(0,\sigma_{\eta}^{2})$. Similar to the diffusion case, the eigenvalues of $R_{zz}$, i.e., $0<\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{D}$, play a pivotal role in the convergence’s study of the algorithm. Applying similar assumptions as in the case of the standard LMS (e.g., independence between $\bm{x}_{n},\bm{x}_{m}$, for $n\neq m$ and between $\bm{x}_{n},\eta_{n}$), we can prove the following results.

Proposition 5**.**

For datasets generated by (19) we have:

If $0<\mu<2/\lambda_{D}$, then RFFKLMS converges in the mean, i.e., $E[\bm{\theta}_{n}-\bm{\theta}_{\textrm{o}}]\rightarrow\bm{0}$. 2. 2.

The optimal MSE is given by

[TABLE]

For large enough $D$, we have $J_{n}^{\textrm{opt}}\approx\sigma_{\eta}^{2}$. 3. 3.

The excess MSE is given by $J_{n}^{\textrm{ex}}=J_{n}-J_{n}^{\textrm{opt}}=\mathop{\operator@font tr}\left(R_{zz}A_{n}\right)$, where $A_{n}=E[(\bm{\theta}_{n}-\bm{\theta}_{\textrm{o}})(\bm{\theta}_{n}-\bm{\theta}_{\textrm{o}})^{T}]$. 4. 4.

If $0<\mu<1/\lambda_{D}$, then $A_{n}$ converges. For large enough $n$ and $D$ we can approximate $A_{n}$’s evolution as $A_{n+1}\approx A_{n}-\mu\left(R_{zz}A_{n}+A_{n}R_{zz}\right)+\mu^{2}\sigma_{\eta}^{2}R_{zz}$. Using this model we can approximate the steady-state MSE ($\approx\mathop{\operator@font tr}\left(R_{zz}A_{n}\right)+\sigma_{\eta}^{2}$).

Proof.

The proofs use standard arguments as in the case of the standard LMS. Hence we do not provide full details. The reader is addressed to any LMS textbook.

1. See Proposition 3.

2. Replacing $\bm{\theta}_{n}$ with $\bm{\theta}_{\textrm{o}}$ in $J_{n}=E[\varepsilon_{n}^{2}]$ gives the result. For large enough $D$, $\epsilon_{n}$ is almost zero, hence we have $J_{n}^{\textrm{opt}}\approx\sigma_{\eta}^{2}$.

3. Here, we use the additional assumptions that $\bm{v}_{n}$ is independent of $\bm{x}_{n}$ and that $\epsilon_{n}$ is independent of $\eta_{n}$. The result follows after replacing $J_{n}$ and $J_{n}^{\textrm{opt}}$ and performing simple algebra calculations.

4. Replacing $\bm{\theta}_{\textrm{o}}$ and dropping out the terms that contain the term $\epsilon_{n}$, the result is obtained. ∎

Remark 3**.**

Observe that, while the first two results can be regarded as special cases of the distributed case (see proposition 3 and the related discussion in section III), the two last ones describe more accurately the evolution of the solution in terms of mean square stability, than the one given in proposition 4, for the general distributed scheme (where no formula for $\bm{B}_{n}$ is given). This becomes possible because the related formulas take a much simpler form, if the graph structure is reduced to a single node.

In order to illustrate the performance of the proposed algorithm and compare its behavior to the other variants of KLMS, we also present some related simulations. We choose the QKLMS [39] as a reference, since this is one of the most effective and fast KLMS pruning methods. In all experiments, that are presented in this section (described below), we use the same kernel parameter, i.e., $\sigma$, for both RFFKLMS and QKLMS as well as the same step-update parameter $\mu$. The quantization parameter $q$ of the QKLMS controls the size of the dictionary. If this is too large, then the dictionary will be small and the achieved MSE at steady state will be large. Typically, however, there is a value for $q$ for which the best possible MSE (which is very close to the MSE of the unsparsified version) is attained at steady state, while any smaller quantization sizes provide negligible improvements (albeit at significantly increased complexity). In all experimental set-ups, we tuned $q$ (using multiple trials) so that it leads to the best performance. On the other hand, the performance of RFFKLMS depends largely on $D$, which controls the quality of the kernel approximation. Similar to the case of QKLMS, there is a value for $D$ so that RFFKLMS attains its lowest steady-state MSE, while larger values provide negligible improvements. For our experiments, the chosen values for $q$ and $D$ provide results so that to trace out the results provided by the original (unsparsified) KLMS. Table V gives the mean training times for QKLMS and RFFKLMS on the same i7-3770 machine using a MatLab implementation (both algorithms were optimized for speed). We note that the complexity of the RFFKLMS is $\mathcal{O}(Dd)$, while the complexity of QKLMS is $\mathcal{O}(Md)$. Our experiments indicate that in order to obtain similar error floors, the required complexity of RFFKLMS is lower than that of QKLMS.

IV-1 Example 5. A Linear expansion in terms of Kernels

Similar to example 1 in section III-C, we generate $5000$ data pairs using (19) and the same parameters (for only one node). Figure 2 shows the evolution of the MSE for $500$ realizations of the experiment over different values of $D$. The algorithm reaches steady-state around $n=3000$. The attained MSE is getting closer to the approximation given in proposition 5 (dashed line in the figure) as $D$ increases. Figure 3(a) compares the performances of RFFKLMS and QKLMS for this particular set-up for $500$ realizations of the experiment using $8000$ data pairs. The quantization size of QKLMS was set to $q=5$ and the number of Fourier features for the RFFKLMS was set to $D=2500$.

IV-2 Example 6

Next, we use the same non-linear model as in example 2 of section III, i.e., $y_{n}=\bm{w}_{0}^{T}\bm{x}_{n}+0.1\cdot(\bm{w}_{1}^{T}\bm{x}_{n})^{2}+\eta_{n}$. The parameters of the model and the RFF-KLMS are the same as in example 1. The quantization size of the QKLMS was set to $q=5$, leading to an average dictionary size $M=100$. Figure 3(b) shows the evolution of the MSE for both QKLMS and RFFKLMS running $1000$ realizations of the experiment over $15000$ samples.

IV-3 Example 7

Here we adopt the same chaotic series model as in example 3 of section III-C, with the same parameters. Figure 3(c) shows the evolution of the MSE for both QKLMS and RFFKLMS running $1000$ realizations of the experiment over $500$ samples. The quantization parameter $q$ for the QKLMS was set to $q=0.01$, leading to an average dictionary size $M=7$.

IV-4 Example 8

For the final example, we use the chaotic series model of example 4 in section III-C with the same parameters. Figure 3(d) shows the evolution of the MSE for both QKLMS and RFFKLMS running $1000$ realizations of the experiment over $1000$ samples. The parameter $q$ was set to $q=0.01$, leading to $M=32$.

V Conclusion

We have presented a complete fixed-budget framework for non-linear online distributed learning in the context of RKHS. The proposed scheme achieves asymptotic consensus under some reasonable assumptions. Furthermore, we showed that the respective regret bound grows sublinearly with time. In the case of a network comprising only one node, the proposed method can be regarded as a fixed budget alternative for online kernel-based learning. The presented simulations validate the theoretical results and demonstrate the effectiveness of the proposed scheme.

Appendix A Proof of Proposition 2

In the following, we will use the notation $\mathcal{L}_{k,n}(\bm{\theta}):=\mathcal{L}(\bm{x}_{k,n},y_{k,n},\bm{\theta})$ to shorten the respective equations. Choose any $\bm{g}\in\mathcal{B}_{[\bm{0}_{D},U_{2}]}$. It holds that

[TABLE]

Moreover, as $\mathcal{L}_{k,n}$ is convex, we have:

[TABLE]

for all $\bm{\theta},\bm{\theta}^{\prime}\in\mathrm{dom}(\mathcal{L}_{k,n})$ where $\bm{h}:=\nabla\mathcal{L}_{k,n}(\bm{\theta})$ is the gradient (for a differentiable cost function) or a subgradient (for the case of a non–differentiable cost function). From (20), (21) and the boundness of the (sub)gradient we take

[TABLE]

where $U$ is an upper bound for the (sub)gradient. Recall that for the whole network we have: $\bm{\underline{\psi}}_{n}=\bm{A}\bm{\underline{\theta}}_{n-1}$ and that for any doubly stochastic matrix, $\bm{A}$, its norm equals to its largest eigenvalue, i.e., $\|\bm{A}\|=\lambda_{\textrm{max}}=1$. A respective eigenvector is $\bm{\underline{g}}=(\bm{g}^{T}\dots,\bm{g}^{T})^{T}\in\mathbb{R}^{DK}$, hence it holds that $\bm{\underline{g}}=\bm{A}\bm{\underline{g}}$ and

[TABLE]

where $\underline{\bm{\psi}}_{n}=(\bm{\psi}_{n}^{T},\dots,\bm{\psi}_{n}^{T})^{T}\in\mathbb{R}^{DK}$. Going back to (22) and summing over all $k\in{\mathcal{N}}$, we have:

[TABLE]

However, for the left hand side of the inequality we obtain $\sum_{k\in{\mathcal{N}}}(\|\bm{\psi}_{k,n}-\bm{g}\|^{2}-\|\bm{\theta}_{k,n}-\bm{g}\|^{2})=\|\underline{\bm{\psi}}_{n}-\underline{\bm{g}}\|^{2}-\|\underline{\bm{\theta}}_{n}-\underline{\bm{g}}\|^{2}.$ If we combine the last relation with (23) and (24) we have

[TABLE]

The last inequality leads to

[TABLE]

where we have taken into consideration, Assumption 3 and the boundeness of $g$. Next, summing over $i=1,\ldots,N+1$, taking into consideration that $\sum_{i=1}^{N}\mu_{i}\leq 2\mu\sqrt{N}$ (Assumption 1) and noticing that some terms telescope, we have:

[TABLE]

Rearranging the terms and omitting the negative ones completes the proof:

[TABLE]

Appendix B Proof of Proposition 3

For the whole network, the step update of RFF-DKLMS can be recasted as

[TABLE]

where $\underline{\bm{\varepsilon}}_{n}=(\varepsilon_{1,n},\varepsilon_{2,n},\dots,\varepsilon_{K,n})^{T}$ and $\varepsilon_{k,n}=y_{k,n}-\bm{\psi}_{k,n}^{T}\bm{z}_{\Omega}(\bm{x}_{k,n})$, or equivalently, $\underline{\bm{\varepsilon}}_{n}=\underline{\bm{y}}_{n}-\bm{V}_{n}^{T}\bm{A}\underline{\bm{\theta}}_{n-1}$. If we define $\bm{U}_{n}=\underline{\bm{\theta}}_{n}-\underline{\bm{\theta}}_{o}$ and take into account that $\bm{A}\underline{\bm{g}}=\underline{\bm{g}}$, for all $\underline{\bm{g}}\in\mathbb{R}^{DK}$, such that $\underline{\bm{g}}=(\bm{g}^{T},\bm{g}^{T},\dots,\bm{g}^{T})^{T}$ for $\bm{g}\in\mathbb{R}^{D}$, we obtain:

[TABLE]

If we take the mean values and assume that $\bm{\theta}_{k,n}$ and $\bm{z}_{\Omega}(\bm{x}_{k,n})$ are independent for all $k=1,\dots,K$, $n=1,2,\dots,$ we have

[TABLE]

Taking into account that $\bm{\eta}_{n}$ and $\bm{V}_{n}$ are independent, that $E[\bm{\eta}_{n}]=\bm{0}$ and that for large enough $D$ we have $E[\bm{V}_{n}\bm{\epsilon}_{n}]\approx\bm{0}$, we can take $E[\bm{U}_{n}]\approx\left((I_{KD}-\mu\bm{R})\bm{A}\right)^{n-1}E[\bm{U}_{1}]$. Hence, if all the eigenvalues of $(I_{KD}-\mu\bm{R})\bm{A}$ have absolute value less than $1$, we have that $E[\bm{U}_{n}]\rightarrow\bm{0}$. However, since $\bm{A}$ is a doubly stochastic matrix we have $\|\bm{A}\|\leq 1$ and

[TABLE]

Moreover, as $I_{KD}-\mu\bm{R}$ is a diagonal block matrix, its eigenvalues are identical to the eigenvalues of its blocks, i.e., the eigenvalues of $I_{D}-\mu R_{zz}$. Hence, a sufficient condition for convergence is $|1-\mu\lambda_{D}(R_{zz})|<1$, which gives the result. $\square$

Remark 4**.**

Observe that $|\lambda_{\max}\left((I_{KD}-\mu\bm{R})\bm{A}\right)|\leq|\lambda_{\max}\left((I_{KD}-\mu\bm{R})I_{KD}\right)|$, which means that the spectral radius of $(I_{KD}-\mu\bm{R})\bm{A}$ is generally smaller than that of $(I_{KD}-\mu\bm{R})I_{KD}$ (which corresponds to the non-cooperative protocol). Hence, cooperation under the diffusion rationale has a stabilizing effect on the network [8].

Appendix C Proof of Proposition 4

Let $\bm{B}_{n}=E[\bm{U}_{n}\bm{U}_{n}^{T}]$, where $\bm{U}_{n}=\bm{A}\bm{U}_{n-1}-\mu\bm{V}_{n}\bm{V}_{n}^{T}\bm{A}\bm{U}_{n-1}+\mu\bm{V}_{n}\underline{\bm{\epsilon}}_{n}+\mu\bm{V}_{n}\underline{\bm{\eta}}_{n}$. Taking into account that the noise is i.i.d., independent from $\bm{U}_{n}$ and $\bm{V}_{n}$ and that $\bm{\epsilon}_{n}$ is close to zero (if $D$ is sufficiently large), we can take that:

[TABLE]

For sufficiently small step-sizes, the rightmost term can be neglected [53, 49], hence we can take the simplified form

[TABLE]

Next, we observe that $\bm{B}_{n}$, $\bm{R}$ and $\bm{A}$ can be regarded as block matrices, that consist of $K\times K$ blocks with size $D\times D$. We will vectorize equation (C) using the $\mathop{\operator@font vecb_{r}}$ operator, as this has been defined in [54]. Assuming a block-matrix $C$:

[TABLE]

the $\mathop{\operator@font vecb_{r}}$ operator applies the following vectorization:

[TABLE]

Moreover, it is closely related to the following block Kronecker product:

[TABLE]

The interested reader can delve into the details of the $\mathop{\operator@font vecb_{r}}$ operator and the unbalanced block Kronecker product in [54]. Here, we limit our interest to the following properties:

$\mathop{\operator@font vecb_{r}}(DCE^{T})=(E\boxtimes D)\mathop{\operator@font vecb_{r}}C$. 2. 2.

$(C\boxtimes D)(E\boxtimes F)=CE\boxtimes DF$.

Thus, applying the $\mathop{\operator@font vecb_{r}}$ operator, on both sizes of (C) we take $\bm{b}_{n}=(\bm{A}\boxtimes\bm{A})\bm{b}_{n-1}-\mu\left((\bm{RA})\boxtimes\bm{A}\right)\bm{b}_{n-1}-\mu\left((\bm{A})\boxtimes\bm{RA}\right)\bm{b}_{n-1}+\mu^{2}\sigma_{\eta}^{2}\bm{r}$, where $\bm{b}_{n}=\mathop{\operator@font vecb_{r}}\bm{B}_{n}$ and $\bm{r}=\mathop{\operator@font vecb_{r}}\bm{R}$. Exploiting the second property, we can take:

[TABLE]

Hence, we finally get:

[TABLE]

which gives the result.