1

Decentralized Online Learning with Kernels
Alec Koppel§ , Santiago Paternain? , Cédric Richard† and Alejandro Ribeiro?

arXiv:1710.04062v1 [math.OC] 11 Oct 2017

Abstract
We consider multi-agent stochastic optimization problems over reproducing kernel Hilbert spaces (RKHS).
In this setting, a network of interconnected agents aims to learn decision functions, i.e., nonlinear statistical
models, that are optimal in terms of a global convex functional that aggregates data across the network, with
only access to locally and sequentially observed samples. We propose solving this problem by allowing each agent
to learn a local regression function while enforcing consensus constraints. We use a penalized variant of functional
stochastic gradient descent operating simultaneously with low-dimensional subspace projections. These subspaces
are constructed greedily by applying orthogonal matching pursuit to the sequence of kernel dictionaries and weights.
By tuning the projection-induced bias, we propose an algorithm that allows for each individual agent to learn, based
upon its locally observed data stream and message passing with its neighbors only, a regression function that is
close to the globally optimal regression function. That is, we establish that with constant step-size selections agents’
functions converge to a neighborhood of the globally optimal one while satisfying the consensus constraints as
the penalty parameter is increased. Moreover, the complexity of the learned regression functions is guaranteed to
remain finite. On both multi-class kernel logistic regression and multi-class kernel support vector classification with
data generated from class-dependent Gaussian mixture models, we observe stable function estimation and state of
the art performance for distributed online multi-class classification. Experiments on the Brodatz textures further
substantiate the empirical validity of this approach.

I. I NTRODUCTION
We consider decentralized online optimization problems: a network G = (V, E) of agents aims to
minimize a global objective that is a sum of local convex objectives available only to each node.
The problem is online and distributed because data samples upon which the local objectives depend
are sequentially and locally observed by each agent. In this setting, agents aim to make inferences
as well as one which has access to all data at a centralized location in advance. Instead of assuming
agents seek a common parameter vector w ∈ Rp , we focus on the case where agents seek to learn a
common decision function f (x) that belong to a reproducing kernel Hilbert space (RKHS). Such functions
represent, e.g., nonlinear statistical models [2] or trajectories in a continuous space [3]. Learning in multiagent settings arises predominately in two technological settings: industrial-scale machine learning, where
optimizing statistical model parameters is decentralized across a parallel processing architecture to attain
computational speedup; and networked intelligent systems such as sensor networks [4], multi-robot teams
[5], [6], and Internet of Things [7], [8]. In the later setting, decentralized processing justified as opposed
to using a fusion center when the communication cost of centralization exceeds the cost of distributed
information protocols. This is true of multi-agent systems with streaming data considered here.
Efforts to develop optimization tools for multi-agent online learning have thus far been restricted to
the case where each agent learns a linear statistical model [9] or a task-driven dictionary [10] that is
as good as one with data aggregated across the network. However, these efforts exclude the state of the
art tools for statistical learning based on nonlinear interpolators: namely, kernel methods [11], [12] and
neural networks [13], [14]. We note that instabilities associated with non-convexity which are only a
minor issue in centralized settings [15] become both theoretically and empirically difficult to overcome
This work in this paper is supported by NSF CCF-1017454, NSF CCF-0952867, ONR N00014-12-1-0997, ARL MAST CTA, and ASEE
SMART.Part of the results in this paper appeared in [1].
§
Computational and Information Sciences Directorate, U.S. Army Research Laboratory, Adelphi, MD, 20783. Email:
alec.e.koppel.civ@mail.mil
?
Department of ESE, University of Pennsylvania, 200 South 33rd Street, Philadelphia, PA 19104. Email: {spater, aribeiro}@seas.upenn.edu
†
Laboratory Lagrange - UMR CNRS 7293, Observatory of the French Riviera University of Nice Sophia-Antipolis, Nice, France, 06108

2

in settings with consensus constraints [10], and therefore efforts to extend neural network learning to
multi-agent online learning likely suffer the same drawbacks.1 Therefore, we focus on extending kernel
methods to decentralized online settings, motivated both by its advantageous empirical performance, as
well as the theoretical and practical benefits of the fact that the optimization problem defined by their
training is convex. This stochastic convex problem, however, is defined over an infinite dimensional space,
and therefore it is not enough to solve the optimization problem, but one must also solve it in an optimally
sparse way. Doing so in multi-agent settings is the goal of this work.
To contextualize our solution methodology, consider centralized vector-valued stochastic convex programming, which has classically been solved with stochastic gradient descent (SGD) [16]. SGD involves
descending along the negative of the stochastic gradient rather than the true gradient to avoid the fact that
computing the gradient of the average objective has complexity comparable to the training sample size,
which could be infinite. In contrast, the setting considered in this work is a stochastic program defined
over a function space, which is in general an intractable variational inference problem. However, when
the function space is a RKHS [17], the Representer Theorem allows us to transform a search over an
infinite space into one over a set of weights and data samples [18]. Unfortunately, the feasible set of the
resulting problem has complexity comparable to the sample size N , and thus is intractable for N → ∞
[19]. Compounding this problem is that the storage required to construct the functional generalization of
SGD is comparable to the iteration index of the algorithm, which is untenable for online settings.
Efforts to mitigate the complexity of the function representation (“the curse of kernelization”) have
been previously developed. These combine functional extensions of stochastic gradient method with
compressions of the function parameterization independently of the optimization problem to which they are
applied [20]–[24] or approximate the kernel during training [25]–[29], and at best converge on average.
In contrast, a method was recently proposed that combines greedily constructed [30] sparse subspace
projections with functional stochastic gradient method and guarantees exact convergence to the minimizer
of the average risk functional. This technique, called parsimonious online learning with kernels (POLK),
tailors the parameterization compression to preserve the descent properties of the underlying RKHS-valued
stochastic process [31], and inspires the approach considered here.
In this work, we extend the ideas in [31] to multi-agent settings. Multiple tools from distributed
optimization may be used to do so; however, we note that the Representer Theorem [18] has not been
established for general stochastic saddle point problems in RKHSs. Therefore, we adopt an approximate
primal-only approach based on penalty methods [32], [33], which in decentralized optimization is known
as distributed gradient descent (DGD). Using functional stochastic extensions of DGD, together with the
greedy Hilbert subspace projections designed in POLK, we develop a method such that each agent, through
its local data stream and message passing with only its neighbors, learns a memory-efficient approximation
to the globally optimal regression function with probability 1. Such global stability guarantees are in
contrast to specialized results for multi-agent kernel learning [34], [35] and alternative distributed online
nonlinear function estimation methods such as dictionary learning [10], [15], [36] or neural networks [14],
which suffer from instability due to the non-convexity of the optimization problem their training defines.
The result of the paper is organized as follows. In Section II we clarify the problem setting of stochastic
programming in RKHSs in the centralized and decentralized case. In Section III, we propose a new
penalty functional that permits deriving a decentralized online method for kernel regression without any
complexity bottleneck by making use of functional stochastic gradient method (Section III-A) combined
with greedy subspace projections (Section III-B). In Section IV we present our main theoretical results,
which establishes that the function sequence of each agent generated by the proposed technique converges
to a neighborhood of the globally optimal function with probability 1. In Section V, we present numerical
examples of decentralized online multi-class kernel logistic regression and kernel support vector machines
with data generated from Gaussian mixtures, and observe a state of the art trade-off between Lyapunov
1
In general, globally convergent decentralized online training of neural networks is an open problem, whose solution requires fundamentally
new approaches to stochastic global optimization.

3

stability and statistical accuracy. We then apply the resulting method to the benchmark Brodatz texture
dataset [37] and observe state of the art decentralized online multi-class classification performance.
II. P ROBLEM F ORMULATION
A. Decentralized Functional Stochastic Programming
Consider the problem of expected risk minimization, where the goal is to learn a regressor that minimizes
a loss function quantifying the merit of a statistical model averaged over a data set. We focus on the case
when the number of training examples N is very large or infinite. In this work, input-output examples,
(xn , yn ), are i.i.d. realizations drawn from a stationary joint distribution over the random pair (x, y) ∈
X × Y, where X ⊂ Rp and Y ⊂ R. Here, we consider finding regressors that are not vector valued
parameters, but rather functions f˜ ∈ H in a hypothesized function class H, which allows for learning
nonlinear statistical models rather than generalized linear models that rarely achieve satisfactory statistical
error rates in practice [12], [38]. The merit of the function f˜ is evaluated by the convex loss function
` : H×X ×Y → R that quantifies the merit of the estimator f˜(x̃) evaluated at feature vector x̃. This loss is
averaged over all possible training examples to define the statistical loss L̃(f˜) := Ex,y [`(f˜(x), y)], which
we combine with a Tikhonov regularizer to construct the regularized loss R̃(f˜) := argminf˜∈H L̃(f˜) +
(λ/2)kf˜k2H [39], [40]. We then define the optimal function as
h
i λ
(1)
f˜∗ = argmin R̃(f˜) := argmin Ex̃,ỹ `(f˜ x̃), ỹ + kf˜k2H
2
f˜∈H
f˜∈H
In this work, we focus on extensions of the formulation in (1) to the case where data is scattered across
an interconnected network that represents, for instance, robotic teams [10], communication systems [41],
or sensor networks [4]. To do so, we define a symmetric, connected, and directed network G = (V, E)
with |V| = V nodes and |E| = E edges and denote as ni := {j : (i, j) ∈ E} the neighborhood of agent
i. For simplicity we assume that the number of edges E is even. Each agent i ∈ V observes a local data
sequence as realizations (xi,n , yi,n ) from random pair (xi , yi ) ∈ X × Y and seeks to learn a common
globally optimal regression function f . This setting may be mathematically captured by associating to
each node i a convex loss functional `i : H × X × Y → R that quantifies the merit of the estimator fi (xi )
evaluated at feature vector xi , and defining the goal for each node as the minimization of the common
global loss

h
X
i λ
∗
2
f = argmin
Exi ,yi `i (f xi ), yi + kf kH
(2)
2
f ∈H
i∈V
Observe that this global loss is a network-wide average (scaled by V ) of all local losses, and therefore the
minimizers of (1) and (2) coincide when (xi , yi ) have a common joint distribution for each i. However,
in multi-agent optimization, this is not generally the case, thus when selecting a regression function f
with only local data, different agents will learn a different decision function fi∗ that it is not optimal
as compared to one selected in a centralized manner, i.e., with the data gathered by all agents. To
overcome this limitation we allow message passing between agents and we impose a consensus constraint
on the regression function among neighbors fi = fj , (i, j) ∈ E. Thus we consider the nonparametric
decentralized stochastic program:

h
X
i λ
∗
2
f = argmin
Exi ,yi `i (fi x), yi + kfi kH
2
{fi }⊂H
i∈V
such that

fi = fj , (i, j) ∈ E

(3)

For further define the product Hilbert space HV of functions aggregated over the network whose elements
are stacked functions f (·) = [f1 (·); · · · ; fV (·)] that yield vectors of length V when evaluated at local
random vectors f (x) = [f1 (x1 ); · · · ; fV (xV )] ∈ RV . Moreover, define the stacked random vectors x =

4

[x1 ; · · · ; xV ] ∈ X V ⊂ RV p and y = [y1 ; · · · yV ] ∈ RV that represents V labels or physical measurements,
for instance.
The goal of this paper is to develop an algorithm to solve (3) in distributed online settings where nodes
do not know the distribution of the random pair (xi , yi ) but observe local independent training examples
(xi,n , yi,n ) sequentially.
B. Function Estimation in Reproducing Kernel Hilbert Spaces
The optimization problem in (1), and hence (3), is intractable in general, since it defines a variational
inference problem integrated over the unknown joint distribution P(x, y). However, when H is equipped
with a reproducing kernel κ : X × X → R (see [12], [42]), a function estimation problem of the form
(1) may be reduced to a parametric form via the Representer Theorem [19], [43]. Thus, we restrict the
Hilbert space in Section II-A to be one equipped with a kernel κ that satisfies for all functions f˜ : X → R
in H:
(i) hf˜, κ(xi , ·))iH = f˜(xi ), (ii) H = span{κ(xi , ·)}
(4)
for all xi ∈ X . Here h·, ·iH denotes the Hilbert inner product for H. Further assume that the kernel
is positive semidefinite, i.e. κ(xi , x0i ) ≥ 0 for all xi , x0i ∈ X . Function spaces of this type are called
reproducing kernel Hilbert spaces (RKHS).
In (4), property (i) is the reproducing property (via Riesz Representation Theorem [43]). Replacing f˜
by κ(x0i , ·) in (4) (i) yields hκ(x0i , ·), κ(xi , ·)iH = κ(xi , x0i ) which is the origin of the term “reproducing
kernel.” This property induces a nonlinear transformation of the input space X : denote by φ(·) a nonlinear
map of the feature space that assigns to each xi the kernel function κ(·, xi ). The reproducing property
yields that the inner product of the image of distinct feature vectors xi and x0i under the map φ requires
only kernel evaluations: hφ(xi ), φ(x0i )iH = κ(xi , x0i ) (the ’kernel trick’).
Moreover, property (4) (ii) states that functions f˜ ∈ H may be written as a linear combination of kernel
evaluations. For kernelized and regularized empirical risk minimization (ERM), the Representer Theorem
[17], [18] establishes that the optimal f˜ in hypothesized function class H admit an expansion in terms of
kernel evaluations only over training examples
N
X
˜
f (xi ) =
wi,n κ(xi,n , xi ) ,
(5)
n=1
T

N

where wi = [wi,1 , · · · , wi,N ] ∈ R denotes a set of weights. The upper index N in (5) is referred
to as the model order, and for ERM the model order and training sample size are
d equal. Common
choices κ include the polynomial and radial basis kernels, i.e., κ(xi , x0i ) = xTi x0i + b and κ(xi , x0i ) =
exp{−kxi − x0i k22 /2d2 }, respectively, where xi , x0i ∈ X .
Suppose, for the moment, that we have access to N i.i.d. realizations of the random pairs (xi , yi ) for
each agent i such that the expectation in (3) is computable, and we further ignore the consensus constraint.
Then the objective in (3) becomes:
N
1 XX
λ
`(fi (xi,n ), yi,n ) + kfi k2H
f = argmin
2
f ∈HV N n=1 i∈V
∗

(6)

Then, by substituting the Representer Theorem [cf. (5)] into (3), we obtain that optimizing in HV reduces
to optimizing over the set of N V weights:
N
1XX
λ
f ∗=argmin
`i (wiTκXi(xi,n),yi,n)+ wiT KXi ,Xiwi ,
2
{wi}∈RN Nn=1 i∈V

(7)

where we have defined the Gram (or kernel) matrix KXi ,Xi ∈ RN ×N , with entries given by the kernel
evaluations between xi,m and xi,n as [KXi ,Xi ]m,n = κ(xi,m , xi,n ). We further define the vector of kernel

5

evaluations κXi (·) = [κ(xi,1 , ·) . . . κ(xi,N , ·)]T , which are related to the kernel matrix as KXi ,Xi =
[κXi (xi,1 ) . . . κXi (xi,N )]. The dictionary of training points associated with the kernel matrix is defined as
Xi = [xi,1 , . . . , xi,N ].
By exploiting the Representer Theorem, we transform a nonparametric infinite dimensional optimization problem in HV (6) into a finite N V -dimensional parametric problem (7). Thus, for empirical risk
minimization, the RKHS provides a principled framework to solve nonparametric regression problems as
a search over RV N for an optimal set of coefficients.
However, to solve problems of the form (6) when training examples (xi,n , yi,n ) become sequentially
available or their total number N is not finite, the objective in (6) becomes an expectation over random
pairs (xi , yi ) as [11]
X
X
f ∗ = argmin
Exi ,yi [`i (
wi,n κ(xi,n , xi ), yi )]
wi ∈RI ,{xi,n }n∈I i∈V

n∈I

λ X
wi,n wi,m κ(xi,m , xi,n )k2H ,
+ k
2 n,m∈I

(8)

where we substitute the Representer Theorem generalized to the infinite sample-size case established in
[19] into the objective (3) with I as some countably infinite indexing set. That is, as the data sample size
N → ∞, the representation of fi becomes infinite as well. Thus, our goal is to solve (8) in an approximate
manner such that each fi admits a finite representation near fi∗ , while satisfying the consensus constraints
fi = fj for (i, j) ∈ E (which were omitted for the sake of discussion between (6) - (8)).
III. A LGORITHM D EVELOPMENT
We turn to developing an online iterative and decentralized solution to solving (3) when the functions
{fi }i∈V are elements of a RKHS, as detailed in Section II-B. To exploit the properties of this function
space, we require the applicability of the Representer Theorem [cf. (5)], but this result holds for any
regularized minimization problem with a convex functional. Thus, we may address the consensus constraint
fi = fj , (i, j) ∈ E in (3) by enforcing approximate consensus on estimates fi (xi ) = fj (xi ) in expectation.
This specification may be met by introducing the penalty functional
h
X
i λ
ψc (f ) =
Exi ,yi `i (fi xi ), yi + kfi k2H
2
i∈V


cX
2
(9)
+
Ex [fi (xi )−fj (xi )]
2j∈n i
i

The reasoning for the definition (9) rather than one that directly addresses the consensus constraint
deterministically is given in Remark 1, motivated by following the algorithm derivation. For future
reference, we also define the local penalty as
h
i λ
ψi,c (fi ) = Exi ,yi `i (fi xi ), yi + kfi k2H
2

cX
+
Exi [fi (xi )−fj (xi )]2
(10)
2j∈n
i
P
and we observe from (9) - (10) that ψc (f ) = i ψi,c (fi ). Further define fc∗ = argminf ∈HV ψc (f ). We
note that in the vector-valued decision variable case, other techniques to address the constraint in (3) are
possible such as primal-dual methods [9] or dual methods [44], but the Representer Theorem has not been
established for RKHS-valued stochastic saddle point problems. It is an open question whether expressions
of the form (5) apply to problems with general functional constraints, but this matter is beyond the scope
of this work. Therefore, these other approaches which make use of Lagrange duality do not readily extend
to the nonparametric setting considered here.

6

Algorithm 1 Greedy Projected Penalty Method
Require: {xt , yt , ηt , t }t=0,1,2,...
initialize fi,0 (·) = 0, Di,0 = [], w0 = [], i.e. initial dictionary, coefficients are empty for each i ∈ V
for t = 0, 1, 2, . . . do
loop in parallel for agent i ∈ V
Observe local training example realization (xi,t , yi,t )
Send obs. xi,t to nodes j ∈ ni , receive scalar fj,t (xi,t )
Receive obs. xj,t from nodes j ∈ ni , send fi,t (xj,t )
Compute unconstrained stochastic grad. step [cf. (22)]
f˜i,t+1 (·) = (1 − ηt λ)fi,t − ηt ∇fi ψ̂i,c (fi (xi,t ), yi,t ) .
Update params: D̃i,t+1 = [Di,t , xi,t ], w̃i,t+1 [cf. (23)]
Greedily compress function using matching pursuit
(fi,t+1, Di,t+1,wi,t+1) = KOMP(f˜i,t+1 ,D̃i,t+1 ,w̃i,t+1 ,t)
end loop
end for

A. Functional Stochastic Gradient Method
Given that the data distribution P(x, y) is unknown, minimizing ψc (f ) directly via variational inference is not possible. Rather than postulate a specific distribution for (x, y), we only assume access to
sequentially available (streaming) independent and identically distributed samples (xt , yt ) from their joint
density. Then, we may wield tools from stochastic approximation to minimize (9), which in turn yields
a solution to (3). Begin by defining, ψ̂c (f (xt ), yt ), the stochastic approximation of the penalty function
ψc (f ), evaluated at a realization (xt , yt ) of the stacked random pair (x, y):
X
 λ
`i (fi xi,t ), yi,t + kfi k2H
ψ̂c (f (xt ), yt ) =
2
i∈V

X
c
+
(fi (xi,t )−fj (xi,t ))2
(11)
2 j∈n
j

and the local instantaneous penalty function ψ̂i,c (fi (xi,t ), yi,t ) similarly. To compute the functional stochastic gradient of ψc (f ) evaluated at a sample point (xt , yt ), we first address the local loss `i (fi xi,t ), yi,t )
in (11) as [22], [31]:
∇fi `i (fi (xi,t ), yi,t )(·) =

∂`i (fi (xi,t ), yi,t ) ∂fi (xi,t )
(·)
∂fi (xi,t )
∂fi

(12)

where we have applied the chain rule. Now, define the short-hand notation
`0i (fi (xi,t ), yi,t ) := ∂`i (fi (xi,t ), yi,t )/∂fi (xi,t )
for the derivative of `i (f (xi,t ), yi,t ) with respect to its first scalar argument fi (xi,t ) evaluated at xi,t . To
evaluate the second term on the right-hand side of (12), differentiate both sides of the expression defining
the reproducing property of the kernel [cf. (4)(i)] with respect to fi to obtain
∂fi (xi,t )
∂hfi , κ(xi,t , ·)iH
=
= κ(xi,t , ·)
∂fi
∂fi

(13)

7

Then, given (12) - (13), we may compute the overall gradient of the instantaneous penalty function
ψ̂c (f (xt ), yt ) in (11) as
h
∇f ψ̂c (f (xt ), yt ) = vec `0i (fi (xi,t ), yi,t )κ(xi,t , ·)+ λfi
(14)
i
X
+c (fi (xi,t )−fj (xi,t ))κ(xi,t , ·)
j∈ni

where on the right-hand side of (14), we have defined the vector stacking notation vec[·] to denote the
stacking of V component-wise functional gradients, each associated with function fi , i ∈ V, and used the
fact that the variation of the instantaneous approximate of the cross-node term, [fi (xi )−fj (xi )]2 , by the
same reasoning as (12) - (13), is 2[fi (xi,t )−fj (xi,t )]κ(xi,t , ·). With this computation in hand, we present
the stochastic gradient method for the λ-regularized multi-agent expected risk minimization problem in
(3) as
h
ft+1 = (1 − ηt λ)ft − ηt vec `0i (fi,t (xi,t ), yi,t )κ(xi,t , ·)
i
X
+c (fi,t (xi,t )−fj,t (xi,t ))κ(xi,t , ·) ,
(15)
j∈ni

where ηt > 0 is an algorithm step-size either chosen as diminishing with O(1/t) or a small constant –
see Section IV. We may glean from (15) that the update for the network-wide P
function ft decouples into
ones for each agent i ∈ V, using the node-separability of the penalty ψc (f ) = i ψi,c (fi ), i.e.,
h
fi,t+1 = (1 − ηt λ)fi,t − ηt `0i (fi,t (xi,t ), yi,t )κ(xi,t , ·)
i
X
+c (fi,t (xi,t )−fj,t (xi,t ))κ(xi,t , ·) .
(16)
j∈ni

We further require that, given λ > 0, the step-size satisfies ηt < 1/λ and the global sequence is initialized
as f0 = 0 ∈ HV . With this initialization, the Representer Theorem (5) implies that, at time t, the function
fi,t admits an expansion in terms of feature vectors xi,t observed thus far as
fi,t (x) =

t−1
X

T
wi,n κ(xi,n , x) = wi,t
κXi,t (x) .

(17)

n=1

On the right-hand side of (17) we have introduced the notation Xi,t = [xi,1 , . . . , xi,t−1 ] ∈ Rp×(t−1) ,
κXi,t (·) = [κ(xi,1 , ·), . . . , κ(xi,t−1 , ·)]T , and wi,t = [wi,1 , . . . wi,t−1 ] ∈ Rt−1 . Moreover, observe that the
kernel expansion in (17), taken together with the functional update (15), yields the fact that performing the
stochastic gradient method in HV amounts to the following V parallel parametric updates on the kernel
dictionaries Xi and coefficients wi :
Xi,t+1 = [Xi,t , xi,t ] ,
(
(1 
− ηt λ)[wi,t ]u for 0 ≤ u ≤ t − 1

P
[wi,t+1 ]u=
0
−ηt `i (fi,t (xi,t),yi,t)+c j∈ni(fi,t (xi,t)−fj,t (xi,t)) ,

(18)

where the second case on the last line of (18) is for u = t. Observe that this update causes Xi,t+1 to have
one more column than Xi,t . We define the model order as number of data points Mi,t in the dictionary
of agent i at time t (the number of columns of Xt ). FSGD is such that Mi,t = t − 1, and hence grows
unbounded with iteration index t. Next we address this intractable memory growth such that we may
execute stochastic descent through low-dimensional projections of the stochastic gradient, inspired by
[31]. First, we clarify the motivation for the choice of the penalty function (9).

8

Remark 1 In principle, it is possible to address the RKHS-valued consensus constraint in (3) directly,
through primal-only stochastic methods, by introducing the penalty function

h
X
i λ
cX
kfi −fj k2H
(19)
ψ̃c (f ) =
Exi ,yi `i (fi xi ), yi + kfi k2H +
2
2
j∈n
i∈V
i

Observe, however, that FSGD applied to (19), using comparable reasoning to that which leads to (16)
from (9), yields
h
fi,t+1 = (1 − ηt λ)fi,t − ηt ∇fi `0i (fi,t (xi,t ), yi,t )κ(xi,t , ·)
i
X
+c (fi,t −fj,t ) .
(20)
j∈ni

Unfortunately, we cannot inductively define a parametric representation of (20) for node i in terms of its
own kernel dictionaries and weights independently of the entire function associated to node j, since the
last term in (20) lives directly in the Hilbert space. Thus, to implement (20) each agent would need to store
the entire kernel dictionary and weights of all its neighbors at each step, which is impractically costly. The
use of (9) rather than (19) is further justified that under a hypothesis regarding the mean transformation
of the local data spaces, Exi [κ(xi , ·)], consensus with respect to the Hilbert norm, in addition to the mean
square sense, is achieved when the penalty coefficient is c → ∞ (see Section IV for details).
B. Sparse Subspace Projections
To mitigate the complexity growth noted in Section III-A, we approximate the function sequence
(15) by one that is orthogonally projected onto subspaces HD ⊆ H that consist only of functions that
can
some dictionary D = [d1 , . . . , dM ] ∈ Rp×M , i.e., HD = {f : f (·) =
PM be represented using
T
M
n=1 wn κ(dn , ·) = w κD (·)} = span{κ(dn , ·)}n=1 , and {dn } ⊂ {xu }u≤t . For convenience we define
[κD (·) = κ(d1 , ·) . . . κ(dM , ·)], and KD,D as the resulting kernel matrix from this dictionary. We enforce
function parsimony by selecting dictionaries Di with Mi,t << O(t) for each i [31].
To be specific, we propose replacing the local update (16) in which the dictionary grows at each iteration
Mt+1
by its projection onto subspace HDi,t+1 = span{κ(di,n , ·)}n=1
as
2

fi,t+1 = argmin f − fi,t −ηt ∇fi ψ̂i,c (fi (xi,t ), yi,t )
H

f ∈HDi,t+1

h

:= PHDi,t+1 (1 − ηt λ)fi,t − ηt ∇fi `i (fi,t (xi,t ), yi,t )
i
X
+c (fi,t (xi,t )−fj,t (xi,t ))κ(xi,t , ·) .

(21)

j∈ni

where we define the projection operator P onto subspace HDi,t+1 ⊂ H by the update (21).
Coefficient update The update (21), for a fixed dictionary Di,t+1 ∈ Rp×Mt+1 , yields one in the coefficient
space only. This fact may be observed by defining the un-projected stochastic gradient step starting at
function fi,t parameterized by dictionary Di,t and coefficients wi,t :
f˜i,t+1 = fi,t − ηt ∇fi ψ̂i,c (fi (xi,t ), yi,t ) .

(22)

This update may be represented using dictionary and weights
D̃i,t+1 = [Di,t , xi,t ] ,
(
(1 
− ηt λ)[wi,t ]u for 0 ≤ u ≤ t − 1

P
[w̃i,t+1 ]u=
0
−ηt `i (fi,t (xi,t),yi,t)+c j∈ni(fi,t (xi,t)−fj,t (xi,t)) ,

(23)

9

where the last coefficient is for u = t. Note that D̃i,t+1 has M̃ = Mi,t +1 columns, which is also the length
of w̃i,t+1 . For a fixed Di,t+1 , the stochastic projection (21) is a least-squares update on the coefficient
vector: the Representer Theorem allows us to rewrite (21) in terms of kernel expansions as in Section 3.2
of [31], which yields
wi,t+1 = K−1
(24)
Di,t+1 Di,t+1 KDi,t+1 D̃i,t+1 w̃i,t+1 ,
where we define the cross-kernel matrix KDi,t+1 ,D̃i,t+1 whose (n, m)th entry is given by κ(di,n , d̃i,m ). The
other kernel matrices KD̃i,t+1 ,D̃i,t+1 and KDi,t+1 ,Di,t+1 are defined similarly. Observe that Mi,t+1 is the
number of columns in Di,t+1 , while M̃i = Mi,t + 1 is the number of columns in D̃t+1 [cf. (23)]. Given
that the local projections of f˜i,t+1 onto stochastic subspaces HDi,t+1 , for a fixed node-specific dictionaries
Di,t+1 , is a least-squares problem, we now detail the kernel dictionary Di,t+1 selection from past data
{xi,u , yi,u }u≤t .
Dictionary Update The selection procedure for the kernel dictionary Di,t+1 is based upon greedy compression [45]: function f˜i,t+1 defined by the stochastic gradient method without projection is parameterized
by dictionary D̃i,t+1 [cf. (23)] of model order M̃i = Mi,t + 1. We form Di,t+1 by selecting a subset of
Mi,t+1 columns from D̃i,t+1 that best approximate f˜i,t+1 in terms of Hilbert norm error, which may be
done by executing kernel orthogonal matching pursuit (KOMP) [30], [46] with error tolerance t to find
a kernel dictionary matrix Di,t+1 based on the one which adds the latest sample point D̃i,t+1 . This choice
is due to the fact that we can tune its stopping criterion to guarantee stochastic descent, and guarantee
the model order of the learned function remains finite – see Section IV for details.
We now describe the variant of KOMP we propose using, called Destructive KOMP with Pre-Fitting
(see [46], Section 2.3). Begin with an input a candidate function f˜ of model order M̃ parameterized by
kernel dictionary D̃ ∈ Rp×M̃ and coefficients w̃ ∈ RM̃ . The method then approximates f˜ by a function
f ∈ H with a lower model order. Initially, this sparse approximation is the original function f = f˜ so
that its dictionary is initialized with that of the original function D = D̃, with corresponding coefficients
w = w̃. Then, the algorithm sequentially removes dictionary elements from the initial dictionary D̃,
yielding a sparse approximation f of f˜, until the error threshold kf − f˜kH ≤ t is violated, in which case
it terminates. See Appendix A for further details.
We summarize the key steps of the proposed method in Algorithm 1 for solving (3) while maintaining a
finite model order, thus allowing for the memory-efficient learning of nonparametric regression functions
online in multi-agent systems. The method, Greedy Projected Penalty Method, executes the stochastic
projection of the functional stochastic gradient iterates onto sparse subspaces HDi,t+1 stated in (21). Initial
functions are set to null fi,0 = 0, i.e., it has empty dictionary Di,0 = [] and coefficient vector wi,0 = [].
The notation [] is used to denote the empty matrix or vector respective size p × 0 or 0. Then, at each
step, given an independent training example (xi,t , yi,t ) and step-size ηt , we compute the unconstrained
functional stochastic gradient iterate (22) with respect to the instantaneous penalty function (11) which
admits the parameterization D̃i,t+1 and w̃i,t+1 as stated in (23). These parameters are then fed into KOMP
with approximation budget t , such that (fi,t+1 , Di,t+1 , wi,t+1 ) = KOMP(f˜i,t+1 , D̃i,t+1 , w̃i,t+1 , t ).
IV. C ONVERGENCE A NALYSIS
We turn to establishing that the method presented in Algorithm 1 converges with probability 1 to the
minimizer of the penalty function ψc (f ) [cf. (9)] when attenuating algorithm step-sizes are used, and to
a neighborhood of the minimizer along a subsequence when constant step-sizes are used. Moreover, for
the later case, the kernel dictionary that parameterizes the regression function fi for each agent i remains
finite in the worst case. This analysis is an application of Section IV of [31], but these results, together
with the properties of the penalty function ψc (f ) allow us to establish bounds on the deviation for each
individual in the network from the common globally optimal regression function.
Before analyzing the proposed method developed in Section III, we define key quantities to simplify
the analysis and introduce standard assumptions which are necessary to establish convergence. Define the

10

local projected stochastic functional gradient associated with the update in (21) as
˜ f ψ̂i,c (fi,t (xi,t), yi,t ) =
∇
i

h
i
fi,t − PHDi,t+1 fi,t − ηt∇fi ψ̂i,c (fi,t (xi,t),yi,t) /ηt

(25)

such that the local update of Algorithm 1 [cf. (21)] may be expressed as a stochastic descent using
˜ f ψ̂i,c (fi,t (xi,t), yi,t ) . The definitions of (25) and the local
projected functional gradients fi,t+1 = fi,t − ηt ∇
i
stochastic gradient ∇fiψ̂i,c (fi,t (xi,t), yi,t ) may be stacked to analyze the global convergence behavior of
the algorithm. For further reference, we define the stacked projected functional stochastic gradient of the
˜ f ψ̂c (ft (xt ), yt ) = [∇
˜ f1 ψ̂1,c (f1,t (x1,t), y1,t ); · · · ; ∇
˜ f ψ̂V,c (fV,t (xV,t), yV,t )]. Then the
penalty function as ∇
V
stacked global update of the algorithm is
˜ f ψ̂c (ft (xt ), yt ) .
ft+1 = ft − ηt ∇

(26)

Moreover, observe that the stochastic functional gradient in (14), based upon the fact that (xt , yt ) are
independent and identically distributed realizations of the random pair (x, y), is an unbiased estimator of
the true functional gradient of the penalty function ψc (f ) in (9), i.e.
E[∇f ψ̂c (f (xt ), yt ) Ft ] = ∇f ψc (f )

(27)

for all t. In (27), we denote as Ft the sigma algebra which measures the algorithm history for times
u ≤ t, i.e. Ft = {xu , yu , uu }tu=1 . Next, we formally state technical conditions on the loss functions, data
domain, and stochastic approximation errors that are necessary to establish convergence.
Assumption 1 The feature space X ⊂ Rp and target domain Y ⊂ R are compact, and the kernel map
may be bounded as
p
sup κ(x, x) = X < ∞
(28)
x∈X

Assumption 2 The local losses `i (fi (x), y) are convex and differentiable with respect to the first (scalar)
argument fi (x) on R for all x ∈ X and y ∈ Y. Moreover, the instantaneous losses `i : H × X × Y → R
are Ci -Lipschitz continuous for all z ∈ R for a fixed y ∈ Y
|`i (z, y) − `i (z 0 , y)| ≤ Ci |z − z 0 |

(29)

with C := maxi Ci as the largest modulus of continuity.
Assumption 3 The projected functional gradient of the instantaneous penalty function defined by stacking
(25) has finite conditional second moments:
˜ f ψ̂c (ft (xt ), yt )k2 | Ft ] ≤ σ 2
E[k∇
H

(30)

Assumption 1 holds in most settings by the data domain itself, and justifies the bounding of the loss.
Taken together, these conditions permit bounding the optimal function fc∗ in the Hilbert norm, and imply
that the worst-case model order is guaranteed to be finite. Variants of Assumption 2 appear in the analysis
of stochastic descent methods in the kernelized setting [47], [48], and is satisfied for supervised learning
problems such as logistic regression, support vector machines with the square-hinge-loss, the square loss,
among others. Moreover, it is standard in the analysis of descent methods (see [49]). Assumption 3 is
common in stochastic methods, and ensures that the stochastic approximation error has finite variance.
Next we establish a few auxiliary results needed in the proof of the main results. Specifically, we
introduce a proposition which quantifies the error due to sparse projections in terms of the ratio of the
compression budget to the learning rate.

11

Proposition 1 Given independent realizations (xt , yt ) of the random pair (x, y), the difference between
the stacked projected stochastic functional gradient and the its un-projected variant defined by (25) and
(14), respectively, is bounded as
t V
(31)
ηt
where ηt > 0 denotes the algorithm step-size and t > 0 is the approximation budget parameter of
Algorithm 2.
˜ f ψ̂c (ft (xt ), yt ) − ∇f ψ̂c (f (xt ), yt )kH ≤
k∇



Proof: See Appendix B.

With the error induced by sparse projections quantified, we may now shift focus to analyzing the
Hilbert-norm sub-optimality of the stacked iterates generated by Algorithm 1. Specifically, we have a
descent property of the sequence {ft }.
Lemma 1 (Stochastic Descent) Consider the sequence generated {ft } by Algorithm 1 with f0 = 0. Under
Assumptions 1-3, the following expected descent relation holds.


E kft+1 − fc∗ k2H Ft ≤ kft −fc∗ k2H −2ηt [ψc (ft )−ψc (fc∗ )]
+ 2t V kft − fc∗ kH +ηt2 σ 2
(32)


Proof: See Appendix B.

Now that Lemma 1 establishes a descent-like property, we may apply the proof of Theorem 1 in [31]
to kft − fc∗ kH with diminishing step-sizes. Thus we have the following corollary.
Corollary 1 Consider the sequence {ft } generated by Algorithm 1 with f0 = 0 and regularizer λ > 0.
Under Assumptions 1-3 and the hypothesis that the projection sets HDi,t in (21) are intersected with some
finite Hilbert-norm ball kf kH ≤ D for all t, with diminishing step-sizes and compression budget, i.e.,
∞
X
t=0

ηt = ∞ ,

∞
X

ηt2 < ∞ ,

t = ηt2 ,

(33)

t=0

such that ηt < 1/λ, the sequence converges exactly to the minimizer of the penalty [cf. (9)]: ft → fc∗ with
probability 1.
To attain exact convergence to the minimizer of the penalty, fc∗ , we require the compression budget
determining the error t incurred by sparse projections to approach null. This means that to have exact
convergence, we require the function representation to require an increasing amount of memory which is,
in the limit, of infinite complexity. In contrast, when constant step-size and compression budget are used,
then the algorithm settles to a neighborhood, as we state next.
Theorem 1 The sequence {ft } generated by Algorithm 1 with f0 = 0 and regularizer λ > 0, under
Assumptions 1-3, with constant step-size selection ηt = η < 1/λ and constant compression budget t =
 = Kη 3/2 for a positive constant K, converges to a neighborhood of fc∗ with probability 1:
√ h
i
√
η
√
∗
lim inf kft −fc kH ≤
KV + K 2 V 2 +λσ 2 = O( η) a.s.
(34)
t
λ
Proof: See Appendix D.

Empirically, the use of constant step-sizes has the effect of maintaining consistent algorithm adaptivity in
the face of new data, at the cost of losing exact convergence. But this drawback is more than compensated
for by the fact that in this case we may apply Theorem 3 of [31], which guarantees the model order of
the function sequence remains finite, and in the worst case, is related to the covering number of the data
domain

12

2

2

2

1

1

1

0

0

0

-1

-1

-1

-2

-2

-2
-3

-2

-1

0

1

2

(a) Gaussian Mixtures data.

3

-2

0

2

(b) Logistic Decision surface.

-3

-2

-1

0

1

2

3

(c) Hinge Decision surface.

Fig. 1: Visualizations of the Gaussian mixture data set (Figure 1a) as in [24] and the learned low-memory multi-class kernel
logistic regressor of a randomly chosen agent in the network (Figure 1b), which attains 95.2% classification accuracy on a
hold-out test set. Curved black lines denote decision boundaries between classes; dotted lines denote confidence intervals; bold
black dots denote kernel dictionary elements associated to an arbitrary i ∈ V. Kernel dictionary elements concentrate at peaks
of the Gaussian clusters and near points of overlap between classes. In Figure 1c we plot the resulting decision surface learned
by kernel SVM which attains 95.7% accuracy – the state of the art.

Corollary 2 Denote ft ∈ HV as the stacked function sequence defined by Algorithm 1 with constant
step-size ηt = η < 1/λ and approximation budget  = Kη 3/2 where K > 0 is an arbitrary positive scalar.
Let Mt be the model order of the stacked function ft i.e., the number of columns of the dictionary Dt
which parameterizes ft . Then there exists a finite upper bound M ∞ such that, for all t ≥ 0, the model
order is always bounded as Mt ≤ M ∞ .
Thus, only constant step-sizes attain a reasonable tradeoff between performance relative to fc∗ and the
complexity of storing the function sequence {ft }: in this setting, we obtain approximate convergence to
fc∗ while ensuring the memory requirements are always finite, as stated in Corollary 2.
We are left to analyze the goodness of the solution fc∗ as an approximation of the solution of the original
problem (3). In particular, we establish consensus in the mean square sense. Let us start by establishing
that the penalty term is bounded by a p∗ /c, where p∗ is the primal value of the optimization problem (3)
and c is the barrier parameter introduced in (9).
Proposition 2 Let Assumptions 1 - 3 hold. Let fc∗ be the minimizer of the penalty function (9) and let p∗
be the primal optimal value of (3). Then, it holds that
 ∗
1 XX
p∗
∗
Exi [fc,i
(xi )−fc,j
(xi )]2 ≤ .
(35)
2 i∈Vj∈n
c
i



Proof: See Appendix E.

Proposition 2 establishes a relationship between the choice of penalty parameter c and constraint
satisfaction. This result may be used to attain convergence in mean square of each individual agent’s
regression function to ones which coincide with one another. Under an additional hypothesis, we obtain
exact consensus, as we state next.
Theorem 2 Let Assumptions 1 - 3 hold. Let fc∗ be the minimizer of the penalty function (9). Then, suppose
∗
∗
the penalty parameter c in (9) approaches infinity c → ∞, and that the node-pair differences fi,c
− fj,c
are not orthogonal to mean transformation Exi [κ(xi , ·)] of the local input spaces xi for all (i, j) ∈ E.
∗
∗
Then fi,c
= fj,c
for all (i, j) ∈ E.
Proof: As a consequence, the limit of (35) when c tends to infinity yields consensus in L2 sense, i.e.,
 ∗
1 XX
∗
Exi [fc,i
(xi )−fc,j
(xi )]2 = 0,
(36)
lim
c→∞ 2
i∈V j∈n
i

13

0.6
0.4

0.2
0

500 1000 1500 2000 2500 3000 3500 4000 4500 5000

Mi,t , Model Order

Network Disagreement

Global Objective

20

10-11

1.6
1.4
1.2
1
0.8

10-12
10-13
10-14
10-15
0

1000

2000

3000

4000

5000

t, number of samples processed

t, number of samples processed

15

10

5

0
0

500 1000 1500 2000 2500 3000 3500 4000 4500 5000

t, number of samples processed

(c) Model Order Mi,t vs. samples pro(a) Global objective vs. samples processed (b) Disagreement vs. samples processed
cessed

P
Fig. 2: In Fig. 2a, we plot the global objective i∈V (Exi ,yi [`i (fi,t x), yi ]) versus P
the number of samples processed, and
observe convergence. In Fig. 2b we display the Hilbert-norm network disagreement (i,j)∈E kfi,t − fj,t k2H with a penalty
parameter c that doubles every 200 samples. As c increases, agents attain consensus. In Fig. 2c, we plot the model order of a
randomly chosen agent’s regression function, which stabilizes to 18 after 162 samples.

which, by pulling the limit outside the sum in (36), yields
 ∗
∗
lim Exi [fc,i
(xi )−fc,j
(xi )]2 = 0 ,
c→∞

(37)

for all (i, j) ∈ E. Consensus in the mean square sense is a less stringent constraint that equality in the
Hilbert norm as desired in (3). In particular, for any (i, j) ∈ E, if fi = fj , then consensus in the mean
square sense is satisfied as well. Then, apply the reproducing property of the kernel (4)(i), to write

∗
∗
−fc,j
, k(xi , ·) >
(38)
0 = lim Exi < fc,i
c→∞

∗
∗
≥ lim Exi < fc,i
−fc,j
, k(xi , ·) >
c→∞

∗
∗
= lim < fc,i
−fc,j
, Exi k(xi , ·) >
c→∞

where in the previous expression we pull the absolute value outside the expectation, and in the later we
apply linearity of the expectation. Thus, (38) implies consensus is achieved with respect to the Hilbert
∗
∗
are not orthogonal to Exi [κ(xi , ·)], the mean of the
− fc,j
norm, whenever the function differences fc,i
transformation of the local input data xi .

V. N UMERICAL E XPERIMENTS
We consider the task of kernel logistic regression (KLR) (Section V-A) from multi-class training data
scattered across a multi-agent system in two settings: classification of data from a Gaussian mixture model
and texture classification. In Section V-B, we consider kernel support vector machines (KSVM).2
A. Kernel Logistic Regression
For KLR, the merit of a particular regressor for agent i is quantified by its contribution to the classconditional probability. We define a set of class-specific functions fi,k : X → R, and denote them jointly
as fi ∈ HD , where {1, . . . , D} denotes the set of classes. Then, define the probabilistic model
exp(fi,d (xi ))
P (yi = d | xi ) := P
.
d0 exp(fi,d0 (xi ))

(39)

2
We thank Garrett Warnell and Ethan Stump of the U.S. Army Research Laboratory for invaluable assistance in the algorithm
implementation.

14

which models the odds ratio of a sample being in class d versus all others. The negative log likelihood
defined by (39) is the instantaneous loss (see, e.g., [50]) at sample (xi,n , yi,n ):
`i (fi , xi,n , yi,n ) = −log P (yi = yi,n |xi,n ).

(40)

For a given set of activation functions, classification decisions d˜ for xi is given by the maximum likelihood
estimate, i.e., d˜ = argmaxd∈{1,...,D} fi,d (x).
Gaussian Mixture Model Following [24], [31], we generate a data set from Gaussian mixture models,
which consists N = 5000 feature-label pairs for training and 2500 for testing. Each label yn was drawn
uniformly at random from the label set. The corresponding feature vector xn ∈ Rp was
P then drawn from
2
a planar (p = 2), equitably-weighted Gaussian mixture model, i.e., x y ∼ (1/3) 3j=1 N (µy,j , σy,j
I)
2
where σy,j = 0.2 for all values of y and j. The means µy,j are themselves realizations of their own
Gaussian distribution with class-dependent parameters, i.e., µy,j ∼ N (θ y , σy2 I), where {θ 1 , . . . , θ D } are
equitably spaced around the unit circle, one for each class label, and σy2 = 1.0. We fix the number of
classes D = 5, meaning that the feature distribution has, in total, 15 distinct modes. The data is plotted
in Figure 1a.
Each agent in a V = 20 network observes a unique stream of training examples from this common data
set. Here the communications graph is a random network with edges generated randomly between nodes
with probability 1/5 repeatedly until we obtain one that is connected, and then symmetrize it. We run
Algorithm 1 when the entire training set is fed to each agent in a streaming fashion, a Gaussian kernel is
used with bandwidth d = 0.6, with constant learning rate η = 3, compression budget chosen as  = η 3/2
with parsimony constant K = 0.04, mini-batch size 32, and regularizer λ = 10−6 . The penalty coefficient
is initialized as c = 0.01 and doubled after every 200 training examples.
We plotPthe results of this implementation
in Figures 1b and 2. In Figure 2a, we plot the global

objective i∈V (Exi ,yi [`i (fi,t x), yi ]) relative to the number of training examples processed, and observe
stable
convergence to a global minimum. In Figure 2b we display Hilbert-norm network disagreement
P
2
(i,j)∈E kfi,t − fj,t kH versus observed sample points. Since each regression function is initialized as null,
initially the disagreement is trivially null, but it remains small over the function sample path as model
training occurs. Moreover, the model order of an arbitrarily chosen agent i = 15 versus samples processed
is given in Figure 2c: observe that the model order stabilizes after only a couple hundred training examples
to 18, which is only a couple more than 15, the number of modes of the joint data density function. The
resulting decision surface of node 15 is given in Figure 1b, which achieves 95.2% classification accuracy
on the test set which is comparable to existing centralized batch approaches (see Table 2 of [31]) to kernel
logistic regression.
Texture Classification We generated the brodatz data set using a subset of the images provided in
[37]. Specifically, we used 13 texture images (i.e. D=13), and from them generated a set of 256 textons
[51]. Next, for each overlapping patch of size 24-pixels-by-24-pixels within these images, we took the
feature to be the associated p = 256-dimensional texton histogram. The corresponding label was given
by the index of the image from which the patch was selected. We then randomly selected N = 10000
feature-label pairs for training and 5000 for testing. Each agent in network with V = 5 observes a unique
stream of training examples from this common data set. Here the communication graph is a random
network with edges generated randomly between nodes with probability 1/5 repeatedly until we obtain
one that is connected, and then symmetrize it. To train the classifier we run Algorithm 1 ten epoches:
in each epoch we fed the entire training set to each agent in a streaming fashion. A Gaussian kernel
is used with bandwith σ 2 = 0.1, with constant learning rate η = 4, compression budget  = η 3/2 with
parsimony constant K = 0.04, mini-batch size 32 and regularizer λ = 10−5 . The penalty coefficient is set
to c = 0.02.
PWe plot the results of this experiment in Figure 3. In Figure 3a we display the global objective
convergence to
i∈V (Exi ,yi [`i (fi,t x), yi ]) relative to the number of observed examples, and observe
P
a global minimum. In Figure 3b we plot the Hilbert norm network disagreement (i,j)∈E kfi,t − fj,t k2H .
Since the initial regression function is null for all agents the disagreement is zero and as observed in Figure

15

3

5000

10 -10

2

2.5

4000
1.5

2

3000

1.5
1

2000

1
0.5

0.5
0
0

2

4

6

8

10
10 4

1000

0

0
0

2

4

6

8

10
10

0

4

5

10
10 4

(a) Global Objective vs. samples pro- (b) Disagreement vs. samples processed (c) Model Order Mi,t vs. samples processed
cessed

P
Fig. 3: In Fig. 3a, we plot the global objective i∈V (Exi ,yi [`i (fi,t x), yi ]) versus P
the number of samples processed, and
observe convergence. In Fig. 3b we display the Hilbert-norm network disagreement (i,j)∈E kfi,t − fj,t k2H with a penalty
parameter c = 0.02. In Fig. 3c, we plot the model order of a randomly chosen agent’s regression function, which stabilizes to
4299.

3b it remains small over the training. Moreover, the model order of an agent chosen at random versus
samples processed is given in Figure 3c. The resulting decission function achives 93.5% classification
accuracy over the test set which is comparable with the accuracy of the centralized version (95.6%) [31].
However the model order requiered is more than twice the model order in the centralized case (4358
in average v.s. 1833 [31]). Compared to other distributed classification algorithms the current algorithm
outperforms them. For instance D4L achieves around 75% classification accuracy [10].
B. Kernel Support Vector Machines
Now we address the problem of training a multi-class kernel support vector machine online in a multiagent systems. The merit of a particular regressor is defined by its ability to maximize its classification
margin, which may be formulated by first defining a set of class-specific activation functions fi,d : X → R,
and denote them jointly as fi ∈ HD . In Multi-KSVM, points are assigned the class label of the activation
function that yields the maximum response. KSVM is trained by taking the instantaneous loss ` to be
the multi-class hinge function which defines the margin separating hyperplane in the kernelized feature
space, i.e.,
`i (fi , xn , yn ) = max(0, 1 + fi,r (xn ) − fi,yn (xn )),

(41)

where r = argmaxd0 6=y fi,d0 (x). See [50] for further details.
We consider an implementation where each agent in a V = 20 network observes a unique stream of
training examples from the Gaussian mixtures data set (see Figure 1a). Moreover, the communications
graph is fixed as a random network with edges generated randomly between nodes with probability 1/5
repeatedly until we obtain one that is connected, and then symmetrize it. We run Algorithm 1 when the
entire training set is fed to each agent in a streaming fashion, a Gaussian kernel is used with bandwidth
σ̃ 2 = 0.6, with constant learning rate η = 3, compression budget chosen as  = η 3/2 with parsimony
constant K = 0.04, mini-batch size 32, and regularizer λ = 10−6 . The penalty coefficient is initialized as
c = 0.01 and doubled after every 200 training examples.
We plot the results
in Figures 1c and 4. In Figure 4a, we observe that the

P of this implementation
global objective
(E
[`
(f
x),
y
])
converges
stably to a global minimum as P
the number of
xi ,yi i i,t
i
i∈V
samples processed increases. In Figure 4b we display Hilbert-norm network disagreement (i,j)∈E kfi,t −
fj,t k2H versus observed sample points. Since each regression function is initialized as null, initially the
disagreement is trivially null, but it remains small over the function sample path as model training occurs,
and periodically spikes when the penalty parameter is increased. Moreover, the model order of an arbitrarily

0.3
0.2

500 1000 1500 2000 2500 3000 3500 4000 4500 5000

t, number of samples processed

Mi,t , Model Order

25

0.4

0

30

100

0.8
0.7
0.6
0.5

Network Disagreement

Global Objective

16

10-5

10

-10

10-15
0

500 1000 1500 2000 2500 3000 3500 4000 4500 5000

t, number of samples processed

20
15
10
5
0
0

500 1000 1500 2000 2500 3000 3500 4000 4500 5000

t, number of samples processed

(b) Disagreement vs. samples
(a) Global objective vs. samples
(c) Model Order Mi,t vs. samples

P
Fig. 4: In Fig. 4a, we plot the global objective i∈V (Exi ,yi [`i (fi,t x), yi ]) versus the number of samples processed, and
observe convergence,
albeit more noisily than for the differentiable logistic loss. In Fig. 4b we display the Hilbert-norm network
P
disagreement (i,j)∈E kfi,t − fj,t k2H with a penalty parameter c that doubles every 200 samples. As c increases, agents attain
consensus with respect to the Hilbert norm. In Fig. 4c, we plot the model order of a randomly chosen agent’s regression
function, which stabilizes to 22 after 354 samples. Here we obtain a slightly higher complexity classifier that achieves slightly
better accuracy.

chosen agent i = 6 versus samples processed is given in Figure 4c: the model order stabilizes after only
a couple hundred training examples to 22, which is only a couple more than 15, the number of modes
of the joint data density function. The resulting decision surface of node 6 is given in Figure 1c, which
achieves 95.7% classification accuracy, which is approximately state of the art.
VI. C ONCLUSION
In this paper, we extended the ideas in [31] to multi-agent settings with the intent of developing a
method such that a network of autonomous agents, based on their local data stream, may learn a kernelized
statistical model which is optimal with respect to information aggregated across the entire network. To do
so, we proposed an unusual penalty function whose structure is amenable to efficient parameterizations
when developing stochastic approximation-based updates. By applying functional stochastic gradient
method to this node-separable penalty combined with greedily constructed subspace projections, we
obtain a decentralized online algorithm for memory-efficient nonparametric function approximation that
is globally convergent. We obtain a controllable trade-off between optimality and memory requirements
through the design of the greedy subspace projections. Moreover, for large penalty parameter selections,
agents achieve consensus.
The empirical performance of this protocol, the Greedy Projected Penalty Method, yields state of the art
statistical accuracy for a team of interconnected agents learning from streaming data for both multi-class
kernel logistic regression and multi-class kernel support vector machines problems. These results provide
a mathematical and empirical foundation for accurate and stable multi-agent statistical inference in online
settings while preserving memory-efficiency.
A PPENDIX A: D ETAILS OF M ATCHING P URSUIT
The removal procedure is as follows: at each step, a single dictionary element j of D is selected to
be removed which contributes the least to the Hilbert-norm error minf ∈HD\{j} kf˜ − f kH of the original
function f˜, when dictionary D is used. Since at each stage the kernel dictionary is fixed, this amounts
to a computation involving weights w ∈ RM −1 only;P
that is, the error of removing dictionary point dj is
˜
computed for each j as γj = minwI\{j} ∈RM −1 kf (·) − k∈I\{j} wk κ(dk , ·)k. We use the notation wI\{j} to
denote the entries of w ∈ RM restricted to the sub-vector associated with indices I \ {j}. Then, we define
the dictionary element which contributes the least to the approximation error as j ∗ = argminj γj . If the
error incurred by removing this kernel dictionary element exceeds the given compression budget γj ∗ > t ,

17

Algorithm 2 Kernel Orthogonal Matching Pursuit (KOMP)
Require: function f˜ defined by dict. D̃ ∈ Rp×M̃ , coeffs. w̃ ∈ RM̃ , approx. budget t > 0
initialize f = f˜, dictionary D = D̃ with indices I, model order M = M̃ , coeffs. w = w̃.
while candidate dictionary is non-empty I =
6 ∅ do
for j = 1, . . . , M̃ do
Find minimal approximation error with dictionary element dj removed
X
γj =
min
kf˜(·) −
wk κ(dk , ·)kH .
wI\{j} ∈RM −1

k∈I\{j}

end for
Find index minimizing approx. error: j ∗ = argminj∈I γj
if minimal approx. error exceeds threshold γj ∗ > t
stop
else
Prune dictionary D ← DI\{j ∗ }
Revise set I ← I \ {j ∗ }, model order M ← M − 1.
Update weights w defined by current dictionary D
w = argminkf˜(·) − wT κD (·)kH
w∈RM

end
end while
return f, D, w of model order M ≤ M̃ such that kf − f˜kH ≤ t

the algorithm terminates. Otherwise, this dictionary element dj ∗ is removed, the weights w are revised
based on the pruned dictionary as w = argminw∈RM kf˜(·) − wT κD (·)kH , and the process repeats as long
as the current function approximation is defined by a nonempty dictionary. This procedure is summarized
in Algorithm 2.
A PPENDIX B: P ROOF OF P ROPOSITION 1
˜ f ψ̂c (ft (xt), yt )
Consider the square-Hilbert-norm difference of the stacked projected stochastic gradient ∇
and its un-projected variant ∇f ψ̂c (ft (xt ), yt ) defined in (25) and (14), respectively,
˜ f ψ̂c (ft (xt ), yt ) − ∇f ψ̂c (f (xt ), yt )k2H
k∇

h
i
= vec fi,t −PHDi,t+1 fi,t −ηt∇fi ψ̂i,c (fi,t (xi,t),yi,t) /ηt

 2
− vec ∇fi ψ̂i,c (fi,t (xi,t ), yi,t )
H

h
i
2
≤ V max fi,t −PHDi,t+1 fi,t −ηt∇fi ψ̂i,c (fi,t (xi,t),yi,t) /ηt

(42)

i∈V

2

− ∇fi ψ̂i,c (fi,t (xi,t ), yi,t )

H

where we apply the fact that the functional gradient is a concatenation of functional gradients associated
with each agent in (42) for the first equality, and for the second inequality we consider the worst-case
estimate across the network. Now, let’s focus on the term inside the Hilbert-norm on the right-hand side.

18

Multiply and divide ∇fi ψ̂i,c (fi,t (xi,t ), yi,t ), the last term, by ηt , and reorder terms to write
h
i

fi,t −PHDi,t+1 fi,t −ηt∇fi ψ̂i,c (fi,t (xi,t),yi,t) /ηt
2

− ∇fi ψ̂i,c (fi,t (xi,t ), yi,t )
H

1
=
fi,t −ηt ∇fi ψ̂i,c (fi,t (xi,t ), yi,t )
ηt
h
i
1
− PHDi,t+1 fi,t −ηt ψ̂i,c (fi,t (xi,t ), yi,t )
ηt
1
= 2 kf˜i,t+1 − fi,t+1 k2H
ηt

2
H

(43)

where we have substituted the definition of f˜i,t+1 and fi,t+1 in (22) and (21), respectively, and pulled
the nonnegative scalar ηt outside the norm. Now, observe that the KOMP residual stopping criterion in
Algorithm 2 is kf˜i,t+1 − fi,t+1 kH ≤ t , which we may apply to the last term on the right-hand side of
(43). This result with the inequality (42) yields (31).

A PPENDIX C: P ROOF OF L EMMA 1
Begin by considering the square of the Hilbert-norm difference between ft+1 and fc∗ = argmin ψc (f )
which minimizes (9), and expand the square to write
˜ f ψ̂c (ft (xt ), yt )k2H
kft+1 − fc∗ k2H = kft − ηt ∇
˜ f ψ̂c (ft (xt ), yt)iH
= kft −f ∗ k2 −2ηt hft −f ∗ , ∇
+

H
c
2 ˜
ηt k∇f ψ̂c (ft (xt ), yt )k2H

(44)

Add and subtract the functional stochastic gradient of the penalty function ∇f ψ̂c (ft (xt ), yt) defined in (14)
to the second term on the right-hand side of (44) to obtain
kft+1 − fc∗ k2H = kft −fc∗ k2H −2ηt hft −fc∗ ,∇f ψ̂c (ft (xt ), yt)iH
˜ f ψ̂c (ft (xt ), yt)−∇f ψ̂c (ft (xt ), yt)iH
−2ηthft −f ∗ ,∇
+

c
2 ˜
ηt k∇f `(ft (xt ), yt )k2H

(45)

We deal with the third term on the right-hand side of (45), which represents the directional error
associated with the sparse stochastic projections, by applying the Cauchy-Schwartz inequality together
with Proposition 1 to obtain
kft+1 −fc∗ k2H = kft −fc∗ k2H −2ηt hft −fc∗ ,∇f ψ̂c (ft (xt ), yt)iH
˜ f `(ft (xt ), yt )k2H
+2t V kft −fc∗ kH +ηt2 k∇
Now compute the expectation of (46) conditional on the algorithm history Ft


E kft+1 −fc∗ k2H Ft = kft −fc∗ k2H + 2t V kft − fc∗ kH + ηt2 σ 2
−2ηthft −fc∗, ∇f ψc (ft )iH

(46)

(47)

where we have applied the fact that the stochastic functional gradient in (14) is an unbiased estimator [cf.
(27)] for the functional gradient of the penalty function in (9), as well as the fact that the variance of the
functional projected stochastic gradient is finite stated in (30) (Assumption 3). Observe that since ψc (f )
is an expectation of a convex function, it is also convex, which allows us to write
ψc (ft ) − ψc (fc∗ ) ≤ hft − fc∗ , ∇f ψc (ft )iH ,

(48)

19

which we substitute into the second term on the right-hand side of the relation given in (47) to obtain


E kft+1 − fc∗ k2H Ft ≤ kft −fc∗ k2H −2ηt [ψc (ft )−ψc (fc∗ )]
(49)
+ 2t V kft − fc∗ kH + ηt2 σ 2 .
Thus the claim in Lemma 1 is valid.



A PPENDIX D: P ROOF OF T HEOREM 1
The use of the regularizer (λ/2)kf k2H in (9) implies that the penalty is λ-strongly convex in f ∈ H,
yielding
λ
(50)
kft − fc∗ k2H ≤ ψc (ft ) − ψc (fc∗ )
2
Substituting the relation (50) into the second term on the right-hand side of the expected descent relation
stated in Lemma 1, with constant step-size ηt = η and budget t = , yields
E[kft+1 − fc∗ k2H Ft ]
≤ (1 − ηλ)kft − fc∗ k2H + 2V kft − fc∗ kH + η 2 σ 2 .

(51)

The expression in (51) may be used to construct a stopping stochastic process , which tracks the suboptimality of kft − fc∗ k2H until it reaches a specific threshold, as in the proof of Theorem 2 of [31]. In doing
so, we obtain convergence
 to a neighborhood. We may define a stochastic process δt that qualifies as a
supermartingale, i.e. E δt+1 Ft ≤ δt by considering (51) and solving for the appropriate threshold by
analyzing when the following holds true
E[kft+1 − fc∗ k2H Ft ]
≤ (1 − ηλ)kft − fc∗ k2H + 2V kft − fc∗ kH + η 2 σ 2
≤ kft − fc∗ k2H .

(52)

which may be rearranged to obtain the sufficient condition
−ηλkft − fc∗ k2H + 2V kft − fc∗ kH + η 2 σ 2 ≤ 0 .

(53)

Note that (53) defines a quadratic polynomial in kft − fc∗ kH , which, using the quadratic formula, has roots
p
V
±
2 V 2 + λη 3 σ 2
kft − fc∗ kH =
(54)
λη
Observe (53) is a downward-opening polynomial in kft − fc∗ kH which is nonnegative. Thus, focus on
the positive root, substituting the approximation budget selection  = Kη 3/2 to define the radius of
convergence as
p
√ 

√
η
V + 2 V 2 +λη 3 σ 2
=
KV + K 2 V 2 +λσ 2
(55)
∆ :=
λη
λ
(55) allows us to construct a stopping process: define δt as
δt = kft − fc∗ kH
o
n
× 1 min −ηλkfu − fc∗ k2H + 2V kfu − fc∗ kH + η 2 σ 2 > ∆

(56)

u≤t

where 1{E} denotes the indicator process of event E ∈ Ft . Note that δt ≥ 0 for all t, since both kft −f ∗ kH
and the indicator function are nonnegative. The rest of the proof applies the same reasoning as that of

20

Theorem 2 in [31]: in particular, given the definition (56), either minu≤t −ηλkfu − fc∗ k2H + 2V kfu −
fc∗ kH + η 2 σ 2 > ∆ holds, in which case we may compute the square root of the condition in (52) to write
E[δt+1 Ft ] ≤ δt

(57)

Alternatively, minu≤t −ηλkfu − fc∗ k2H + 2V kfu − fc∗ kH + η 2 σ 2 ≤ ∆, in which case the indicator function
is null for all s ≥ t from the use of the minimum inside the indicator in (56). Thus in either case, (57) is
valid, implying δt converges almost surely to null, which, as a consequence we obtain the fact that either
limt→∞ kft − fc∗ kH − ∆ = 0 or the indicator function is null for large t, i.e. limt→∞ 1{minu≤t −ηλkfu −
fc∗ k2H + 2V kfu − fc∗ kH + η 2 σ 2 > ∆} = 0 almost surely. Therefore, we obtain that
√ 

√
η
∗
2
2
KV + K +λσ
a.s. ,
(58)
lim inf kft −fc kH ≤ ∆ =
t→∞
λ
as stated in Theorem 1.

A PPENDIX E: P ROOF OF P ROPOSITION 2
Let
be the minimizer of ψc (f ) defined in (9) and f ∗ be the solution of the problem (3). Since the
former is the minimizer of ψc (f ) it holds that
h
X
i λ
Exi ,yi `i (fi∗ xi ), yi + kfi∗ k2H
ψc (fc∗ ) ≤ ψc (f ∗ ) =
2
i∈V

 ∗
cX
∗
2
+
Ex [f (xi )−fj (xi )]
.
(59)
2j∈n i i
fc∗

i

Where the equality follows from the definition of ψc (f ) in (9). Since f ∗ is solution to the problem (3) it
satisfies that fi = fj for all (i, j) ∈ E, thus

Exi [fi∗ (xi ) − fj∗ (xi )]2 = 0 ,
(60)
for all (i, j) ∈ E. As a consequence, replacing ψc (fc∗ ) by its expression in the first equality in (59) and
rearranging terms yields a bound the constraint violation of fc∗ as
 ∗
1
1 XX
∗
Exi [fc,i
(xi )−fc,j
(xi )]2 ≤ (R(f ∗ ) − R(fc∗ )) ,
(61)
2 i∈Vj∈n
c
i

where R(f ) is the global regularized objective in (2), i.e.,
h

X
i λ
2
R(f ) =
Exi ,yi `i (fi xi ), yi + kfi kH .
2
i∈V

(62)

The fact that by definition p∗ = R(f ∗ ) yields (35).
R EFERENCES
[1] A. Koppel, S. Paternain, C. Richard, and A. Ribeiro, “’decentralized efficient nonparametric stochastic optimization’,” in Signal and
Information Processing (GlobalSIP), 2017 IEEE Global Conference on (to appear). IEEE, 2017.
[2] M. Anthony and P. L. Bartlett, Neural network learning: Theoretical foundations. cambridge university press, 2009.
[3] Z. Marinho, B. Boots, A. Dragan, A. Byravan, G. J. Gordon, and S. Srinivasa, “Functional gradient motion planning in reproducing
kernel hilbert spaces,” in Proceedings of Robotics: Science and Systems, Ann Arbor, MI, July 2016.
[4] R. J. Kozick and B. M. Sadler, “Source localization with distributed sensor arrays and partial spatial coherence,” IEEE Transactions
on Signal Processing, vol. 52, no. 3, pp. 601–616, 2004.
[5] A. Koppel, J. Fink, G. Warnell, E. Stump, and A. Ribeiro, “Online learning for characterizing unknown environments in ground robotic
vehicle models,” in Proc. Int. Conf. Intelligent Robots and Systems.
[6] M. Schwager, P. Dames, D. Rus, and V. Kumar, “A multi-robot control policy for information gathering in the presence of unknown
hazards,” in Robotics Research. Springer, 2017, pp. 455–472.
[7] J. Liu, Q. Chen, and H. D. Sherali, “Algorithm design for femtocell base station placement in commercial building environments,” in
INFOCOM, 2012 Proceedings IEEE. IEEE, 2012, pp. 2951–2955.

21

[8] A. Ghosh and S. Sarkar, “Pricing for profit in internet of things,” in Information Theory (ISIT), 2015 IEEE International Symposium
on. IEEE, 2015, pp. 2211–2215.
[9] A. Koppel, F. Jakubiec, and A. Ribeiro, “A saddle point algorithm for networked online convex optimization,” IEEE Trans. Signal
Process., p. 15, Oct 2015.
[10] A. Koppel, G. Warnell, E. Stump, and A. Ribeiro, “D4l: Decentralized dynamic discriminative dictionary learning,” IEEE Trans. Signal
and Info. Process. over Networks, vol. (submitted), June 2017, available at http://www.seas.upenn.edu/ aribeiro/wiki.
[11] K. Slavakis, P. Bouboulis, and S. Theodoridis, “Online learning in reproducing kernel hilbert spaces,” Signal Processing Theory and
Machine Learning, pp. 883–987, 2013.
[12] J.-B. Li, S.-C. Chu, and J.-S. Pan, Kernel Learning Algorithms for Face Recognition. Springer, 2014.
[13] S. Haykin, “Neural networks: A comprehensive foundation,” 1994.
[14] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in Advances in neural
information processing systems, 2012, pp. 1097–1105.
[15] J. Mairal, F. Bach, and J. Ponce, “Task-driven dictionary learning,” Pattern Analysis and Machine Intelligence, IEEE Transactions on,
vol. 34, no. 4, pp. 791–804, 2012.
[16] H. Robbins and S. Monro, “A stochastic approximation method,” Ann. Math. Statist., vol. 22, no. 3, pp. 400–407, 09 1951.
[17] G. Kimeldorf and G. Wahba, “Some results on tchebycheffian spline functions,” Journal of mathematical analysis and applications,
vol. 33, no. 1, pp. 82–95, 1971.
[18] B. Schölkopf, R. Herbrich, and A. J. Smola, “A generalized representer theorem,” Subseries of Lecture Notes in Computer Science
Edited by JG Carbonell and J. Siekmann, p. 416.
[19] V. Norkin and M. Keyzer, “On stochastic optimization and statistical learning in reproducing kernel hilbert spaces by support vector
machines (svm),” Informatica, vol. 20, no. 2, pp. 273–292, 2009.
[20] Y. Engel, S. Mannor, and R. Meir, “The kernel recursive least-squares algorithm,” IEEE Transactions on Signal Processing, vol. 52,
no. 8, pp. 2275–2285, Aug 2004.
[21] W. Liu, P. P. Pokharel, and J. C. Principe, “The kernel least-mean-square algorithm,” Signal Processing, IEEE Transactions on, vol. 56,
no. 2, pp. 543–554, 2008.
[22] J. Kivinen, A. J. Smola, and R. C. Williamson, “Online Learning with Kernels,” IEEE Transactions on Signal Processing, vol. 52, pp.
2165–2176, August 2004.
[23] O. Dekel, S. Shalev-Shwartz, and Y. Singer, “The forgetron: A kernel-based perceptron on a fixed budget,” in Advances in Neural
Information Processing Systems 18. MIT Press, 2006, p. 259266. [Online]. Available: http://research.microsoft.com/apps/pubs/default.
aspx?id=78226
[24] J. Zhu and T. Hastie, “Kernel Logistic Regression and the Import Vector Machine,” Journal of Computational and Graphical Statistics,
vol. 14, no. 1, pp. 185–205, 2005.
[25] B. Dai, B. Xie, N. He, Y. Liang, A. Raj, M.-F. F. Balcan, and L. Song, “Scalable kernel methods via doubly stochastic gradients,” in
Advances in Neural Information Processing Systems, 2014, pp. 3041–3049.
[26] T. Le, V. Nguyen, T. D. Nguyen, and D. Phung, “Nonparametric budgeted stochastic gradient descent,” in Artificial Intelligence and
Statistics, 2016, pp. 654–572.
[27] T. Le, T. Nguyen, V. Nguyen, and D. Phung, “Dual space gradient descent for online learning,” in Advances in Neural Information
Processing Systems, 2016, pp. 4583–4591.
[28] J. Lu, S. C. Hoi, J. Wang, P. Zhao, and Z.-Y. Liu, “Large scale online kernel learning,” Journal of Machine Learning Research, vol. 17,
no. 47, p. 1, 2016.
[29] D. Calandriello, A. Lazaric, and M. Valko, “Second-order kernel online convex optimization with adaptive sketching,” in International
Conference on Machine Learning, 2017.
[30] Y. Pati, R. Rezaiifar, and P. Krishnaprasad, “Orthogonal Matching Pursuit: Recursive Function Approximation with Applications to
Wavelet Decomposition,” in Proceedings of the Asilomar Conference on Signals, Systems and Computers, 1993.
[31] A. Koppel, G. Warnell, E. Stump, and A. Ribeiro, “Parsimonious online learning with kernels via sparse projections in function space,”
arXiv preprint arXiv:1612.04111, 2016.
[32] B. Johansson, T. Keviczky, M. Johansson, and K. Johansson, “Subgradient methods and consensus algorithms for solving convex
optimization problems,” in Proc. of the 47th IEEE Conference on Decision and Control, Cancun, Mexico, 2008, pp. 4185–4190.
[33] S. Ram, A. Nedic, and V. Veeravalli, “Distributed stochastic subgradient projection algorithms for convex optimization,” J Optimiz.
Theory App., vol. 147, no. 3, pp. 516–545, Sep. 2010.
[34] X. Nguyen, M. J. Wainwright, and M. I. Jordan, “Nonparametric decentralized detection using kernel methods,” IEEE Transactions on
Signal Processing, vol. 53, no. 11, pp. 4053–4066, 2005.
[35] P. A. Forero, A. Cano, and G. B. Giannakis, “Consensus-based distributed support vector machines,” Journal of Machine Learning
Research, vol. 11, no. May, pp. 1663–1707, 2010.
[36] P. Chainais and C. Richard, “Learning a common dictionary over a sensor network,” in Computational Advances in Multi-Sensor
Adaptive Processing (CAMSAP), 2013 IEEE 5th International Workshop on. IEEE, 2013, pp. 133–136.
[37] P. Brodatz, Textures: A Photographic Album for Artists and Designers. Dover, 1966.
[38] S. Mukherjee and S. K. Nayar, “Automatic generation of rbf networks using wavelets,” Pattern Recognition, vol. 29, no. 8, pp. 1369–
1383, 1996.
[39] S. Shalev-Shwartz, O. Shamir, N. Srebro, and K. Sridharan, “Learnability, stability and uniform convergence,” Journal of Machine
Learning Research, vol. 11, no. Oct, pp. 2635–2670, 2010.
[40] T. Evgeniou, M. Pontil, and T. Poggio, “Regularization networks and support vector machines,” Advances in computational mathematics,
vol. 13, no. 1, pp. 1–50, 2000.
[41] A. Ribeiro, “Ergodic stochastic optimization algorithms for wireless communication and networking,” IEEE Transactions on Signal
Processing, vol. 58, no. 12, pp. 6369–6386, 2010.

22

[42] K. Müller, T. Adali, K. Fukumizu, J. C. Principe, and S. Theodoridis, “Special issue on advances in kernel-based learning for
signal processing [from the guest editors],” IEEE Signal Process. Mag., vol. 30, no. 4, pp. 14–15, 2013. [Online]. Available:
http://dx.doi.org/10.1109/MSP.2013.2253031
[43] R. Wheeden, R. Wheeden, and A. Zygmund, Measure and Integral: An Introduction to Real Analysis, ser. Chapman & Hall/CRC
Pure and Applied Mathematics. Taylor & Francis, 1977. [Online]. Available: https://books.google.com/books?id=YDkDmQ hdmcC
[44] T. Suzuki, “Dual averaging and proximal gradient descent for online alternating direction multiplier method,” in Proc. 30th Int. Conf.
Machine Learning, vol. 28, no. 1, Atlanta, GA, USA, Jun. 16-21 2013, pp. 392–400.
[45] D. Needell, J. Tropp, and R. Vershynin, “Greedy signal recovery review,” in Signals, Systems and Computers, 2008 42nd Asilomar
Conference on. IEEE, 2008, pp. 1048–1050.
[46] P. Vincent and Y. Bengio, “Kernel matching pursuit,” Machine Learning, vol. 48, no. 1, pp. 165–187, 2002.
[47] M. Pontil, Y. Ying, and D. xuan Zhou, “Error analysis for online gradient descent algorithms in reproducing kernel hilbert spaces,”
Tech. Rep., 2005.
[48] Y. Ying and D. X. Zhou, “Online regularized classification algorithms,” IEEE Transactions on Information Theory, vol. 52, no. 11, pp.
4775–4788, Nov 2006.
[49] Y. Nesterov, “Introductory lectures on convex programming volume i: Basic course,” 1998.
[50] K. Murphy, Machine Learning: A Probabilistic Perspective. MIT press, 2012.
[51] T. Leung and J. Malik, “Representing and Recognizing the Visual Appearence of Materials using Three-dimensional Textons,”
International Journal of Computer Vision, vol. 43, no. 1, pp. 29–44, 1999.

