DMFSGD: A Decentralized Matrix Factorization Algorithm for Network Distance Prediction

0.97 MB·English

by  Yongjun Liao

SUBMITTEDTOIEEE/ACMTRANSACTIONSONNETWORKING 1
DMFSGD: A Decentralized Matrix Factorization
Algorithm for Network Distance Prediction
Yongjun Liao, Wei Du, Pierre Geurts and Guy Leduc
Abstract—The knowledge of end-to-end network distances is motivated numerous research on Network Coordinate System
essential to many Internet applications. As active probing of all (NCS) [8]–[12]. For instance, approaches based on Euclidean
pairwisedistancesisinfeasibleinlarge-scalenetworks,anatural
embedding have been widely studied and achieved good
idea is to measure a few pairs and to predict the other ones
performance in interesting scenarios [9], [13]. Realizing that
withoutactuallymeasuringthem.Thispaperformulatesthedis-
tance prediction problem as matrix completion where unknown the assumption of Euclidean distance properties (symmetry
2
entries of an incomplete matrix of pairwise distances are to be and triangle inequality) are often violated in practice, as
1
0 predicted. The problem is solvable because strong correlations observedinvariousstudies[9],[14]–[18],matrixfactorization
amongnetworkdistancesexistandcausetheconstructeddistance
2 has recently drawn increasing attention of the networking
matrixtobelowrank.Thenewformulationcircumventsthewell-
community [10], [11].
n known drawbacks of existing approaches based on Euclidean
a embedding. Inthispaper,weinvestigatematrixfactorizationfornetwork
J Anewalgorithm,so-calledDecentralizedMatrixFactorization distanceprediction.Inparticular,weformulatetheproblemof

5 byStochasticGradientDescent(DMFSGD),isproposedtosolve network distance prediction as a matrix completion problem
the network distance prediction problem. By letting network where a partially observed matrix is to be completed [19]–
] nodes exchange messages with each other, the algorithm is
I fully decentralized and only requires each node to collect and [21]. Here, the matrix contains distance measurements such
N
to process local measurements, with neither explicit matrix as RTTs between network nodes with some of them known
. constructions nor special nodes such as landmarks and central and the others unknown thus to be filled. Matrix completion
s
c servers.Inaddition,wecomparedcomprehensivelymatrixfactor- is only possible if matrix entries are largely correlated, which
[ ization and Euclidean embedding to demonstrate the suitability certainly holds for network distances because Internet paths
oftheformeronnetworkdistanceprediction.Wefurtherstudied
1 theincorporationofarobustlossfunctionandofnon-negativity with nearby end nodes often overlap and share common
v constraints. Extensive experiments on various publicly-available bottleneck links. These redundancies among network paths
4 datasets of network delays show not only the scalability and the cause the constructed distance matrix to be low rank, which
7 accuracy of our approach but also its usability in real Internet will be empirically demonstrated for various RTT datasets.
1 applications.
Although numerous approaches to matrix completion have
1
. IndexTerms—networkdistanceprediction,matrixcompletion, been proposed, many of which are based on low-rank matrix
1
matrix factorization, stochastic gradient descent. factorization [22]–[24], very few are directly applicable to
0
networkapplicationswheredecentralizedprocessingofdatais
2
1 I. INTRODUCTION most of the time a necessity. In this paper, we propose a fully
: decentralized algorithm based on Stochastic Gradient Descent
v On large networks such as the Internet, many applications
i (SGD),whichisfoundedonthestochasticoptimizationtheory
X require the knowledge of end-to-end network distances in
with nice convergence guarantees [25].
order to achieve Quality of Service (QoS) objectives. In the
r The so-called Decentralized Matrix Factorization by
a networking field, the distance between two network nodes is
Stochastic Gradient Descent (DMFSGD) algorithm has two
defined by the delay or latency between them, in the form
distinctfeatures.First,itrequiresneitherexplicitconstructions
of either one-way delay or more often round-trip time (RTT).
of matrices nor special nodes such as landmarks and central
Examples include peer-to-peer file sharing, overlay networks,
servers where measurements are collected and processed. In-
andcontentdistributionsystemswherepeerspreferablyaccess
stead, by letting network nodes exchange messages with each
nodes or servers that are likely to respond fast [1]–[7].
other, matrix factorization is collaboratively and iteratively
Clearly, it is infeasible to probe actively end-to-end dis-
achievedatallnodes,witheachnodeequallyretrievinganum-
tances among all pairs of nodes in large networks as the
berofdistancemeasurements.Second,thealgorithmissimple,
demandintermsofmeasurementsgrowsquadraticallywiththe
with no infrastructure, and is computationally lightweight,
scale of the network. A natural idea is to probe a small set of
containing only vector operations. These features make it
pairsandthenpredictthedistancesbetweenotherpairswhere
suitablefordealingwithpracticalproblems,whendeployedin
there are no direct measurements. This understanding has
real applications, such as measurement dynamics where net-
workmeasurements varylargelyovertime andnetworkchurn
YongjunLiaoandGuyLeducarewithResearchUnitinNetworking(RUN),
UniversityofLie`ge,Belgium.Email:{yongjun.liao,guy.leduc}@ulg.ac.be. where nodes join and leave a network frequently. Extensive
Wei Du is with Intelligent and Interactive Systems (IIS), University of experiments on various publicly-available RTT datasets show
Innsbruck,Austria.Email:[email protected].
not only the scalability and the accuracy of our approach but
PierreGeurtsiswithSystemsandModeling,UniversityofLie`ge,Belgium.
Email:[email protected]. also its usability in real Internet applications.
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 2
Our preliminary work on decentralized matrix factorization coordinate systems. Then, the ordinary nodes calculate their
for network distance prediction was published in [11]. In the coordinates with respect to the landmarks. Vivaldi extended
present paper, we make the following distinct contributions: GNP in a decentralized manner by eliminating the landmarks.
• Our previous approach in [11] was based on Alternating It simulates the network by a physical system of spring and
Least Squares (ALS), requiring each node to probe all minimizes its energy according to Hooke’s law to find an
local measurements simultaneously. In contrast, the new optimal embedding.
SGD-based approach allows each node to probe one In all metric spaces, distances undergo two important prop-
measurementatatime,makingthesystemmoreflexible. erties:
The new approach also addresses several issues arising
• Symmetry: d(A,B)=d(B,A);
when applied practically, including the difficult choice • Triangle Inequality: d(A,B)+d(B,C)(cid:62)d(A,C).
of the learning rate parameter and the consideration of
However, network distances are not necessarily symmetric
passivedistanceacquisitionsanddynamicmeasurements.
especiallywhenrepresentedbyone-waydelays[26],[27].The
• We compare comprehensively matrix factorization and
biggerissueisthepropertyoftriangleinequality.Manystudies
Euclidean embedding to reveal the suitability of matrix
haveshownthattheviolationsoftriangleinequality(TIV)are
factorization. A unified view is provided which leads to
widespread and persistent in current Internet [9], [14]–[18].
a unified optimization framework to solve both of them.
In the presence of TIVs, metric space embedding shrinks the
• Two extensions of the current matrix factorization model
long edges and stretches the short ones, degrading heavily the
areproposed,includingtheincorporationofarobustloss
accuracy of the embedding. Figure 1 illustrates the idea of
function and the introduction of constraints in the model
Euclidean embedding for network distance prediction and the
to ensure the nonnegativity of the predicted distances.
impact of TIVs on the accuracy.
These extensions are found helpful in improving the
Without loss of generality, we focus on the simplest metric
accuracy of the prediction and require little modification
space,namelyEuclideancoordinatesystems,intherestofthe
to the algorithm with no additional computational cost.
paper.
• Inaddition,moreextensiveevaluationshavebeencarried
outtostudynotonlytheimpactoftheparametersbutalso
the accuracy of our approach. In particular, we highlight
B. Matrix Factorization
the usability of our approach by simulations on real
dynamic data collected from a peer-to-peer file sharing Alternatively, matrix factorization has also been used for
application, namely Azureus [2], [13]. network distance prediction (see Figure 2 for an illustration).
The biggest advantage of matrix factorization is that it makes
The rest of the paper is organized as follows. Section II
no assumption of Euclidean distance properties and thus can
summarizes the related work on network distance prediction
tolerate the widespread TIVs and the possible asymmetry in
based on Euclidean embedding and matrix factorization. Sec-
network distance spaces.
tion III introduces the formulation of network performance
The first system based on matrix factorization was Inter-
predictionasmatrixcompletionanditsresolutionbylow-rank
net Distance Estimation Service (IDES) [10] which has the
matrix factorization. Section IV describes the decentralized
same landmark-based architecture as GNP. IDES factorizes a
matrix factorization algorithm based on Stochastic Gradient
small but full inter-landmark distance matrix, at a so-called
Descent.SectionVdiscussespossibleextensionstothecurrent
information server, by using Singular Value Decomposition
matrixfactorizationmodel.SectionVIevaluatesourapproach
(SVD). Similarly, Phoenix treated the early-entered nodes
onvariouspubliclyavailabledatasetsofRTT.Conclusionsand
as landmarks and allowed an ordinary node to select any
future work are given in Section VII.
existing nodes in the system which already have coordinates
II. RELATEDWORK assigned [28]. Landmark-based systems suffers from several
Among numerous works on network distance prediction, drawbacksincludingsingle-pointfailures,landmarkoverloads,
we only discuss and compare approaches based on Euclidean and potential security problems. The selection of landmarks
embedding and on matrix factorization due to their simplicity can also affect the accuracy of the prediction. Moreover, in
and generality. We refer the interested readers to [12] for a Section III-D, we will show that landmark-based approaches
more detailed review of this field. are actually special cases of a general decentralized matrix
factorization model and thus can also be solved by our
A. Euclidean Embedding approach.
Astraightforwardapproachtonetworkdistancepredictionis
to embed network nodes into a metric space where each node
III. NETWORKDISTANCEPREDICTIONBYMATRIX
is assigned a coordinate from which distances can be directly
FACTORIZATION
computed.TworepresentativesareGlobeNetworkPositioning
(GNP) [8] and Vivaldi [9]. This section formulates the problem of network distance
GNP firstly proposed the idea of network embedding that prediction as matrix completion and describes its resolution
relies on a small number of landmarks. Based on inter- bymatrixfactorization.Wealsoprovideaunifiedviewofdif-
landmark distance measurements, the landmarks are first em- ferent approaches to network distance prediction, the insights
bedded into a metric space such as Euclidean or spherical of which lead to a unified optimization framework.
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 3
Euclidean Space
Internet Delay Space
C(35,26)
C
Embedding
30 3
10 1 0
0
A B A(10,10) 50 B(60,10)
100
Fig.1. NetworkdistancepredictionbyEuclideanEmbedding.
D r columns Dˆ
25 20 32 23 26 39 19 35 25 28 41
25 27 20 25 31 z}|{ 24 29 6 18 6 23 30
23 25 27 33 39 29 27 19 33 13 30
20 27 20 10 19 6 24 19 6 18 34
X YT =
20 18 21 29 ≈ × 35 18 19 20 20 25 16
27 21 33 23 6 25 6 20 19 34
31 18 19 39 28 23 13 16 25 22 44
43 31 33 39 41 30 30 34 16 34 44
Fig.2. Networkdistancepredictionbymatrixfactorization.NotethatthediagonalentriesofD andDˆ areempty.
A. Problem Formulation B. Low-Rank Approximation and Matrix Factorization
Additional constraints are needed to solve the matrix com-
Assuming n nodes in the network, a n n distance matrix
× pletion problem in eq. 1. A common approach is to constrain
is constructed with some distances between nodes measured
the rank of the approximate matrix Dˆ so that
and the others unmeasured. Denote D the measured distance
matrix with dij the measured distance from node i to node Rank(Dˆ)=r, (3)
j and Dˆ the predicted distance matrix with dˆ the predicted
ij where r n for D of size n n
distance computed from some function. (cid:28) ×
The assumption in this low-rank approximation is that the
Given the above notations, network distance prediction can
entries of D are largely correlated, which causes D to have
be viewed as a matrix completion problem that estimates
a low effective rank. To show that it holds for our problem,
the missing entries in D from a small number of known
Figure 3 plots the singular values of two RTT matrices. It
entries [20]. Its resolution generally amounts to minimizing
can be seen that the singular values of both matrices decrease
a loss function of the following form
fast as the 10th singular values are 5.7% and 2.9% of the
largest ones respectively, indicating strong correlations in
n
L(D,Dˆ,W)= w l(d ,dˆ ), (1) them. The low-rank nature of many other RTT datasets have
ij ij ij
been previously reported in [29].
i,j=1
(cid:88) Directly finding Dˆ by minimizing eq. 1 subject to eq. 3 is
where W is a weight matrix with w taking values between considerably difficult due to the rank constraint. However, as
ij
0 and 1. In a simple case, w = 1 if d is measured Dˆ is of low rank, we can factorize it into the product of two
ij ij
and 0 otherwise. Note that if the distance measurements are smaller matrices, i.e.,
RTTs, then dji = dij as RTTs are approximately symmetric. Dˆ =XYT, (4)
Consequently,w =w asd andd areeitherbothknown
ji ij ji ij
or both unknown. where X and Y are of size n r. Therefore, we can get rid
l is a loss function that penalizes the difference between an of the rank constraint by repla×cing Dˆ by XYT in eq. 1, and
estimate and its desired or true value. The most commonly- then look for X and Y instead by minimizing
used loss function is the L or square loss function, n
2
L(D,X,Y,W)= w l(d ,x yT), (5)
ij ij i j
l(d,dˆ)=(d dˆ)2. (2) i,j=1
(cid:88)
−
where x is the ith row of X, y is the ith row of Y, and
i i
We will discuss other loss functions in Section V. x yT = dˆ is the estimate of d . Note that the factorization
i j ij ij
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 4
1
P2psim525
s0.8 Meridian2255
e
u
al0.6
v

r
a
ul0.4
g
n
si0.2
Fig.4. Architecturesoflandmark-based,theleftplot,anddecentralized,the
0 rightplot,systemsfornetworkdistanceprediction.Thesquaresarelandmarks
1 5 10 15 20
andthecirclesareordinarynodes.Thedirectedpathfromnodeitonodej
# singular value meansthatnodeiprobesnodej andthereforewij =1.
Fig.3. ThesingularvaluesofaRTTmatrixof2255×2255,extractedfrom
theMeridiandataset[30]andcalled“Meridian2255”,andofaRTTmatrixof of X and Y which produce the same Dˆ, the incorporation
525×525,extractedfromtheP2psimdataset[30]andcalled“P2psim525”.
of the regularization will force to choose the pair with the
Thesingularvaluesarenormalizedsothatthelargestsingularvaluesofboth
matricesareequalto1. smallest norm.
D. A Unified View of Approaches to Network Distance Pre-
in eq. 4 has no unique solution as
diction
Dˆ =XYT =XGG−1YT, (6) Although popular approaches to network distance predic-
tion vary by adopting various models including Euclidean
whereGisanyarbitraryr rinvertiblematrix.Thus,replacing
embedding and matrix factorization and by adopting different
X by XG and YT by G−×1YT results in the same Dˆ.
architectures of either landmark-based or landmark-less and
Generally, the class of techniques to solve the low-rank
thus decentralized, these seemingly different approaches all
approximation is matrix factorization. When D is complete, optimize the same function in eq. 1 but differ only in the
analytic solutions can be found by using singular value de-
setting of w and in the associated distance functions to
ij
composition (SVD) [31]. With missing entries, the factoriza- calculate dˆ .
ij
tion is usually done by iterative optimization methods such
as Gradient Descent or Newton algorithms [32]. Note that • Setting of wij: For landmark-based methods, as all paths
between landmarks are measured and ordinary nodes
additionalconstraintscanbeimposedineq.5.Forinstance,the
probe only the landmarks,
entriesofX andY canberequiredtobenonnegativeinorder
to recover a nonnegative matrix, leading to the nonnegative 1 if node j is a landmark
w = .
matrix factorization (NMF) [33]. ij (cid:40)0 otherwise
For decentralized methods, as each node equally probes
C. Incorporation of the Regularization
a number of nodes,
Matrix completion by matrix factorization suffers from a
1 if node i probes node j
well-known problem called overfitting in the field of machine w = .
ij
learning [34]. In words, directly optimizing eq. 5 often leads (cid:40)0 otherwise
to a “perfect” model with no or small errors on the training Figure 4 illustrates the architectures of landmark-based
data while having large errors on the unseen data which are and decentralized systems.
not used in learning. The problem is more severe when D is • Distance functions to calculate dˆij: For matrix factoriza-
sparse or when r is large. tion, as described above,
Acommonwaytoavoidoverfittingisthroughregularization
dˆ =x yT, (8)
that penalizes the norms of the solutions, resulting in the ij i j
following regularized loss function, For Euclidean embedding, the Euclidean distance is de-
fined as
L(D,X,Y,W,λ)= (7)
dˆ = (x x )T(x x ), (9)
n n n ij i− j i− j
wijl(dij,xiyjT)+λ xixTi +λ yiyiT, where xi and xj are(cid:113)the Euclidean coordinates of node i
i(cid:88),j=1 (cid:88)i=1 (cid:88)i=1 and node j.
whereλistheregularizationcoefficientthatcontrolstheextent Theaboveinsightssuggestaunifiedframeworktotreatand
of regularization. to solve equally network distance prediction under different
Besides avoiding overfitting, the regularization also helps modelsanddifferentarchitectures.Forinstance,thedecentral-
overcome the drift of the solutions due to the non-uniqueness izedmatrixfactorizationalgorithmsproposedinthefollowing
of the factorization (see eq. 6), which often leads to the sections can be used to solve both Euclidean embedding and
overflowsofthesolutions.Amongtheinfinitenumberofpairs landmark-based systems with little modification.
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 5
IV. DECENTRALIZEDMATRIXFACTORIZATIONFOR applications,asmeasurementscanbeacquiredondemandand
NETWORKDISTANCEPREDICTION processedlocallyateachnode.Italsohassimpleupdaterules
The goal is to find X and Y such that XYT best approxi- that involve only vector operations and is able to deal with
large-scale dynamic measurements.
matesDbyminimizingEq.7.Commonly,itrequirestocollect
1) StochasticUpdates: WhenusingSGD,eachnodeprobes
andtoprocessanumberofdistancemeasurementsatacentral
one neighbor at a time, measures its distance with respect to
node, which is a major obstacle to network applications.
that node and retrieves that node’s coordinates. Let node j be
Below, we introduce algorithms to minimize eq. 7 in a fully
the chosen neighbor by node i at the current time. Then, the
decentralized manner.
regularized losses that node i seeks to reduce with respect to
node j are
A. Problem Formulation
l = l(d ,x yT)+λx x T, (12)
A decentralized resolution of eq. 7 forbids the explicit ij ij i j i i
constructions of large matrices. Thus, each row of X and of lji = l(dji,xjyiT)+λyiyiT. (13)
Y,denotedbyx andy andcalledthexandy coordinatesof
i i The gradients of l and l are
node i in the sequel, are stored distributively at each node. To ij ji
calculate xi and yi locally, we let node i probe and exchange ∂lij = ∂l(dij,xiyjT) +λx , (14)
messages with a number of nodes in the network, called ∂x ∂x i
i i
neighbors in the sequel, and optimize the following losses ∂l ∂l(d ,x yT)
ji = ji j i +λy . (15)
n ∂y ∂y i
i i
l = wil(d ,x yT)+λx xT, (10)
i j ij i j i i In particular, the gradients of the L loss function are
2
j=1
(cid:88)
n ∂l
= (d x yT)y , (16)
li = wjil(dji,xjyiT)+λyiyiT, (11) ∂xi − ij − i j j
j=1 ∂l
(cid:88) = (d x yT)x . (17)
where wi represents the neighboring relationship of node j ∂yi − ji− j i j
j
to node i, i.e., wi = 1 if node j is a neighbor of node i or Note that we dropped the factor 2 from the derivatives of
j
equivalently if node i probes node j, and 0 otherwise. Note the regularization terms and of the L loss function for
2
that as wi is not necessarily equal to wj, measurements d mathematical convenience.
j i ij
and d may only be known by one node of either node i or Then, node i updates its coordinates along the negative
ji
node j but not by the other. gradient directions, given by
Essentially,l istheregularizedlossoftheedgesfromnode
i x = (1 ηλ)x +η(d x yT)y , (18)
itoothernodesandli isthatoftheedgesfromothernodesto i − i ij − i j j
nodei.Thus,weaddressthelarge-scaleoptimizationproblem yi = (1−ηλ)yi+η(dji−xjyiT)xj, (19)
in eq. 7 by decomposing it into a number of subproblems in
where η, called learning rate or step size, controls the speed
eqs. 10 and 11 which can be solved locally at each node by
of the updates.
using only local measurements.
2) Minibatch and Line Search: The SGD algorithm is
In seeing that eqs. 10 and 11 are standard least-squares
sensitive to the learning rate η, where a too large η results in
problems where analytical solutions exist, our previous work
largestepsofupdatesandmayoverflowthesolution,whereas
in [11] solved the matrix factorization problem by Alternat-
atoosmallη makestheconvergenceslow.Thissensitivitycan
ing Least Squares (ALS), which alternatively and iteratively
be relieved by using more training samples at the same time,
solves the small least-squares problems in eqs. 10 and 11.
leading to minibatch SGD with the following update rules
While the ALS-based algorithm performed well in simula-
n
tions on datasets containing static measurements, it requires
x =(1 ηλ)x +η wi(d x yT)y , (20)
each node to probe measurements with a number of nodes i − i j ij − i j j
j=1
simultaneously, which is impractical when deployed in real (cid:88)
n
applications. Below, we propose a different algorithm based y =(1 ηλ)y +η wi(d x yT)x . (21)
on Stochastic Gradient Descent (SGD) that processes, at each i − i j ji− j i j
j=1
node, measurements one by one and one at a time. (cid:88)
To completely get rid of η, a line search strategy can be
incorporated to determine η adaptively [35]. In particular, in
B. Stochastic Gradient Descent (SGD)
eachupdate,η startswithalargeinitialvalueandisgradually
SGD is a variation of traditional Batch Gradient Descent decreaseduntilthelossesineqs.10or11arereducedafterthe
which is often used for online machine learning [25]. Instead update. The line search algorithm for updating x is given in
i
of collecting all training samples beforehand and computing Algorithm 1. The same algorithm can be used for updating y
i
the gradients over them, each iteration of SGD chooses one by replacing eq. 10 by eq. 11 and eq. 20 by eq. 21. Note that
training sample at random and updates the parameters being δ in Line 6 is a small positive constant that helps overcome
estimated along the negative gradients computed over that poor local optimums. We will demonstrate the effectiveness
chosen sample. SGD is particularly appropriate for network of adapting η by line search using Algorithm 1 in Section VI.
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 6
Algorithm 1 Line Search (for updating x ) Algorithm 2 DMFSGD(i, j)
i
1: compute li0 by eq. 10; 1: node i retrieves dij,dji,xj,yj actively or passively;
2: initialize η with a large value; 2: node i updates the weights of its neighbors by eq. 22;
3: for i=1 to maxNumberLineSearch do 3: update xi by eq. 20 with η set by line search;
4: compute xi by eq. 20; 4: update yi by eq. 21 with η set by line search;
5: compute li by eq. 10;
6: if li <li0+δ then
7: return D. Algorithm
8: end if
We denote the SGD-based decentralized matrix factoriza-
9: η η/2;
←− tion algorithm as DMFSGD, given in Algorithm 2. Like
10: end for
Vivaldi[9],ourDMFSGDalgorithmhasnoinfrastructureand
employsthesameprocessatallnodes.Itissimple,withupdate
rules containing only vector operations.
C. Neighbor Decay and Neighbor Selection In the implementation, the coordinates of each node are
initialized with random numbers uniformly distributed be-
As mentioned earlier, it is preferable to have a system that
tween 0 and 1. Empirically, the algorithm is insensitive to
probes and processes measurements one by one. Thus, we let
the random initializations of the coordinates. We would like
each node maintain the information (distance measurements
to point out that the algorithm is one of those randomized
and coordinates) of its neighbors, i.e., the nodes with which
gossip algorithms where each node exchanges messages with
it communicates. In minibatch SGD, each node probes one
a number of other nodes randomly [36].
neighbor at a time but updates its coordinates with respect to
As mentioned earlier, the algorithm is generic and can also
allneighborsintheneighborsetusingtheirrecordedhistorical
deal with landmark-based architectures, by letting each node
information.
only select landmarks as its neighbors, and with Euclidean
A neighbor decay strategy is incorporated that scales the embedding, by adopting the Euclidean distance defined as
weight of each node in the neighbor set by its age so that eq. 9 when optimizing eq. 1, which leads to the update rule
older information receives less weight, i.e., of Vivaldi [9], given by
wji = j∈NeighbaomrSaext(−i)(aajmax−aj), (22) xi =xi−η∂l(d∂ijx,idˆij) =xi+η(dij −dˆij)xidˆ−ijxj.
(cid:80) Note that Vivaldi adopted the L2 loss function in eq. 1 with
where a is the age of the information of node j and a is
j max noregularizationincorporated,andthelearningrateη,termed
theageoftheoldestinformationintheneighborset.Notethat
differently as timestep, was adapted by taking into account
this neighbor decay strategy was firstly proposed by [13] to
someconfidencemeasureofeachnodetoitscoordinate.Thus,
overcome the problem of skewed neighbor update in Vivaldi.
VivaldicanbeviewedasaSGD-baseddecentralizedEuclidean
In words, some nodes may be probed at far greater frequency
embedding algorithm, instead of the simulation of a spring
than others due simply to their longer life cycles and a direct
system in [9].
consequence is that the optimization will become skewed
toward these nodes.
V. EXTENDEDMATRIXFACTORIZATIONMODELS
Conventionally, the neighbors of a node are selected ran-
domly and the distances between a node and its neighbors are Thissectiondiscussespossiblewaystoextendthecommon
probed by active measurements [9]. However, in practice, it matrix factorization model.
is more attractive to perform the updates of the coordinates
passively without generating any extra traffic. In some appli-
A. Robust Matrix Factorization
cations such as Azureus, passivity is enforced, as we have
no control over the selection of neighbors with which a node The widely-used L2 loss function is known to be sensitive
communicates and when it communicates with them [13]. to outliers which often occur in network measurements due
to network anomaly such as sudden traffic bursts and attacks
Therefore, we differentiate the situations where distances
from malicious nodes. Other loss functions such as the L
are probed by active and passive measurements. For the 1
loss function, the (cid:15)-insensitive loss function and the Huber
former,theconventionalrandomneighborselectionprocedure
loss function are more robust and can tolerate outliers [37],
is adopted, i.e., each node randomly selects k nodes as its
[38]. For example, the L loss function is defined as
neighbors and actively probes one of them from time to time. 1
For the latter, no neighbor selection is performed explicitly l(d,dˆ)= d dˆ. (23)
and each node maintains a relatively small set of active | − |
neighborswithwhichitrecentlycommunicatedandupdatesits Thus, we can potentially enhance the robustness of matrix
coordinates whenever a new measurement is made available. factorization by replacing the L loss function by e.g. the L
2 1
Note that this difference has no impact on the update rules in loss function, and the same SGD procedure can be applied to
eqs 18 and 19 or in eqs 20 and 21. solve the robust matrix factorization problem. Note that the
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 7
L loss function is non-differentiable and the gradients have Algorithm 3 Extended DMFSGD(i, j)
1
to be approximated by the subgradients 1, given by 1: node i retrieves dij,dji,xj,yj actively or passively;
2: node i updates the weights of its neighbors by eq. 22;
∂l
∂x =−sign(dij −xiyjT)yj, (24) 3: if use L2 loss function then
∂li 4: update xi by eq. 20 with η set by line search;
∂y =−sign(dji−xjyiT)xj. (25) 5: update yi by eq. 21 with η set by line search;
i 6: else // use L1 loss function
Replacing the gradient functions of eqs. 14 and 15 by eqs. 24 7: update xi by eq. 26 with η set by line search;
and 25, the update rules of minibatch SGD become 8: update yi by eq. 27 with η set by line search;
n 9: end if
xi =(1−ηλ)xi+η sign(dij −xiyjT)wjiyj, (26) 10: if force nonnegativity then
(cid:88)j=1 11: turn the negative entries in xi and yi into 0;
n 12: end if
y =(1 ηλ)y +η sign(d x yT)wix , (27)
i − i ji− j i j j
j=1
(cid:88)
Comparing eqs. 26 and 27 with eqs. 20 and 21, the only As distances are defined as in eq. 29, a systematic solution
difference is that for the L loss function, the magnitudes of is to factorize D by optimizing
2
the updates are proportional to the fitting errors (d xyT), n n n n
whereas for the L1 loss function, only the signs of th−e fitting wijl(dij,dˆsij)+λ xixTi +λ yiyiT. (30)
errors are taken into consideration and decide the directions i=1j=1 i=1 i=1
(cid:88)(cid:88) (cid:88) (cid:88)
of the updates. Similar SGD update rules can be derived.
B. NonNegativity Constraint D. Height Model
Conventional matrix factorization techniques do not pre- The height model in Vivaldi [9] can also be incorporated.
serve the nonnegativity of the distances. Empirically, only This model augments the x and y coordinates of a node
a very small portion of the predicted distances were found with a height. Similarly, the x and y coordinates model the
negative by our DMFSGD algorithm, and a direct solution is high-speed Internet core, while the height models the time
to turn dˆ into a small positive value if dˆ =x yT <0. packets take to travel the access link from the node to the
ij ij i j
A systematic solution is to incorporate the nonnegativity Internet core. The cause of the access link distance includes
constraint in matrix factorization, leading to the nonnegative queuing delay and low bandwidth [9]. The height augmented
matrix factorization (NMF) that optimizes symmetric distance is defined as
n n n x yT +x yT
wijl(dij,xiyjT)+λ xixTi +λ yiyiT, (28) dˆhijs = i j 2 j i +hi+hj. (31)
i(cid:88),j=1 (cid:88)i=1 (cid:88)i=1 Correspondingly, the loss function to be optimized becomes
subject to x (cid:62)0, y (cid:62)0, i=1,...,n.
i i n n n n
The optimization of NMF is not fundamentally different from wjil(dij,dˆhijs)+λ xixTi +λ yiyiT. (32)
thatoftheunconstrainedmatrixfactorization,addingonlyone (cid:88)i=1(cid:88)j=1 (cid:88)i=1 (cid:88)i=1
projectionstepthatturns thenegativeentriesinx and y into Similar SGD update rules can be derived.
i i
zeroaftereachSGDupdatewhichcausesnonoticeableimpact
on the speed of the algorithm. The technique is also known E. Extended DMFSGD Algorithm
as projected gradient descent [39]. Empirically,wefoundnoorlittleimprovementsbyincorpo-
Note that the nonnegativity constraint has been previously rating the symmetric or height-augmented symmetric distance
studied in [10], [28], both of which adopted a more heavy- function in eq. 29 or 31, thus include neither of them in
weight nonnegative least-squares solver. our system. However, the incorporation of the nonnegativity
constraint and the robust loss function not only improved
C. Symmetric Distance Matrix Factorization the accuracy but also made the results more stable and less
sensitive to parameter settings, which will be demonstrated
Also note that network distances are symmetric if repre-
in Section VI. The extended DMFSGD algorithm is given in
sented by RTT and that this symmetry is not preserved either.
Algorithm 3. Note that as the basic version in Algorithm 2 is
A direct solution is to turn the predicted distances symmetric
aspecialcaseoftheextendedversioninAlgorithm3,wewill
by defining a symmetric distance function as
refer to Algorithm 3 simply as DMFSGD in the sequel.
dˆs = dˆij +dˆji = xiyjT +xjyiT. (29)
ij 2 2 VI. EXPERIMENTSANDEVALUATIONS
In this section, we evaluate our DMFSGD algorithm2 and
1Analogously, the subgradient-based technique that optimizes non-
differentiable functions is called subgradient descent [35]. Following the compare it with state-of-the-art approaches.
conventionin[25],weusethetermSGDtorefertobothStochasticGradient
andSubGradientDescent. 2Thesourcecodeofthealgorithmwillbepubliclyavailablesoon.
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 8
TABLEI
A. Evaluation Methodology
PROPERTIESOFTHEDATASETS
Theevaluationswereperformedunderthefollowingcriteria
Dataset Nodes Symmetry TIVpercentage Dynamic
and on the following datasets.
Harvard226 226 / / Yes
1) Evaluation Criteria:
P2PSim1740 1740 Yes 85.53% No
• Cumulative Distribution of Relative Estimation Error Meridian2500 2500 Yes 96.55% No
Relative Estimation Error (REE) is defined as P2PSim525 525 Yes 76.17% No
Meridian2255 2255 Yes 96.25% No
dˆ d
ij ij Synthetic1000 1000 Yes No No
REE = | − |.
d
ij
• Stress Stress measures the overall fitness and is used to
active neighbor set for each node in the scenarios where mea-
illustrate the convergence of the algorithm, defined as
surements are obtained passively. Thus, for the Harvard226
n (d dˆ )2 dataset, we let each node maintain an active neighbor set
stress= i,j=1 ij − ij . containing the nodes it has contacted within the past 30
(cid:118)(cid:117)(cid:80) ni,j=1dij2 minutes and the timestamped measurements are processed
(cid:117)
(cid:116) in time order. For the other datasets, the random neighbor
• Median Absolute Error Me(cid:80)dian Absolute Error (MAE)
selection is used and the measurements are processed in
is defined as
randomorderwithnoneighbordecay(Line2inAlgorithm3)
MAE =median (d dˆ ). as they are static.
ij ij ij
| − | TohandlethedynamicsofthemeasurementsinHarvard226,
2) Datasets: the distance filter in [13] is adopted that smooths the streams
• Harvard226 contains dynamic and passive measure- of measurements within a moving time window, 30 minutes
ments of application-level RTTs, with timestamps, be- in this paper, by a median filter. In the evaluation, we built a
tween 226 Azureus clients collected in 4 hours [13]. static distance matrix by extracting the median values of the
• P2PSim1740 was obtained from the P2PSim project that streamsofmeasurementsbetweeneachpairofnodesandused
containsstaticRTTmeasurementsbetween1740Internet it as the ground truth.
DNS servers [40], [41].
• Meridian2500 was obtained from the Cornell Meridian B. Euclidean Embedding vs. Matrix Factorization
project that contains static RTT measurements between
Euclidean embedding and matrix factorization both solve
2500 nodes [30].
the same problem in eq. 1 but subject to different constraints.
• P2PSim525 is a complete submatrix between 525 nodes
Euclidean embedding requires Dˆ to be symmetric and to
derived from P2psim1740.
satisfy the triangle inequality, whereas matrix factorization
• Meridian2255 is a complete submatrix between 2255
only requires Dˆ to be low rank. Below, we compare em-
nodes derived from Meridian2500.
pirically Euclidean embedding and matrix factorization to
• Synthetic1000 contains the pairwise distances between
show whether this difference in constraints makes matrix
1000 nodes that are randomly generated in a 10-
factorization more suitable for network distance prediction.
dimensional Euclidean space.
1) Algorithms: To make the model comparison fair, we
Thefirstfivedatasetswereobtainedfromreal-worldnetworks
chose the state-of-the-art algorithms to solve the Euclidean
and contain a large percentage of TIV edges, whereas the last
embedding and matrix factorization problems so that both
one was synthesized and is TIV free. Here, an edge AB is
are solved to their limits. For Euclidean embedding, Multi-
claimed to be a TIV if there exists a triangle ABC where
(cid:52) Dimensional Scaling (MDS) is the most popular technique
AB > BC + AC. The last three datasets were only used
that searches the optimal embedding using an iterative algo-
in section VI-B for the purpose of comparing the models of
rithm. We adopted the MDS implementation, mdscale, in the
Euclidean embedding and matrix factorization.
statistical toolbox of matlab [42].
TableIsummarizesthesedatasets.Notethatwecanneither
Formatrixfactorization,SVDprovidestheanalyticsolution
tell the symmetry nor calculate the TIV percentage of the
which is globally optimal [31]. Generally, SVD factorizes a
Harvard226 dataset, as the measurements between network
given matrix D into three matrices of the form
nodes vary over time largely, sometimes in several orders of
magnitudes.TheHarvard226datasetisratherdensewithabout D =USVT,
3.9% pairwise paths unmeasured in 4 hours. The other paths
are measured in uneven frequencies with one measured the where U and V are unitary matrices, and S is a diagonal
maximal662timesandonetheminimal2times.About94.0% matrix with nonnegative real numbers on the diagonal. The
of the paths are measured between 40 and 60 times. positive diagonal entries are called the singular values and
3) Implementations for Different Datasets: As mentioned their number is equal to the rank of D.
earlier, the DMFSGD algorithm adopts the conventional ran- Toobtainalow-rankfactorization,wekeeponlyther large
domneighborselectionprocedureinthescenarioswheremea- singular values in S and replace the other small ones by zero.
1 1
surements are probed actively and maintains dynamically an Let S be the new S, X = US2 and YT = S2VT, where
r r r
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 9
Synthetic1000 P2PSim525 Meridian2255
0.3 0.2 0.3
SVD
0.25 MDS
0.25
0.2 0.15
s s s
s s s
e 0.15 e e 0.2
Str Str Str
0.1 0.1
0.15
0.05 SVD SVD
MDS MDS
03 5 8 10 0.053 10 50 100 0.13 10 50 100
r r r
ms)250 ms)14 ms)10
SVD SVD SVD
or (200 MDS or ( MDS or ( 9 MDS
Err Err12 Err 8
e 150 e e
ut ut10 ut 7
ol ol ol
s100 s s
b b b 6
A A A
an 50 an 8 an 5
di di di
Me 03 5 8 10 Me 63 10 50 100 Me 43 10 50 100
r r r
Fig.5. ComparisonofMDS-basedEuclideanembeddingandSVD-basedmatrixfactorizationonsynthetic1000,P2psim525andMeridian2255.Thestresses
and the median absolute errors by both methods in different dimensions/ranks are shown on the first and second rows respectively. Note that a perfect
embeddingwithnoerrorswasgeneratedforSynthetic1000inthe10dimensionalEuclideanspacebyMDS.
TABLEII
MATRIXFACTORIZATIONVS.EUCLIDEANEMBEDDING
MatrixFactorization EuclideanEmbedding
nodecoordinate xyii==((xyii11,,······,,yxiirr)) xi=(xi1,···,xir)
distancefunction dˆij =xiyjT dˆij =(cid:112)(xi−xj)T(xi−xj)
constraints lowrank Symmetry:dˆij =dˆji
TriangleInequality:dˆij <dˆik+dkj
Sr(i,i)21 = Sr(i,i). Then, Dˆ = XYT is the optimal low- thecaseswheremissingdataispresent.Wewillshowlaterthat
rank approximation to D. increasing the rank beyond some value in matrix factorization
(cid:112)
2) Evaluations: SinceSVDcannothandlemissingdata,we for a large amount of missing data will not further improve
only compare MDS and SVD on the three complete datasets the accuracy.
including Synthetic1000, P2psim525 and Meridian2255. On This comparative study reveals the model advantages of
each dataset, we ran MDS and SVD in different dimensions Matrix Factorization over Euclidean embedding. Overall, Eu-
and ranks and computed the stresses and MAE, shown in clidean embedding has a geometric interpretation which is
figure 5. It can be seen that the accuracies by SVD mono- usefulforvisualization.However,duetotheexistenceofTIVs
tonically improve on all three datasets as the rank increases, and the possible asymmetry in network distance spaces, low-
whereas consistent improvements by MDS are only found rank matrix factorization is more suitable for modeling the
on Synthetic1000 which is TIV-free. On P2psim525 and network distance spaces. Table II compares the main features
Meridian2255 where severe TIVs exist, MDS achieves no or of matrix factorization and Euclidean embedding.
little visible improvement after 10 dimensions.
These evaluations demonstrate the influences of different C. Impact of Parameters
constraintsimposedonthetwotechniques.ForEuclideanem- This section discusses and demonstrates the impact of the
bedding, the symmetry constraint doesn’t cause any problem parameters of our DMFSGD algorithm.
as the RTTs in all datasets are symmetric. However, the con- 1) k, r and λ: Our DMFSGD algorithm has two main pa-
straint of triangle inequality is strong and can’t be relieved by rameters,theregularizationcoefficientλandtherankr.Active
increasing dimensions. In contrast, matrix factorization makes probing introduces one additional parameter, the number of
no assumptions of triangle inequality, thus is not affected by neighbors k that are selected for each node. Intuitively, r is
the TIVs in the data. Note that the accuracy improvement the number of unknown variables in each coordinate and k
by increasing the rank is guaranteed for SVD-based matrix is the amount of known data that is used to estimate each
factorization. However, this conclusion cannot be extended to unknown variable in a coordinate. λ controls the extent of the
LIAOetal.:DMFSGD:ADECENTRALIZEDMATRIXFACTORIZATIONALGORITHMFORNETWORKDISTANCEPREDICTION 10
regularization which avoids both overfitting and drift of the number per node, i.e. the total number of measurements used
coordinates. by all nodes divided by the number of nodes3. It can be seen
Clearly, increasing k is equivalent to adding more data that the DMFSGD algorithm converges fast after each node
and thus always helps improve the accuracy. However, a probe,onaverage,10 k measurementsfromitsk neighbors.
×
larger k also means more probe traffic and consequently Although no k is set for Harvard226, we treat it as k =226.
higher overheads. On the other hand, only a certain number 4) Discussions: By incorporating the line search strategy,
of unknown variables can be accurately calculated from a the L loss function and the nonnegativity constraint, our
1
certainamountofknowndata.Thus,increasingrbeyondsome DMFSGD algorithm is left with two tunable parameters: the
value for a fixed k will only lead to severe overfitting and rankr andtheregularizationcoefficientλ.Thedefaultconfig-
consequently, a large λ is needed to address it. uration of λ = 1 and r = 10 is not guaranteed to be optimal
In practice, k should be fixed according to the require- in different situations and on different datasets. However,
ments of the applications by trading off between accura- fine tuning of parameters is difficult, if not impossible, for
cies and measurement overheads. Following the suggestion network applications due to the measurement dynamics and
in Vivaldi [9], we set k = 32 for P2psim1740 and for the decentralized processing where local measurements are
Meridian2500 in the rest of the paper. Note that k = 32 processedlocallyateachnodewithnocentralnodesgathering
makes the available measurements considerably sparse. For information of the entire network. Empirically, the default
instance, 32/1740 = 1.84% measurements are available for parametersettingleadstogood,thoughnotthebest,prediction
each node in P2PSim1740 and 32/2500 = 1.28% for each accuracy to a large variety of data.
nodeinMeridian2500.Recallthatnok issetforHarvard226. The setting of k = 32 has been commonly adopted in
2) Experiments under Different Configurations: We then manysystemssuchasVivaldi.However,mostsystemscontain
experimentedwithdifferentconfigurationsofr = 3,10,100 network nodes of a few thousands or less. For large systems
{ }
and λ = 0.01,0.1,1,10 , with different loss functions and of more nodes, k has to be scaled with the number of nodes
{ }
whether to incorporate the nonnegativity constraint, shown in n. According to the theory of matrix completion [19]–[21],
Figure6.η isadaptedbythelinesearch,withtheinitialvalue one can recover an unknown n n matrix of low rank r from
of 10−3 for the L loss function and of 10−2 for the L loss just about O(nrlogn) noisy e×ntries with an error which is
2 1
function. proportional to the noise level. Thus, k rlogn to guarantee
∝
In particular, we made the following observations. First, a decent prediction accuracy.
the DMFSGD algorithm is generally more accurate when the
robust L loss function and the nonnegativity constraint are
1 D. Comparisons with Vivaldi
incorporated.ThelikelyreasonsarethattheL lossfunctionis
1
Among numerous approaches on network distance pre-
insensitivetolargefittingerrorssomeofwhichareintroduced
diction, we consider Vivaldi [9] as the state of the art be-
by measurement outliers and that the nonnegativity constraint
cause of its accuracy and its practicability. To the best of
reduces the searching space which makes it easier to find
our knowledge, Vivaldi is the only system that has been
a stable solution. Thus, the robust L loss function and the
1
actually adopted in a real application, namely Azureus [2].
nonnegativity constraint are incorporated in the DMFSGD
Other approaches such as GNP [8] and IDES [9] are less
algorithm by default.
convenient due to the usage of landmarks, which makes their
Second, λ = 1 seems to be a good choice under most
application impossible in the context of passive probing of
configurationsandisthusadoptedbydefault.Third,theimpact
distancemeasurements(thusimpossibletobeevaluatedonthe
of r depends on the properties of the dataset. In Harvard226
Harvard226 dataset). Due to the insights in Section III-D, we
where available measurements are dense, the prediction ac-
consider these landmark-based systems as a special variation
curacy improves monotonically with r, whereas in the other
of a generic decentralized model.
two datasets where available measurements are sparse due to
In this paper, we only compare our DMFSGD algorithm
the setting of a small k, better performance is achieved with
r (cid:54) 10 and a large λ is needed to overcome the overfitting with Vivaldi. To address the measurement dynamics and the
skewed neighbor updates, we adopted the Vivaldi implemen-
caused by larger r’s, which confirms our analysis in the
tation in [13] 4 when dealing with the Harvard226 dataset.
previoussection.Thus,bytradingoffbetweentheperformance
The conventional Vivaldi in [9] was adopted to deal with the
on all three datasets, r =10 is adopted by default.
other two datasets. We refer to the former as Harvard Vivaldi
3) η: Asmentionedearlier,SGDissensitivetothelearning
to make the distinction. Inaddition, despite the impracticality,
rateηwhereatoolargeηleadstotheoverflowofthesolutions
wealsodemonstratetheflexibilityoftheDMFSGDalgorithm
and a too small η slows down the convergence. Although this
in dealing with the landmark-based architecture, referred to
sensitivity is reduced by minibatch SGD, it is still difficult to
as DMFSGD Landmark, by forcing each node to only select
find an appropriate constant that works for all datasets and in
allsituations.Weexperimentedwithdifferentconstantη’sand
3For P2PSim1740 and Meridian2500, at any time, the number of mea-
withthelinesearchtoadaptη dynamically.Resultsareshown
surements used by each node is statistically the same for all nodes due to
inFigure7.Itcanbeseenthatthelinesearchstrategyperforms therandomselectionsofthesourceandthetargetnodesintheupdates.For
best in terms of both accuracy and convergence speed. Note Harvard226,thisnumberissignificantlydifferentfordifferentnodesbecause
thepathswerepassivelyprobedwithunevenfrequencies.
that the convergence speed is illustrated by the stress and
4The source code was downloaded from http://www.eecs.harvard.edu/
MAE improvements with respect to the average measurement ∼syrah/nc/.