Online Graph-Adaptive Learning with Scalability and Privacy

TL;DR

This paper introduces a scalable, privacy-preserving online learning method for estimating nodal attributes in large, dynamic graphs using a multikernel approach that enables real-time updates without batch re-computation.

Contribution

It proposes a novel multikernel-based online learning algorithm that efficiently updates nodal attribute estimates in large, evolving graphs while maintaining privacy through encrypted connectivity data.

Findings

Effective on synthetic datasets

Demonstrates scalability to large networks

Maintains privacy with encrypted connectivity

Abstract

Graphs are widely adopted for modeling complex systems, including financial, biological, and social networks. Nodes in networks usually entail attributes, such as the age or gender of users in a social network. However, real-world networks can have very large size, and nodal attributes can be unavailable to a number of nodes, e.g., due to privacy concerns. Moreover, new nodes can emerge over time, which can necessitate real-time evaluation of their nodal attributes. In this context, the present paper deals with scalable learning of nodal attributes by estimating a nodal function based on noisy observations at a subset of nodes. A multikernel-based approach is developed which is scalable to large-size networks. Unlike most existing methods that re-solve the function estimation problem over all existing nodes whenever a new node joins the network, the novel method is capable of providing real-time evaluation of the function values on newly-joining nodes without resorting to a batch solver. Interestingly, the novel scheme only relies on an encrypted version of each node's connectivity in order to learn the nodal attributes, which promotes privacy. Experiments on both synthetic and real datasets corroborate the effectiveness of the proposed methods.

Figures

Equations

$$\displaystyle{\mathbf{y}}={\mbox{\boldmath$\Psi$}}{\mathbf{x}}+{\mathbf{e}}$$

$$\mathcal{H}:=\{f:f(v)=\sum _{n=1}^N{\alpha }_n\kappa (v,v_n),{\alpha }_n\in \mathbb{R}\}$$

$$\hat{f}=\arg \underset{f\in \mathcal{H}}{\min }\text{leavevmode}\frac{1}{M}\sum _{m=1}^M\mathcal{C}(f(v_m),y_m)+\mu \Omega \left(\Vert f{\Vert }_{\mathcal{H}}^2\right)$$

$$\underset{\bar{\mathbf{\alpha }}\in {\mathbb{R}}^N}{\min }\frac{1}{M}\sum _{m=1}^M\mathcal{C}({\bar{\mathbf{\alpha }}}^{\top }\bar{\mathbf{k}}(v_m),y_m)+\mu \Omega \left({\bar{\mathbf{\alpha }}}^{\top }\bar{\mathbf{K}}\bar{\mathbf{\alpha }}\right).$$

$$\hat{f}(v)=\sum _{m=1}^M{\alpha }_m\kappa (v,v_m):={\mathbf{\alpha }}^{\top }\mathbf{k}(v).$$

$$\min_{{\color[rgb]{0,0,0}\bm{\alpha}\in\mathbb{R}^{M}}}\leavevmode\nobreak\ \frac{1}{M}\sum_{m=1}^{M}{\cal C}({\bm{\alpha}}^{\top}{\mathbf{k}}(v_{m}),y_{m})+\mu\Omega\left({\bm{\alpha}}^{\top}{\mathbf{K}}{\bm{\alpha}}\right)$$

$$\displaystyle\bar{{\mathbf{K}}}:=r^{\dagger}({\mathbf{L}}):={\mathbf{U}}r^{\dagger}({\mbox{\boldmath$\Lambda$}}){\mathbf{U}}^{\top}$$

$$[\bar{\mathbf{K}}{]}_{n,n^{\prime }}=\kappa (v_n,v_{n^{\prime }})=\kappa ({\mathbf{a}}_n,{\mathbf{a}}_{n^{\prime }}).$$

$$\hat{f}(v_n)=\hat{f}({\mathbf{a}}_n)=\sum _{m=1}^M{\alpha }_m\kappa ({\mathbf{a}}_m,{\mathbf{a}}_n):={\mathbf{\alpha }}^{\top }\mathbf{k}({\mathbf{a}}_n)$$

$$\kappa ({\mathbf{a}}_n-{\mathbf{a}}_{n^{\prime }})$$

$$\displaystyle:=\mathbb{E}_{{\mathbf{v}}}\big{[}e^{j{\mathbf{v}}^{\top}({\mathbf{a}}_{n}-{\mathbf{a}}_{n^{\prime}})}\big{]}$$

$$\hat{\kappa }({\mathbf{a}}_n,{\mathbf{a}}_{n^{\prime }})={\mathbf{z}}_{\mathbf{V}}^{\top }({\mathbf{a}}_n){\mathbf{z}}_{\mathbf{V}}({\mathbf{a}}_{n^{\prime }})$$

$${\mathbf{z}}_{\mathbf{V}}(\mathbf{a})$$

$$\times {\left[\sin ({\mathbf{v}}_1^{\top }\mathbf{a}),\dots ,\sin ({\mathbf{v}}_D^{\top }\mathbf{a}),\cos ({\mathbf{v}}_1^{\top }\mathbf{a}),\dots ,\cos ({\mathbf{v}}_D^{\top }\mathbf{a})\right]}^{\top }.$$

$$\displaystyle\hat{f}^{\rm RF}({\mathbf{a}})=\sum_{m=1}^{M}\alpha_{m}{\mathbf{z}}_{{\mathbf{V}}}^{\top}({\mathbf{a}}_{m}){\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}):={\mbox{\boldmath$\theta$}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}})$$

$$\min_{{\mbox{\boldmath$\theta$}}\in\mathbb{R}^{2D}}\leavevmode\nobreak\ \frac{1}{M}\sum_{m=1}^{M}{\cal C}({\mbox{\boldmath$\theta$}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{m}),y_{m})+\mu\Omega\left(\|{\mbox{\boldmath$\theta$}}\|^{2}\right)$$

$$\displaystyle\min_{\bm{\theta}\in\mathbb{R}^{2D}}\,\sum_{\tau=1}^{t}{\cal L}\left({\mbox{\boldmath$\theta$}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{\tau}),y_{\tau}\right)$$

$$\displaystyle{\cal L}\big{(}{\mbox{\boldmath$\theta$}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t}\big{)}:={\cal C}\big{(}{\mbox{\boldmath$\theta$}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t}\big{)}+\mu\Omega\big{(}\|{\mbox{\boldmath$\theta$}}\|^{2}\big{)}.$$

$$\displaystyle{\mbox{\boldmath$\theta$}}_{t+1}={\mbox{\boldmath$\theta$}}_{t}-\eta_{t}\nabla{\cal L}({\mbox{\boldmath$\theta$}}_{t}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t})$$

$$\displaystyle\hat{f}(v_{i})=\hat{{\mbox{\boldmath$\theta$}}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{i}),\leavevmode\nobreak\ \leavevmode\nobreak\ \forall i\in{\cal S}^{c}$$

$$\displaystyle\nabla{\cal L}\left({\mbox{\boldmath$\theta$}}_{t}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t}\right)=2[{\mbox{\boldmath$\theta$}}_{t}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t})-y_{t}]{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t})+2\mu{\mbox{\boldmath$\theta$}}_{t}.$$

$$\nabla \mathcal{L}$$

$$\displaystyle{\color[rgb]{0,0,0}\nabla{\cal L}\!\!\left({\mbox{\boldmath$\theta$}}_{t}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t}\right)\!=\!\frac{2y_{t}\exp(y_{t}{\mbox{\boldmath$\theta$}}_{t}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}))}{1+\exp(y_{t}{\mbox{\boldmath$\theta$}}_{t}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}))}{\mathbf{z}}_{p}({\mathbf{x}}_{t})+2\mu{\mbox{\boldmath$\theta$}}_{t}.}$$

$$f(v)=f(\mathbf{a}):=\sum _{p=1}^P{\bar{w}}_p{f}_p(\mathbf{a})\in \bar{\mathcal{H}}$$

$$\underset{\{{\bar{w}}_p\},\{{\hat{f}}_p^{\mathrm{RF}}\}}{\min }\sum _{t=1}^T\sum _{p=1}^P{\bar{w}}_p{\mathcal{L}}_t\left({\hat{f}}_p^{\mathrm{RF}}({\mathbf{a}}_t)\right)$$

$$s.\text{leavevmode}\mathrm{to}$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \hat{f}_{p}^{\rm RF}\!\in\!\left\{\hat{f}_{p}({\mathbf{a}}_{t})\!=\!{\mbox{\boldmath$\theta$}}^{\top}{\mathbf{z}}_{{\mathbf{V}}_{p}}({\mathbf{a}}_{t})\right\},\leavevmode\nobreak\ p\in{\cal P}$$

$$\displaystyle\hat{f}_{p,t}^{\rm RF}({\mathbf{a}}_{t})={\mbox{\boldmath$\theta$}}_{p,t}^{\top}{\mathbf{z}}_{p}({\mathbf{a}}_{t})$$

$$\displaystyle{\mbox{\boldmath$\theta$}}_{p,t+1}={\mbox{\boldmath$\theta$}}_{p,t}-\eta\nabla{\cal L}({\mbox{\boldmath$\theta$}}_{p,t}^{\top}{\mathbf{z}}_{p}({\mathbf{a}}_{t}),y_{t})$$

$$w_{p,t+1}=w_{p,t}\exp \left(-\eta {\mathcal{L}}_t\left({\hat{f}}_{p,t}^{\mathrm{RF}}({\mathbf{a}}_t)\right)\right).$$

$$\displaystyle f(v_{n})=\sum_{i=1}^{I}\beta_{i}f_{i}({\mbox{\boldmath$\phi$}}_{i,n})$$

$${\mathrm{Reg}}_{\mathcal{A}}^{\mathrm{s}}(T):=\sum _{t=1}^T{\mathcal{L}}_t({\hat{f}}_t({\mathbf{a}}_t))-\sum _{t=1}^T{\mathcal{L}}_t(f^{\ast }({\mathbf{a}}_t))$$

$$f^{\ast }(\cdot )$$

Full text

000 Work in this paper was supported by grants NSF 1711471, 1500713 and NIH 1R01GM104975-01.000∗Y. Shen and G. B. Giannakis are with the Dept. of ECE and the Digital Technology Center, University of Minnesota, 117 Pleasant Str. SE, Minneapolis, MN 55455. Tel: (612)625-4287; Emails: {shenx513,georgios}@umn.edu 000*†* G. Leus is with the Electrical Engineering, Mathematics and Computer Science, Delft University of Technology, CD, Delft 2826, The Netherlands E-mail:, g.j.t.leus@tudelft.nl

I Introduction

Estimating nodal functions/signals over networks is a task emerging in various domains, such as social, brain, and power networks, to name a few. Functions of nodes can represent certain attributes or classes of these nodes. In Facebook for instance, each node represents a person, and the presence of an edge indicates that two persons are friends, while nodal attributes can be age, gender or movie ratings of each person. In financial networks, where each node is a company, with links denoting trade between two companies, the function of the node can represent the category that each company belongs to, e.g., technology-, fashion-, or education-related.

In real-world networks, there are often unavailable nodal function values, due to, e.g., privacy issues. Hence, a topic of great practical importance is to interpolate missing nodal values (class, ranking or function), based on the function values at a subset of observed nodes. Interpolation of nodal function values often relies on the assumption of “smoothness” over the graphs, which implies that neighboring nodes will have similar nodal function values. For example, in social networks, people tend to rate e.g., movies similar to their friends, and in financial networks, companies that trade with each other usually belong to the same category. From this point of view, function estimation over graphs based on partial observations has been investigated extensively, [1, 2, 3, 4, 5, 6]. Function estimation has been also pursued in the context of semi-supervised learning, e.g., for transductive regression or classification, see e.g., [7, 8, 9, 10]. The same task has been studied recently as signal reconstruction over graphs, see e.g., [11, 12, 13, 14, 15], where signal values on unobserved nodes can be estimated by properly introducing a graph-aware prior. Kernel-based methods for learning over graphs offer a unifying framework that includes linear and nonlinear function estimators [13, 16, 17]. The nonlinear methods outperform the linear ones but suffer from the curse of dimensionality [18], rendering them less attractive for large-scale networks.

To alleviate this limitation, a scalable kernel-based approach will be introduced in the present paper, which leverages the random feature approximation to ensure scalability while also allowing real-time evaluation of the functions over large-scale dynamic networks. In addition, the novel approach incorporates a data-driven scheme for adaptive kernel selection.

Adaptive learning over graphs has been also investigated for tracking and learning over possibly dynamic networks, e.g., [19, 17]. Least mean-squares and recursive least-squares adaptive schemes have been developed in [19], without explicitly accounting for evolving network topologies. In contrast, [17] proposed a kernel-based reconstruction scheme to track time-varying signals over time-evolving topologies, but assumed that the kernel function is selected a priori. All these prior works assume that the network size is fixed.

In certain applications however, new nodes may join the network over time. For example, hundreds of new users are joining Facebook or Netflix every day, and new companies are founded in financial networks regularly. Real-time and scalable estimation of the desired functions on these newly-joining nodes is of great importance. While simple schemes such as averaging over one- or multi-hop neighborhoods are scalable to network size by predicting the value on each newly-coming node as a weighted combination of its multi-hop neighborhoods [20], they do not capture global information over the network. In addition, existing rigorous approaches are in general less efficient in accounting for newly-joining nodes, and need to solve the problem over all nodes, every time new nodes join the network, which incurs complexity $\mathcal{O}(N^{3})$, where $N$ denotes the network size [13, 16]. As a result, these methods are not amenable to real-time evaluation over newly-joining nodes. To this end, the present paper develops a scalable online graph-adaptive algorithm that can efficiently estimate nodal functions on newly-joining nodes ‘on the fly.’

Besides scalability and adaptivity, nodes may have firm privacy requirements, and may therefore not be willing to reveal who their neighbors are. However, most graph-based learning methods require knowing the entire connectivity pattern, and thus cannot meet the privacy requirements. The novel random feature based approach on the other hand, only requires an encrypted version of each node’s connectivity pattern, which makes it appealing for networks with stringent privacy constraints.

In short, we put forth a novel online multikernel learning (MKL) framework for effectively learning and tracking nonlinear functions over graphs. Our contributions are as follows.

c1) A scalable MKL approach is developed to efficiently estimate the nodal function values both on the observed and un-observed nodes of a graph;

c2) The resultant algorithm is capable of estimating the function value of newly incoming nodes with high accuracy without having to solve the batch problem over all nodes, making it highly scalable as the network size grows, and suitable for nodal function estimation in dynamic networks;

c3) Unlike most existing methods that rely on nodal feature vectors in order to learn the function, the proposed scheme simply capitalizes on the connectivity pattern of each node, while at the same time, nodal feature vectors can be easily incorporated if available; and,

c4) The proposed algorithm does not require nodes to share connectivity patterns. Instead, a privacy-preserving scheme is developed for estimating the nodal function values based on an encrypted version of the nodal connectivity patterns, hence respecting node privacy.

The rest of this paper is organized as follows. Preliminaries are in Section II, while Section III presents an online kernel-based algorithm that allows sequential processing of nodal samples. Section IV develops an online MKL scheme for sequential data-driven kernel selection, which allows graph-adaptive selection of kernel functions to best fit the learning task of interest. Finally, results of corroborating numerical tests on both synthetic and real data are presented in Section VI, while concluding remarks along with a discussion of ongoing and future directions are given in Section VII.

Notation. Bold uppercase (lowercase) letters denote matrices (column vectors), while $(\cdot)^{\top}$ and $\lambda_{i}(.)$ stand for matrix transposition, and the $i$th leading eigenvalue of the matrix argument, respectively. The identity matrix will be represented by $\mathbf{I}$, while $\mathbf{0}$ will denote the matrix of all zeros, and their dimensions will be clear from the context. Finally, the $\ell_{p}$ and Frobenius norms will be denoted by $\|\cdot\|_{p}$, and $\|\cdot\|_{F}$, respectively.

II Kernel-based learning over graphs

Consider a graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ of $N$ nodes, whose topology is captured by a known adjacency matrix ${\mathbf{A}}\in\mathbb{R}^{N\times N}$. Let $a_{nn^{\prime}}\in\mathbb{R}$ denote the $(n,n^{\prime})$ entry of ${\mathbf{A}}$, which is nonzero only if an edge is present from node $n^{\prime}$ to $n$. A real-valued function (or signal) on a graph is a mapping $f:{\cal V}\rightarrow\mathbb{R}$, where ${\cal V}$ is the set of vertices. The value $f(v)=x_{v}$ represents an attribute of $v\in{\cal V}$, e.g., in the context of brain networks, $x_{v_{n}}$ could represent the sample of an electroencephalogram (EEG), or functional magnetic resonance imaging (fMRI) measurement at region $n$. In a social network, $x_{v_{n}}$ could denote the age, political alignment, or annual income of the $n$th person. Suppose that a collection of noisy samples $\{y_{m}=x_{v_{n_{m}}}+e_{m}\}_{m=1}^{M}$ is available, where $e_{m}$ models noise, and $M\leq N$ represents the number of measurements. Given $\{y_{m}\}_{m=1}^{M}$, and with the graph topology known, the goal is to estimate $f(v)$, and thus reconstruct the graph signal at unobserved vertices. Letting ${\mathbf{y}}:=[y_{1},\ldots,y_{M}]^{\top}$, the observation vector obeys

[TABLE]

where ${\mathbf{x}}:=[x_{v_{1}},\dots,x_{v_{N}}]^{\top}$, ${\mathbf{e}}:=[e_{1},\dots,e_{M}]^{\top}$, and {\mbox{\boldmath\Psi}}\in\{0,1\}^{M\times N} is a sampling matrix with binary entries [{\mbox{\boldmath\Psi}}]_{m,n_{m}}=1 for $m=1,\dots,M$, and [math], elsewhere.

Given $\Psi$, ${\mathbf{y}}$, and ${\mathbf{A}}$, the goal is to estimate ${\mathbf{x}}$ over the entire network. To tackle the under-determined system (1), consider function $f$ belonging to a reproducing kernel Hilbert space (RKHS) defined as [16, 13]

[TABLE]

where $\kappa:\mathcal{V}\times\mathcal{V}\rightarrow\mathbb{R}$ is a pre-selected kernel function. Hereafter, we will let $n_{m}=m$ for notational convenience, and without loss of generality (wlog). Given ${\mathbf{y}}$, the RKHS-based estimate is formed as

[TABLE]

where the cost ${\cal C}(\cdot,\cdot)$ can be selected depending on the learning task, e.g., the least-squares (LS) for regression, or the logistic loss for classification; $\|f\|_{\mathcal{H}}^{2}:=\sum_{n}\sum_{n^{\prime}}\alpha_{n}\alpha_{n^{\prime}}\kappa(v_{n},v_{n^{\prime}})$ is the RKHS norm; $\Omega(\cdot)$ is an increasing function; and, $\mu>0$ is a regularization parameter that copes with overfitting. According to the definition of graph RKHS in (2), the function estimate can be written as \hat{f}(v)=\sum_{n=1}^{N}\alpha_{n}\kappa(v,v_{n}):=\bar{{\mbox{\boldmath\alpha}}}^{\top}\bar{\mathbf{k}}(v), where $\bar{\bm{\alpha}}:=[\alpha_{1},\ldots,\alpha_{N}]^{\top}\!\in\mathbb{R}^{N}$ collects the basis coefficients, and $\bar{\mathbf{k}}({v}):=[\kappa(v,v_{1}),\ldots,\kappa(v,v_{N})]^{\top}\!$. Substituting into the RKHS norm, we find $\|f\|_{\mathcal{H}}^{2}:=\sum_{n}\sum_{n^{\prime}}\alpha_{n}\alpha_{n^{\prime}}\kappa(v_{n},v_{n^{\prime}})=\bar{\bm{\alpha}}^{\top}\bar{\mathbf{K}}\bar{\bm{\alpha}}$, where the $N\times N$ kernel matrix $\bar{\mathbf{K}}$ has entries $[\bar{\mathbf{K}}]_{n,n^{\prime}}:=\kappa(v_{n},v_{n^{\prime}})$; thus, the functional problem (3) boils down to

[TABLE]

According to the representer theorem, the optimal solution of (3) admits the finite-dimensional form given by [16, 13]

[TABLE]

where $\bm{\alpha}:=[\alpha_{1},\ldots,\alpha_{M}]^{\top}\!\in\mathbb{R}^{M}$, and $\mathbf{k}(v):=[\kappa(v,v_{1}),\ldots,\kappa(v,v_{M})]^{\top}\!$. This means that the coefficients corresponding to the unobserved nodes are all zeros. This implies that the function over the graph can be estimated by optimizing over the $M\times 1$ vector $\alpha$ [cf. (3)]

[TABLE]

where {\mathbf{K}}:={\mbox{\boldmath\Psi}}^{\top}\bar{\mathbf{K}}{\mbox{\boldmath\Psi}}. For general kernel-based learning tasks, $\bar{{\mathbf{K}}}$ is formed using the nonlinear functions of pairwise correlations \kappa(v_{n},v_{n^{\prime}})={\mbox{\boldmath\phi}}_{n}^{\top}{\mbox{\boldmath\phi}}_{n^{\prime}}, where {\mbox{\boldmath\phi}}_{n} denotes the feature vector of node $n$, which can collect, for example, the buying history of users on Amazon, or the trading history of companies in financial networks. However, such information may not be available in practice, due to, e.g., privacy concerns. This has motivated the graph-kernel based approaches in [13] and [16], to reconstruct the graph signal when only the network structure is available, and the kernel matrix is selected as a nonlinear function of the graph Laplacian matrix. Specifically, these works mostly consider undirected networks, ${\mathbf{A}}={\mathbf{A}}^{\top}$.

Given the normalized Laplacian matrix ${\mathbf{L}}:={\mathbf{I}}-{\mathbf{D}}^{-1/2}{\mathbf{A}}{\mathbf{D}}^{-1/2}$, with ${\mathbf{D}}:={\rm diag}({\mathbf{A}}\mathbf{1})$, and letting {\mathbf{L}}:={\mathbf{U}}{\mbox{\boldmath\Lambda}}{\mathbf{U}}^{\top}, the family of graphical kernels is

[TABLE]

where $r(.)$ is a non-decreasing scalar function of the eigenvalues, and † denotes pseudo-inverse. By selecting $r(.)$, different graph properties can be accounted for, including smoothness, band-limitedness, the random walk [16], and diffusion [2].

Although graph-kernel based methods are effective in reconstructing signals over graphs, it can be observed from (7) that formulating $\bar{{\mathbf{K}}}$ generally requires an eigenvalue decomposition of ${\mathbf{L}}$, which incurs complexity ${\cal O}(N^{3})$ that can be prohibitive for large-scale networks. Moreover, even though nodal feature vectors $\{{\mathbf{\phi}}_{n}\}$ are not necessary to form $\bar{{\mathbf{K}}}$, the graph-kernel-based scheme requires knowledge of the topology, meaning ${\mathbf{A}}$, in order to estimate the nodal function of each node. However, in networks with strict privacy requirements, nodes may not be willing to share such information with others. In Facebook, for example, most people do not make their friend list public. In addition, solving (4) assumes that all sampled nodes are available in batch, which may not be true in scenarios where nodes are sampled in a sequential fashion.

In response to these challenges, an online scalable kernel-based method will be developed in the ensuing section to deal with sequentially obtained data samples, over generally dynamic networks. The resultant algorithm only requires encrypted versions of the nodal connectivity patterns of other nodes, and hence it offers privacy.

III Online kernel-based learning over graphs

Instead of resorting to a graph kernel that requires an eigenvalue decomposition of ${\mathbf{L}}$ in (7), the present section advocates treating the connectivity pattern of each node as its feature vector, which can be the $n$th column ${\mathbf{a}}_{n}^{(c)}$ and possibly the $n$th row $({\mathbf{a}}_{n}^{(r)})^{\top}$ of the adjacency (if ${\mathbf{A}}$ is nonsymmetric). We will henceforth term this the connectivity pattern of $v_{n}$, and denote it as ${\mathbf{a}}_{n}$, for brevity. Given ${\mathbf{a}}_{n}$, we will interpolate unavailable nodal function values $\hat{f}(v_{n})$ using a nonparametric approach, that is different and scalable relative to [16] and [13]. The kernel matrix is now

[TABLE]

Again, with $M$ nodes sampled, the representer theorem asserts that the sought function estimator has the form [18]

[TABLE]

where ${\mathbf{k}}({\mathbf{a}}_{n}):=[\kappa({\mathbf{a}}_{n},{\mathbf{a}}_{1})\dots\kappa({\mathbf{a}}_{n},{\mathbf{a}}_{M})]^{\top}$. It can be observed from (9) that $\hat{f}(v_{n})$ involves the adjacency of the entire network, namely $\{{\mathbf{a}}_{m}\}_{m=1}^{M}$, which leads to potentially growing complexity ${\cal O}(M^{3})$ as the number of sampled nodes increases [18].

III-A Batch RF-based learning over graphs

To bypass this growing complexity, we will resort to the so-called random feature approximation [21] in order to reduce the original functional learning task in (4) to a problem with the number of unknown parameters not growing with $M$. We first approximate $\kappa$ in (5) using random features (RFs) [21, 22] that are obtained from a shift-invariant kernel satisfying $\kappa({\mathbf{a}}_{n},{\mathbf{a}}_{n^{\prime}})=\kappa({\mathbf{a}}_{n}-{\mathbf{a}}_{n^{\prime}})$. For $\kappa({\mathbf{a}}_{n}-{\mathbf{a}}_{n^{\prime}})$ absolutely integrable, its Fourier transform $\pi_{\kappa}(\bf v)$ exists and represents the power spectral density, which upon normalizing to ensure $\kappa(\mathbf{0})=1$, can also be viewed as a probability density function (pdf); hence,

[TABLE]

where the last equality is due to the definition of the expected value. Drawing a sufficient number of $D$ independent and identically distributed samples $\{{\mathbf{v}}_{i}\}_{i=1}^{D}$ from $\pi_{\kappa}({\mathbf{v}})$, the ensemble mean (III-A) can be approximated by the sample average

[TABLE]

where ${\mathbf{V}}:=[{\mathbf{v}}_{1},\dots,{\mathbf{v}}_{D}]^{\top}\in\mathbb{R}^{D\times N}$, and ${\mathbf{z}}_{{\mathbf{V}}}$ denotes the $2D\times 1$ real-valued RF vector

[TABLE]

Taking expectations in (11) and using (III-A), one can verify that $\mathbb{E}_{{\mathbf{v}}}[\hat{\kappa}({\mathbf{a}}_{n},{\mathbf{a}}_{n^{\prime}})]=\kappa({\mathbf{a}}_{n},{\mathbf{a}}_{n^{\prime}})$, which means ${\hat{\kappa}}$ is unbiased. Note that finding $\pi_{\kappa}({\mathbf{v}})$ requires an $N$-dimensional Fourier transform of $\kappa$, which in general requires numerical integration. Nevertheless, it has been shown that for a number of popular kernels, $\pi_{\kappa}({\mathbf{v}})$ is available in closed form [21]. Taking the Gaussian kernel as an example, where \kappa({\mathbf{a}}_{n},{\mathbf{a}}_{n^{\prime}})=\exp\big{(}\|{\mathbf{a}}_{n}-{\mathbf{a}}_{n^{\prime}}\|_{2}^{2}/(2\sigma^{2})\big{)}, it has a Fourier transform corresponding to the pdf $\mathcal{N}(0,\sigma^{-2}{\mathbf{I}})$.

Hence, the function that is optimal in the sense of (3) can be cast to a function approximant over the $2D$-dimensional RF space (cf. (9) and (11))

[TABLE]

where {\mbox{\boldmath\theta}}^{\top}:=\sum_{m=1}^{M}\alpha_{m}{\mathbf{z}}_{{\mathbf{V}}}^{\top}({\mathbf{a}}_{m}). While $\hat{f}$ in (5) is the superposition of nonlinear functions $\kappa$, its RF approximant $\hat{f}^{\rm RF}$ in (13) is a linear function of ${\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{i})$. As a result, (3) reduces to

[TABLE]

where \|{\mbox{\boldmath\theta}}\|^{2}:=\sum_{t}\sum_{\tau}\alpha_{t}\alpha_{\tau}{\mathbf{z}}_{{\mathbf{V}}}^{\top}({\mathbf{a}}_{t}){\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{\tau}):=\|f\|_{\cal H}^{2}. A batch solver of (14) has complexity $\mathcal{O}(MD^{3})$ that does not grow with $N$. This batch RF-based approach scales linearly with the number of measured nodes $M$, and the number of variables is $2D$, which does not depend on $M$. This allows us to pursue an online implementation as elaborated next.

III-B Online RF-based learning over graphs

Here, we will further leverage RF-based learning over graphs to enable real-time learning and reconstruction of signals evolving over possibly dynamic networks. A scalable online algorithm will be introduced, which can adaptively handle sequentially sampled nodal features and update the sought function estimates.

Training sequentially. In the training phase, we are given a network of $N$ nodes, and the nodal function is sampled in a sequential fashion. Letting $v_{t}$ denote the node sampled at the $t$th time slot, and having available $\{{\mathbf{a}}_{t},y_{t}\}$ at $v_{t}$, the online inference task can be written as [cf. (14)]

[TABLE]

We will solve (15) using online gradient descent [23]. Obtaining $v_{t}$ per slot $t$, the RF of its connectivity pattern ${\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t})$ is formed as in (12), and {\mbox{\boldmath\theta}}_{t+1} is updated ‘on the fly,’ as

[TABLE]

where $\{\eta_{t}\}$ is the sequence of stepsizes that can tune learning rates. In this paper, we will adopt $\eta_{t}=\eta$ for simplicity. Iteration (16) provides a functional update since \hat{f}^{\rm RF}_{t}({\mathbf{a}})={\mbox{\boldmath\theta}}_{t}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}). The per-iteration complexity of (16) is $\mathcal{O}(D)$, and $\mathcal{O}(MD)$ for the entire training process, which scales better than $\mathcal{O}(MD^{3})$ that is required for a batch solver of (14).

Inferring unavailable nodal values. After the training phase, the nodal function value over the un-sampled nodes can be readily estimated by [cf. (13)]

[TABLE]

where \hat{{\mbox{\boldmath\theta}}} is the final estimate after the training phase, i.e., \hat{{\mbox{\boldmath\theta}}}={\mbox{\boldmath\theta}}_{M+1}, and ${\cal S}^{c}$ denotes the index set of the nodes whose signal values have not been sampled in the training phase.

Newly-joining nodes. When new nodes join the network, batch graph-kernel based approaches must expand $\bar{{\mathbf{K}}}$ in (7) by one row and one column, and re-solve (6) in order to form signal estimates for the newly-joining nodes. Hence, each newly joining node will incur complexity $\mathcal{O}(N^{3})$. The novel online RF method on the other hand, can simply estimate the signal on the newly coming node via \hat{f}(v_{\rm new})=\hat{{\mbox{\boldmath\theta}}}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{\rm new}), where ${\mathbf{a}}_{\rm new}\in\mathbb{R}^{N}$ denotes the connectivity pattern of the new node with the existing nodes in the network. This leads to a complexity of $\mathcal{O}(ND)$ per new node. If in addition, $y_{\rm new}$ is available, then the function estimate can also be efficiently updated via (16) and (13) using ${\mathbf{a}}_{\rm new}$ and $y_{\rm new}$.

The steps of our online RF-based method are summarized in Algorithm 1. A couple of standard learning tasks where Algorithm 1 comes handy are now in order.

Nonlinear regression over graphs. Consider first nonlinear regression over graphs, where the goal is to find a nonlinear function $f\in{\cal H}$, such that $y_{n}=f(v_{n})+e_{n}=f({\mathbf{a}}_{n})+e_{n}$ given the graph adjacency matrix ${\mathbf{A}}$. The criterion is to minimize the regularized prediction error of $y_{n}$, typically using the online LS loss ${\cal L}(f({\mathbf{a}}_{t}),y_{t}):=[y_{t}-f({\mathbf{a}}_{t})]^{2}+\mu\|f\|_{\cal H}^{2}$ in (15), whose gradient is (cf. (21))

[TABLE]

In practice, $y_{t}$ can represent a noisy version of each node’s real-valued attribute, e.g., temperature in a certain city, and the graph can be constructed based on Euclidean distances among cities. For a fair comparison with alternatives, only the regression task will be tested in the numerical section of this paper.

Nonlinear classification over graphs. We can also handle kernel-based perceptron and kernel-based logistic regression, which aim at learning a nonlinear classifier that best approximates either $y_{n}$, or, the pdf of $y_{n}$ conditioned on ${\mathbf{a}}_{n}$. With binary labels $\{\pm 1\}$, the perceptron solves (3) with ${\cal L}(f({\mathbf{a}}_{t}),y_{t})=\max(0,1-y_{t}f({\mathbf{a}}_{t}))+\mu\|f\|_{\cal H}^{2}$, which equals zero if $y_{t}=f({\mathbf{a}}_{t})$, and otherwise equals $1$. In this case, the gradient of the presented online RF-based method is (cf. (21))

[TABLE]

Accordingly, given ${\bf x}_{t}$, logistic regression postulates that ${\rm Pr}(y_{t}=1|{\mathbf{x}}_{t})=1/(1+\exp(f({\mathbf{x}}_{t})))$. Here the gradient takes the form (cf. (21))

[TABLE]

Classification over graphs arises in various scenarios, where $y_{n}$ may represent categorical attributes such as gender, occupation or, nationality of users in a social network.

Remark 1 (Privacy). Note that the update in (16) does not require access to ${\mathbf{a}}_{t}$ directly. Instead, the only information each node needs to reveal is ${\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t})$ for each ${\mathbf{a}}_{t}$, which involves $\{\sin({\mathbf{a}}_{t}^{\top}{\mathbf{v}}_{j}),\leavevmode\nobreak\ \cos({\mathbf{a}}_{t}^{\top}{\mathbf{v}}_{j})\}_{j=1}^{D}$. Being noninvertible, these co-sinusoids functions involved in generating the $\mathbf{z}_{\mathbf{V}}(\mathbf{a}_{t})$ can be viewed as an encryption of the nodal connectivity pattern, which means that given $\mathbf{z}_{\mathbf{V}}(\mathbf{a}_{t})$, vector $\mathbf{a}_{t}$ cannot be uniquely deciphered. Hence, Algorithm 1 preserves privacy.

Remark 2 (Directed graphs). It can be observed from (7) that for $\bar{{\mathbf{K}}}$ to be a valid kernel, graph-kernel based methods require ${\mathbf{A}}$, and henceforth ${\mathbf{L}}$ to be symmetric, which implies they can only directly deal with symmetric/undirected graphs. Such a requirement is not necessary for our RF-based method.

**Remark 3 (Dynamic graphs). Real-world networks may vary over time, as edges may disappear or appear. To cope with such changing topologies, the original graph-kernel method needs to recalculate the kernel matrix, and resolve the batch problem whenever one edge changes. In contrast, our online RF-based method can simply re-estimate the nodal values on the two ends of the (dis)appeared edge using (13) with their current $\{{\mathbf{a}}_{n}\}$. **

Evidently, the performance of Algorithm 1 depends on $\kappa$ that is so far considered known. As the “best” performing $\kappa$ is generally unknown and application dependent, it is prudent to adaptively select kernels by superimposing multiple kernels from a prescribed dictionary, as we elaborate next.

IV Online Graph-adaptive MKL

In the present section, we develop an online graph-adaptive learning approach that relies on random features, and leverages multi-kernel approximation to estimate the desired $f$ based on sequentially obtained nodal samples over the graph. The proposed method is henceforth abbreviated as Gradraker.

The choice of $\kappa$ is critical for the performance of single kernel based learning over graphs, since different kernels capture different properties of the graph, and thus lead to function estimates of variable accuracy [13]. To deal with this, combinations of kernels from a preselected dictionary $\{\kappa_{p}\}_{p=1}^{P}$ can be employed in (3); see also [13, 22]. Each combination belongs to the convex hull $\bar{{\cal K}}:=\{\bar{\kappa}=\sum_{p=1}^{P}\bar{\alpha}_{p}\kappa_{p},\,\bar{\alpha}_{p}\geq 0,\,\sum_{p=1}^{P}\bar{\alpha}_{p}=1\}$. With $\bar{\cal H}$ denoting the RKHS induced by $\bar{\kappa}\in\bar{{\cal K}}$, one then solves (3) with ${\cal H}$ replaced by $\bar{\cal H}:={\cal H}_{1}\bigoplus\cdots\bigoplus{\cal H}_{P}$, where $\{{\cal H}_{p}\}_{p=1}^{P}$ represent the RKHSs corresponding to $\{\kappa_{p}\}_{p=1}^{P}$ [24].

The candidate function ${\bar{f}}\in\bar{\cal H}$ is expressible in a separable form as $\bar{f}({\mathbf{a}}):=\sum_{p=1}^{P}{\bar{f}}_{p}({\mathbf{a}})$, where ${\bar{f}}_{p}({\mathbf{a}})$ belongs to $\mathcal{H}_{p}$, for $p\in{\cal P}:=\{1,\ldots,P\}$. To add flexibility per kernel in our ensuing online MKL scheme, we let wlog $\{{\bar{f}}_{p}={w}_{p}f_{p}\}_{p=1}^{P}$, and seek functions of the form

[TABLE]

where $f:={\bar{f}}/\sum_{p=1}^{P}w_{p}$, and the normalized weights $\{\bar{w}_{p}:=w_{p}/\sum_{p=1}^{P}w_{p}\}_{p=1}^{P}$ satisfy $\bar{w}_{p}\geq 0$, and $\sum_{p=1}^{P}\bar{w}_{p}=1$. Exploiting separability jointly with the RF-based function approximation per kernel, the MKL task can be reformulated, after letting {\cal L}_{t}(\hat{f}_{p}^{\rm RF}({\mathbf{a}}_{t})):={\cal L}\big{(}{\mbox{\boldmath\theta}}^{\top}{\mathbf{z}}_{{\mathbf{V}}_{p}}({\mathbf{a}}_{t}),y_{t}\big{)} in (15), as

[TABLE]

which can be solved efficiently ‘on-the-fly.’ Relative to (14), we replaced $M$ by $T$ to introduce the notion of time, and stress the fact that the nodes are sampled sequentially.

Given the connectivity pattern ${\mathbf{a}}_{t}$ of the $t$th sampled node $v_{t}$, an RF vector ${\mathbf{z}}_{p}({\mathbf{a}}_{t})$ is generated per $p$ from the pdf $\pi_{\kappa_{p}}({\mathbf{v}})$ via (12), where ${\mathbf{z}}_{p}({\mathbf{a}}_{t}):={\mathbf{z}}_{{\mathbf{V}}_{p}}({\mathbf{a}}_{t})$ for notational brevity. Hence, per kernel $\kappa_{p}$ and node sample $t$, we have [cf. (13)]

[TABLE]

and as in (16), {\mbox{\boldmath\theta}}_{p,t} is updated via

[TABLE]

with $\eta\in(0,1)$ chosen constant to effect the adaptation. As far as $\bar{w}_{p,t}$ is concerned, since it resides on the probability simplex, a multiplicative update is well motivated as discussed also in, e.g., [23, 22]. For the un-normalized weights, this update is available in closed form as [22]

[TABLE]

Having found $\{w_{p,t}\}$ as in (22), the normalized weights in (18) are obtained as $\bar{w}_{p,t}:=w_{p,t}/\sum_{p=1}^{P}w_{p,t}$. Note from (22) that when $\hat{f}_{p,t}^{\rm RF}$ has a larger loss relative to other $\hat{f}_{p^{\prime},t}^{\rm RF}$ with $p^{\prime}\neq p$ for the $t$th sampled node, the corresponding $w_{p,t+1}$ decreases more than the other weights. In other words, a more accurate approximant tends to play a more important role in predicting the ensuing sampled node. In summary, our Gradraker for online graph MKL is listed as Algorithm 2.

Remark 4 (Comparison with batch MKL). A batch MKL based approach for signal reconstruction over graphs was developed in [13]. It entails an iterative algorithm whose complexity grows with $N$ in order to jointly estimate the nodal function, and to adaptively select the kernel function. When new nodes join the network, [13] re-calculates the graphical kernels and re-solves the overall batch problem, which does not scale with the network size. In addition, [13] is not privacy preserving in the sense that in order to estimate the function at any node, one needs to have access to the connectivity pattern of the entire network.

Remark 5 (Comparison with $k$-NN). An intuitive yet efficient way to predict function values of a newly joining node is to simply combine the values of its $k$ nearest neighbors ($k$-NN) [20, 25]. Efficient as it is, $k$-NN faces several challenges: a) At least one of the neighbors must be labeled, which does not always hold in practice, and is not required by the Gradraker; and b) $k$-NN can only account for local information, while the Gradraker takes also into account the global information of the graph.

IV-A Generalizations

So far, it is assumed that each node $n$ only has available its own connectivity feature vector ${\mathbf{a}}_{n}$. This allows Gradraker to be applied even when limited information is available about the nodes, which many existing algorithms that rely on nodal features cannot directly cope with.

If additional feature vectors $\{\bm{\phi}_{i,n}\}_{i=1}^{I}$ are actually available per node $n$ other than its own ${\mathbf{a}}_{n}$, it is often not known a priori which set of features is the most informative for estimating the signal of interest on the graph. To this end, the novel Gradraker can be adapted by treating the functions learned from different sets of features as an ensemble of learners, and combine them in a similar fashion as in (18), that is,

[TABLE]

Applications to several practical scenarios are discussed in the following.

Semi-private networks. In practice, a node may tolerate sharing its links to its neighbors, e.g., users of Facebook may share their friends-list with friends. In this scenario, each node not only knows its own neighbors, but also has access to who are its neighbors’ neighbors, i.e., two-hop neighbors. Specifically, node $n$ has access to ${\mathbf{a}}_{n}$, as well as to the $n$th column of ${\mathbf{A}}^{(2)}:={\mathbf{A}}{\mathbf{A}}$ [1], and a learner f_{2}({\mbox{\boldmath\phi}}_{2,n}) can henceforth be introduced and combined in (23). Moreover, when nodes are less strict about privacy, e.g., when a node is willing to share its multi-hop neighbors, more learners can be introduced and combined ‘on the fly’ by selecting {\mbox{\boldmath\phi}}_{i,n} as the $n$th column of ${\mathbf{A}}^{(i)}$ in (23).

Multilayer networks. Despite their popularity, ordinary networks are often inadequate to describe increasingly complex systems. For instance, modeling interactions between two individuals using a single edge can be a gross simplification of reality. Generalizing their single-layer counterparts, multilayer networks allow nodes to belong to $N_{g}$ groups, called layers [26, 27]. These layers could represent different attributes or characteristics of a complex system, such as temporal snapshots of the same network, or different types of groups in social networks (family, soccer club, or work related). Furthermore, multilayer networks are able to model systems that typically cannot be represented by traditional graphs, such as heterogeneous information networks [28, 29]. To this end, Gradraker can readily incorporate the information collected from heterogenous sources, e.g., connectivity patterns $\{{\mathbf{A}}_{i}\}_{i=1}^{N_{g}}$ from different layers, by adopting a kernel based learner $f_{i}({\mathbf{a}}_{i,n})$ on the $i$th layer and combining them as in (23).

Nodal features available. In certain cases, nodes may have nodal features [1] in addition to their $\{{\mathbf{a}}_{n}\}$. For example, in social networks, other than the users’ connectivity patterns, we may also have access to their shopping history on Amazon. In financial networks, in addition to the knowledge of trade relationships with other companies, there may be additional information available per company, e.g., the number of employees, category of products the company sales, or the annual profit. Gradraker can also incorporate this information by introducing additional learners based on the nodal feature vectors, and combine them as in (23).

V Performance analysis

To analyze the performance of the novel Gradraker algorithm, we assume that the following are satisfied.

(as1) For all sampled nodes $\{v_{t}\}_{t=1}^{T}$, the loss function {\cal L}({\mbox{\boldmath\theta}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t}) in (15) is convex w.r.t. $\theta$.

(as2) For $\bm{\theta}$ belonging to a bounded set $\Theta$ with \|{\mbox{\boldmath\theta}}\|\leq C_{\theta}, the loss is bounded; that is, {\cal L}({\mbox{\boldmath\theta}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t})\in[-1,1], and has bounded gradient, meaning, \|\nabla{\cal L}({\mbox{\boldmath\theta}}^{\top}{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{t}),y_{t})\|\leq L.

(as3) The kernels $\{\kappa_{p}\}_{p=1}^{P}$ are shift-invariant, standardized, and bounded, that is, $\kappa_{p}(\mathbf{a}_{n},\mathbf{a}_{n^{\prime}})\!\leq\!1,\,\forall\mathbf{a}_{n},\mathbf{a}_{n^{\prime}}$; and w.l.o.g. they also have bounded entries, meaning $\|\mathbf{a}_{n}\|\leq 1,\forall n$.

Convexity of the loss under (as1) is satisfied by the popular loss functions including the square loss and the logistic loss. As far as (as2), it ensures that the losses, and their gradients are bounded, meaning they are $L$-Lipschitz continuous. While boundedness of the losses commonly holds since \|{\mbox{\boldmath\theta}}\| is bounded, Lipschitz continuity is also not restrictive. Considering kernel-based regression as an example, the gradient is (\bm{\theta}^{\top}\mathbf{z}_{{\mathbf{V}}}(\mathbf{x}_{t})-y_{t})\mathbf{z}_{{\mathbf{V}}}(\mathbf{x}_{t})+\lambda{\mbox{\boldmath\theta}}. Since the loss is bounded, e.g., $\|\bm{\theta}^{\top}\mathbf{z}_{{\mathbf{V}}}(\mathbf{x}_{t})-y_{t}\|\leq 1$, and the RF vector in (12) can be bounded as $\|\mathbf{z}_{{\mathbf{V}}}(\mathbf{x}_{t})\|\leq 1$, the constant is $L:=1+\lambda C_{\theta}$ using the Cauchy-Schwartz inequality. Kernels satisfying the conditions in (as3) include Gaussian, Laplacian, and Cauchy [21]. In general, (as1)-(as3) are standard in online convex optimization (OCO) [30, 23], and in kernel-based learning [31, 21, 32].

In order to quantify the performance of Gradraker, we resort to the static regret metric, which quantifies the difference between the aggregate loss of an OCO algorithm, and that of the best fixed function approximant in hindsight, see also e.g., [30, 23]. Specifically, for a sequence $\{\hat{f}_{t}\}$ obtained by an online algorithm ${\cal A}$, its static regret is

[TABLE]

where $\hat{f}_{t}^{\rm RF}$ will henceforth be replaced by $\hat{f}_{t}$ for notational brevity; and, $f^{*}(\cdot)$ is defined as the batch solution

[TABLE]

where ${\cal F}_{p}:={\cal H}_{p}$, with ${\cal H}_{p}$ representing the RKHS induced by $\kappa_{p}$. We establish the regret of our Gradraker approach in the following lemma.

Lemma 1

*Under (as1), (as2), and with $\hat{f}_{p}^{*}$ defined as $\hat{f}_{p}^{*}(\cdot)\in\arg\min_{f\in\hat{\cal F}_{p}}\,\sum_{t=1}^{T}{\cal L}_{t}(f({\mathbf{a}}_{t}))$,

with \hat{\cal F}_{p}:=\{\hat{f}_{p}|\hat{f}_{p}({\mathbf{a}})={\mbox{\boldmath\theta}}^{\top}\mathbf{z}_{p}({\mathbf{a}}),\,\forall{\mbox{\boldmath\theta}}\in\mathbb{R}^{2D}\}, for any $p$, the sequences $\{\hat{f}_{p,t}\}$ and $\{\bar{w}_{p,t}\}$ generated by Gradraker satisfy the following bound*

[TABLE]

where {\mbox{\boldmath\theta}}_{p}^{*} is associated with the best RF function approximant \hat{f}_{p}^{*}({\mathbf{a}})=\left({\mbox{\boldmath\theta}}_{p}^{*}\right)^{\top}\mathbf{z}_{p}({\mathbf{a}}).

**Proof: ** See Appendix A

In addition to bounding the regret in the RF space, the next theorem compares the Gradraker loss relative to that of the best functional estimator in the original RKHS.

Theorem 2

Under (as1)-(as3), and with $f^{*}$ defined as in (V), for a fixed $\epsilon>0$, the following bound holds with probability at least 1-2^{8}\big{(}\frac{\sigma_{p}}{\epsilon}\big{)}^{2}\exp\big{(}\frac{-D\epsilon^{2}}{4N+8}\big{)}

[TABLE]

where $C$ is a constant, while $\sigma_{p}^{2}:=\mathbb{E}_{\pi_{\kappa_{p}}}[\|{\mathbf{v}}\|^{2}]$ is the second-order moment of the RF vector norm. Setting $\eta=\epsilon={\cal O}(1/\sqrt{T})$ in (2), the static regret in (24) leads to

[TABLE]

**Proof: ** See Appendix B

Observe that the probability of (2) to hold grows as $D$ increases, and one can always find a $D$ to ensure a positive probability for a given $\epsilon$. Theorem 2 establishes that with a proper choice of parameters, the Gradraker achieves sub-linear regret relative to the best static function approximant in (V), which means the novel Gradraker algorithm is capable of capturing the nonlinear relationship among nodal functions accurately, as long as enough nodes are sampled sequentially.

In addition, it is worth noting that Theorem 2 holds true regardless of the sampling order of the nodes $\{v_{1},\dots,v_{T}\}$. However, optimizing over the sampling pattern is possible, and constitutes one of our future research directions.

VI Numerical tests

In this section, Gradraker is tested on both synthetic and real datasets to corroborate its effectiveness. The tests will mainly focus on regression tasks for a fair comparison with existing alternatives.

VI-A Synthetic data test

Data generation. An Erdös-Rényi graph [33] with binary adjacency matrix ${\mathbf{A}}_{0}\in{\mathbf{R}}^{N\times N}$ was generated with probability of edge presence $\pi=0.2$, and its adjacency was symmetrized as ${{\mathbf{A}}}={\mathbf{A}}_{0}+{\mathbf{A}}_{0}^{\top}$. This symmetrization is not required by Gradraker, but it is necessary for alternative graph kernel based methods. A function over this graph was then generated with each entry of the coefficient vector {\mbox{\boldmath\alpha}}\in\mathbb{R}^{N} drawn uniformly from $[0.5,1]$, and each entry of the noise ${\mathbf{e}}$ drawn from $\mathcal{N}(0,0.01{\mathbf{I}})$. In each experiment, the sampling matrix $\Psi$ is randomly generated so that $M=0.05N$ of the nodes are randomly sampled, and the remaining $N-M$ nodes are treated as newly-joining nodes, whose function values and connectivity patterns are both unknown at the training phase, and whose nodal function values are estimated based on their connectivity with existing nodes in the network during the testing phase. All algorithms are carried out on the training set of $M$ nodes, and the obtained model is used to estimate the function value on the newly arriving nodes. The runtime for estimating the function value on the newly-joining nodes, as well as the generalization ${\rm NMSE}:=\frac{1}{|{\cal S}^{c}|}\|\hat{{\mathbf{x}}}_{{\cal S}^{c}}-{\mathbf{x}}_{{\cal S}^{c}}\|_{2}^{2}/\|{\mathbf{x}}_{{\cal S}^{c}}\|_{2}^{2}$ performance is evaluated, with ${\cal S}^{c}$ denoting the index set of new nodes. The Gradraker adopts a dictionary consisting of $2$ Gaussian kernels with parameters $\sigma^{2}=1,5$, using $D=10$ random features, and it is compared with: a) the $k$NN algorithm, with $k$ selected as the maximum number of neighbors a node has in a specific network, and with the combining weights set to $1/k$ in unweighted graphs, and $a_{il}/\sum_{j\in\mathcal{N}_{i}}a_{ij}$ for the $l$th neighbor in weighted graphs; b) the graph kernel (GK) based method using diffusion kernels with different bandwidths (named as GK-DF), or band-limited kernels with different bandwidths (GK-BL); and c) kernel based learning without RF approximation (KL) with a Gaussian kernel of $\sigma^{2}=5$. Results are averaged over $100$ independent runs. The regularization parameter for all algorithms is selected from the set $\mu=\{10^{-7},10^{-6},\dots,10^{0}\}$ via cross validation. Testing results. Figure 1 illustrates the performance in terms of the average runtime and NMSE versus the number of nodes (size) of the network. In this experiment, $\bar{{\mathbf{K}}}$ in (7) is generated from the normalized graph Laplacian ${\mathbf{L}}$, using the diffusion kernel $r(\lambda)=\exp(\sigma^{2}\lambda/2)$. A bandwidth of $\sigma^{2}=5$ was used to generate the data. It is observed that GK attains the best generalization accuracy when the ground-truth model is known, but its computational complexity grows rapidly with the network size. However, GK does not perform as well when a mismatched kernel is applied. The Gradraker method on the other hand, is very efficient, while at the same time it can provide reasonable estimates of the signal on the newly arriving nodes, even without knowledge about the kernels. The k-NN method is very efficient, but does not provide as reliable performance as the Gradraker.

Figure 2 depicts the performance of competitive algorithms. Matrix $\bar{{\mathbf{K}}}$ for data generation is formed based on (8) using the Gaussian kernel $\kappa({\mathbf{a}}_{i}-{\mathbf{a}}_{j})=\exp(\|{\mathbf{a}}_{i}-{\mathbf{a}}_{j}\|^{2}/\sigma^{2})$, with $\sigma^{2}=5$. In this case, KL exactly matches the true model, and hence it achieves the best performance. However, it is the most complex in terms of runtime. Meanwhile, GK-based methods suffer from model mismatch, and are also relatively more complex than Graderaker. The novel Gradraker is capable of estimating the nodal function on the newly joining nodes with high accuracy at very low computational complexity. Note that in real-world scenarios, accurate prior information about the underlying model is often unavailable, in which case Gradraker can be a more reliable and efficient choice.

VI-B Reconstruction of the temperature data

This subsection tests the performance of Gradraker on a real temperature dataset. The dataset comprises $24$ signals corresponding to the average temperature per month in the intervals $1961-1980$ and $1991-2010$ measured by $89$ stations in Switzerland [34]. The training set contains the first $12$ signals, corresponding to the interval $1961-1980$, while the test set contains the remaining $12$. Each station is represented by a node, and the graph was constructed using the algorithm in [35] based on the training signals. Given the test signal on a randomly chosen subset of $M$ vertices, the values at the remaining $N-M$ vertices are estimated as newly-coming nodes. The generalization NMSE over the $N-M$ nodes is averaged across the test signals.

Fig. 3 compares the performance of Gradraker with those of competing alternatives. Gradraker adopts a dictionary consisting of $3$ Gaussian kernels with parameters $\sigma^{2}=1,5,10$, using $D=100$ random features. It is clear from Fig. 3 that Gradraker outperforms GK in both generalization NMSE and runtime. On the other hand, even though KL achieves lower generalization NMSE, it incurs a much higher complexity.

VI-C Reconstruction of the Email-Eu-core data

The Eu-core network was generated using email data from a large European research institution [36], where each node represents a person, and an edge $(i,j)$ is present if person $i$ sent person $j$ at least one email. The e-mails only represent communication between institution members (the core), and the dataset does not contain incoming messages from or outgoing messages to the rest of the world. The dataset also contains “ground-truth” community memberships of the nodes. Each individual belongs to one of 42 departments at the research institute. During the experiment, the department labels are considered to be $y_{n}$ that are to be sampled and estimated. The graph consists of $N=1,005$ nodes, and $25,571$ edges. Gradraker adopts a dictionary consisting of $2$ Gaussian kernels with parameters $\sigma^{2}=1,10$, from which $D=10$ random features are generated. The test results were averaged over $100$ independent runs with randomly sampled nodes.

Fig. 4 compares the performance of Gradraker with those of alternative algorithms when different numbers of nodal labels are observed. It is clear that the RF-based approach outperforms the GK-based method in both reconstruction accuracy and runtime. While the batch KL method without RF approximation outperforms the RF method, it incurs considerably higher computational complexity.

VI-D Reconstruction of the Cora data

This subsection tests the Gradraker algorithm on the Cora citation dataset [5]. Gradraker adopts a dictionary consisting of $2$ Gaussian kernels with parameters $\sigma^{2}=1,10$, using $D=20$ random features. The results were averaged over $100$ independent runs. The Cora dataset consists of $2,708$ scientific publications classified into one of seven classes. The citation network consists of $5,429$ links. The network is constructed so that a link connects node $i$ to node $j$ if paper $i$ cites paper $j$, and the category id the paper belongs to is to be reconstructed. It can be observed again from Figure 5, that the Gradraker markedly outperforms the GK algorithms in terms of generalization NMSE, and is much more computationally efficient than all other algorithms except the kNN method, which however does not perform as well.

It can be readily observed from our numerical results over synthetic and real datasets, that the Gradraker provides reliable performance in terms of NMSE in all tests, while at the same time, it scales much better than all kernel based alternatives. This is because the alternative kernel-based algorithms require re-computing the kernel matrix whenever a new node joins the network. It is worth noting that all kernel-based alternatives require exact knowledge of the entire network topology, which is not necessary for GradRaker that only requires $\{{\mathbf{z}}_{{\mathbf{V}}}({\mathbf{a}}_{n})\}$. These tests corroborate the potential of GradRaker for application settings, where the graphs grow and nodes have privacy constraints.

VII Conclusions

The present paper deals with the problem of reconstructing signals over graphs, from samples over a subset of nodes. An online MKL based algorithm is developed, which is capable of estimating and updating the nodal functions even when samples are collected sequentially. The novel online scheme is highly scalable and can estimate the unknown signals on newly joining nodes. Unlike many existing approaches, it only relies on encrypted nodal connectivity information, which is appealing for networks where nodes have strict privacy constraints.

This work opens up a number of interesting directions for future research, including: a) exploring distributed implementations that are well motivated in large-scale networks; b) graph-adaptive learning when multiple sets of features are available; and c) developing adaptive sampling strategies for Gradraker.

Appendix A Proof of Lemma 1

To prove Lemma 1, we introduce two intermediate lemmata.

Lemma 3

Under (as1), (as2), and $\hat{f}_{p}^{*}$ as in (V) with {\cal F}_{p}:=\{\hat{f}_{p}|\hat{f}_{p}({\mathbf{a}})={\mbox{\boldmath\theta}}^{\top}\mathbf{z}_{p}({\mathbf{a}}),\,\forall{\mbox{\boldmath\theta}}\in\mathbb{R}^{2D}\}, let $\{\hat{f}_{p,t}({\mathbf{a}}_{t})\}$ denote the sequence of estimates generated by Gradraker with a pre-selected kernel $\kappa_{p}$. Then the following bound holds true w.p.1

[TABLE]

where $\eta$ is the learning rate, $L$ is the Lipschitz constant in (as2), and {\mbox{\boldmath\theta}}_{p}^{*} is the corresponding parameter (or weight) vector supporting the best estimator \hat{f}_{p}^{*}({\mathbf{a}})=({\mbox{\boldmath\theta}}_{p}^{*})^{\top}\mathbf{z}_{p}({\mathbf{a}}).

**Proof: ** The proof is similar to the regret analysis of online gradient descent, see e.g., [22].

In addition, we will bound the difference between the loss of the solution obtained from Algorithm 2 and the loss of the best single kernel-based online learning algorithm. Specifically, the following lemma holds.

Lemma 4

Under (as1) and (as2), with $\{\hat{f}_{p,t}\}$ generated from Gradraker, it holds that

[TABLE]

where $\eta$ is the learning rate in (22), and $P$ is the number of kernels in the dictionary.

**Proof: ** Letting $W_{t}:=\sum_{p=1}^{P}w_{p,t}$, the weight recursion in (22) implies that

[TABLE]

where the last inequality holds because $\exp(-\eta x)\leq 1-\eta x+\eta^{2}x^{2}$, for $|\eta|\leq 1$. Furthermore, substituting $\bar{w}_{p,t}:=w_{p,t}/\sum_{p=1}^{P}w_{p,t}=w_{p,t}/W_{t}$ into (A) leads to

[TABLE]

Since $1+x\leq e^{x},\,\forall x$, it follows that

[TABLE]

Telescoping (A) from $t=1$ to $T$ yields

[TABLE]

On the other hand, for any $p$, it holds that

[TABLE]

Combining (A) with (35), we arrive at

[TABLE]

Taking the logarithm on both sides of (A), and recalling that $w_{p,1}=1/P$, we obtain

[TABLE]

Re-organizing the terms leads to

[TABLE]

and the proof is complete, since ${\cal L}_{t}\left(\hat{f}_{p,t}({\mathbf{a}}_{t})\right)^{2}\leq 1$ and $\sum_{p=1}^{P}\bar{w}_{p,t}=1$.

Since ${\cal L}_{t}(\cdot)$ is convex under (as1), Jensen’s inequality implies

[TABLE]

Combining (39) with Lemma 4, one arrives readily at

[TABLE]

where (a) follows due to Lemma 3 and because {\mbox{\boldmath\theta}}_{p}^{*} is the optimal solution for any given $\kappa_{p}$. This proves Lemma 1.

Appendix B Proof of Theorem 2

To bound the performance relative to the best estimator $f^{*}({\mathbf{a}}_{t})$ in the RKHS, the key step is to bound the approximation error. For a given shift-invariant $\kappa_{p}$, the maximum point-wise error of the RF kernel approximant is bounded with probability at least 1-2^{8}\big{(}\frac{\sigma_{p}}{\epsilon}\big{)}^{2}\exp\big{(}\frac{-D\epsilon^{2}}{4N+8}\big{)}, by [21]

[TABLE]

where $\epsilon>0$ is a given constant, $D$ the number of features, while $M$ is the number of nodes already in the network, and $\sigma_{p}^{2}:=\mathbb{E}_{p}[\|{\mathbf{v}}\|^{2}]$ is the second-order moment of the RF vector norm induced by $\kappa_{p}$. Henceforth, for the optimal function estimator (V) in ${\cal H}_{p}$ denoted by $f_{p}^{*}({\mathbf{a}}):=\sum_{t=1}^{T}\alpha_{p,t}^{*}\kappa_{p}({\mathbf{a}},{\mathbf{a}}_{t})$, and its RF-based approximant $\check{f}^{*}_{p}:=\sum_{t=1}^{T}\alpha_{p,t}^{*}{\mathbf{z}}_{p}^{\top}({\mathbf{a}}){\mathbf{z}}_{p}({\mathbf{a}}_{t})\in{\cal F}_{p}$, we have

[TABLE]

where (a) is due to the triangle inequality; (b) uses the Lipschitz continuity of the loss, and (c) is due to the Cauchy-Schwarz inequality. Combining with (41), yields

[TABLE]

where we used that $C:=\max_{p}\sum_{t=1}^{T}|\alpha_{p,t}^{*}|$. Under the kernel bounds in (as3), the uniform convergence in (41) implies that $\sup_{{\mathbf{a}}_{t},{\mathbf{a}}_{t^{\prime}}\in{\cal X}}{\mathbf{z}}_{p}^{\top}({\mathbf{a}}_{t}){\mathbf{z}}_{p}({\mathbf{a}}_{t^{\prime}})\leq 1+\epsilon$, w.h.p., which leads to

[TABLE]

where for the last inequality we used the definition of $C$.

Lemma 1 together with (43) and (B) lead to the regret of the proposed Gradraker algorithm relative to the best static function in ${\cal H}_{p}$, that is given by

[TABLE]

which completes the proof of Theorem 2.