How to learn a graph from smooth signals

0.53 MB·

by  Vassilis Kalofolias

How to learn a graph from smooth signals
Vassilis Kalofolias
Signal Processing Laboratory 2 (LTS2)
Swiss Federal Institute of Technology Lausanne (EPFL), Switzerland
6
Abstract W isbig),theyareexpectedtohaveasmalldistance.
1 ij
0 Using the graph Laplacian matrix L=D−W, where
(cid:80)
2 D is the diagonal degree matrix with D = Wij,
We propose a framework that learns the ii j
n this function can be written in matrix form as
graph structure underlying a set of smooth
a
J signals. Given X ∈ Rm×n whose rows re- tr(cid:0)X(cid:62)LX(cid:1).
sideontheverticesofanunknowngraph, we
1
1 learn the edge weights w ∈Rm(m−1)/2 under The importance of the graph Laplacian has long been
+
the smoothness assumption that tr(cid:0)X(cid:62)LX(cid:1) known as a tool for embedding, manifold learning,
L] is small. We show that the problem is a clusteringandsemisupervisedlearning,seee.g. Belkin
M weighted (cid:96)-1 minimization that leads to nat- and Niyogi (2001); Belkin, Niyogi, and Sindhwani
urally sparse solutions. We point out how (2006);Zhu,Ghahramani,Lafferty,etal.(2003).More
t. known graph learning or construction tech- recently we find an abundance of methods that ex-
a
niques fall within our framework and pro- ploit this notion of smoothness to regularize various
t
s pose a new model that performs better than machine learning tasks, solving problems of the form
[
the state of the art in many settings. We
1 present efficient, scalable primal-dual based minimize g(X)+tr(cid:0)X(cid:62)LX(cid:1). (1)
v algorithmsforbothourmodelandtheprevi- X
3
ous state of the art, and evaluate their per- Zhang, Popescul, and Dom (2006) use it to en-
1
5 formance on artificial and real data. hancewebpagecategorizationwithgraphinformation,
2 Zheng et al. (2011) for graph regularized sparse cod-
0 ing. Cai et al. (2011) use the same term to regularize
. 1 INTRODUCTION
1 NMF, Jiang et al. (2013) for PCA and Kalofolias et
0 al. (2014) for matrix completion. Having good quality
6 We consider a matrix X ∈ Rm×n = [x ,...,x ](cid:62), graphs is key to the success of the above methods.
1 1 m
where each row x ∈ Rn resides on one of m nodes
: i The goal of this paper is to solve the complementary
v of a graph G. Then each of the n columns of X can
i problem of learning a good graph:
X be seen as a signal on the same graph. A simple as-
r sumption about data residing on graphs, but also the minimize tr(cid:0)X(cid:62)LX(cid:1)+f(L), (2)
a most widely used one is that it changes smoothly be- L∈L
tween connected nodes. An easy way to quantify how
where L denotes the set of valid graph Laplacians.
smooth is a set of vectors x ,...,x ∈Rn on a given
1 m
weighted undirected graph is through the function Why is this problem important? Firstly because it en-
ables us to directly learn the hidden graph structure
1(cid:88)
W (cid:107)x −x (cid:107)2, behind our data. Secondly because in most problems
2 ij i j
that can be written in the form of eq. (1), we are of-
i,j
tengivenanoisygraph, ornographatall. Therefore,
where Wij denotes the weight of the edge between startingfromtheinitialgraphandalternatingbetween
nodes i and j. In words, if two vectors xi and xj from solving problems (1) and (2) we can at the same time
a smooth set reside on two well connected nodes (so get a better quality graph and solve the task of the
initial problem.
Preliminary work. Under review by AISTATS 2016. Do Related Work.Dempster (1972) was one of the first
not distribute.
to propose the problem of finding connectivity from
Manuscript under review by AISTATS 2016
measurements,underthename“covarianceselection”. algorithm to solve our proposed model, but also the
Years later, Banerjee, El Ghaoui, and d’Aspremont onebyDongetal.Tothebestofourknowledge,these
(2008)proposedsolvingan(cid:96)-1penalizedlog-likelihood arethefirstscalablesolutionsintheliteraturetolearn
problem to estimate a sparse inverse covariance with agraphundersmoothnessassumption(2)(Section5).
unknown pattern of zeros. However, while a lot of
Toevaluateourmodel, wefirstreviewdifferentdefini-
work has been done on inverse covariance estimation,
tionsofsmoothsignalsintheliterature. Weshowhow
thelatterdifferssubstantiallyfromagraphLaplacian.
they can be unified under the notion of graph filtering
For instance, the off-diagonal elements of a Laplacian
(Section 3). We compare the models under artificial
mustbenon-positive, whileitisnotinvertiblelikethe
and real data settings. We conclude that our model is
inverse covariance.
superior in many cases and achieves better connectiv-
Wang and Zhang (2008) learn a graph with normal- ity when sparse graphs are sought (Section 6).
(cid:80)
ized degrees by minimizing the objective (cid:107)x −
i i
(cid:80)jwijxj(cid:107)2, but they assume a fixed k-NN edge pat- 2 PROPERTIES OF THE
tern. Daitch, Kelner, and Spielman (2009) considered
LAPLACIAN
the similar objective (cid:107)LX(cid:107)2 and they approximately
F
minimizeditwithagreedyalgorithmandarelaxation. Throughoutthispaperweusethecombinatorialgraph
LaplaciandefinedasL=D−W,whereD =diag(W1)
Zhangetal.(2010)alternatebetweenproblems(1)and
and 1=[1,...,1](cid:62). The space of all valid combinato-
a variation of (2). However, while they start from an
rial graph Laplacians, is by definition
initial graph Laplacian L, they finally learn a s.p.s.d.
matrix that is not necessarily a valid Laplacian. (cid:110)
L= L∈Rm×m :
The works most relevant to ours are the ones by Lake
(cid:88) (cid:111)
and Tenenbaum (2010) and by Dong et al. (2015a,b). (∀i(cid:54)=j) L =L ≤0, L =− L .
ij ji ii ij
In the first one, the authors consider a problem sim- j(cid:54)=i
ilar to the one of the inverse covariance estimation,
InordertolearnavalidgraphLaplacian,wemightbe
but impose additional constraints in order to obtain a
tempted to search in the above space, as is done e.g.
valid Laplacian. However, their final objective func-
by Dong et al. (2015b); Lake and Tenenbaum (2010).
tion contains many constraints and a computationally
We argue that it is more intuitive to search for a valid
demanding log-determinant term that makes it diffi-
weighted adjacency matrix W from the space
cult to solve. To the best of our knowledge, there is
no scalable algorithm in the literature to solve their W =(cid:8)W ∈Rm×m : W =W(cid:62), diag(W)=0(cid:9),
model. Dong et al. (2015b) propose a model that out- m +
performstheonebyLakeandTenenbaum,butstilldo leading to simplified problems. Even more, when it
not provide a scalable algorithm. This work is com- comestoactuallysolvingtheproblembyoptimization
plementary to theirs, as we not only compare against techniques, we should consider the space of all valid
theirmodel,butalsoprovideananalysisandascalable edge weights for a graph
algorithm to solve it.
(cid:110) (cid:111)
W = w ∈Rm(m−1)/2 ,
v +
Contributions. In this paper we make the link
between smoothness and sparsity. We show that so that we do not have to deal with the symmetricity
the smoothness term can be equivalently seen as a of W explicitly. The spaces L, Wm and Wv are equiv-
weighted (cid:96)-1 norm of the adjacency matrix, and mini- alent, and connected by bijective linear mappings. In
mizing it leads to naturally sparse graphs (Section 2). this paper we use Wm to analyze the problem in hand
Basedonthis,weformulateourobjectiveasaweighted and Wv when we solve the problem. Table 1 exhibits
(cid:96)-1 problem that we propose as a general framework some of the equivalent forms in the three spaces.
forsolvingproblem(2). Usingthisframeworkwepro-
poseanewmodelforlearningagraph. Weprovethat 2.1 Smooth manifold means graph sparsity
our model has effectively one parameter that controls
Letusdefinethepairwisedistancesmatrix Z ∈Rm×m:
how sparse is the learnt graph (Section 4). +
We show how our framework includes the standard Z =(cid:107)x −x (cid:107)2.
i,j i j
Gaussian kernel weight construction, but also the
modelbyDongetal.(2015b).Wesimplifytheirmodel Using this, we can rewrite the trace term as
and prove fundamental properties (Section 4).
tr(cid:0)X(cid:62)LX(cid:1)= 1tr(WZ)= 1(cid:107)W ◦Z(cid:107) , (3)
Weprovideafast,scalableandconvergentprimal-dual 2 2 1,1
Manuscript under review by AISTATS 2016
Table1: Equivalenttermsforrepresentationsfromsets containing its graph frequencies xˆ ∈R. Low frequen-
i
L,Wm,Wv. Weusez =vectorform(Z),andlinearop- cies correspond to small eigenvalues, and low-pass or
erator S that performs summation in the vector form. smooth filters correspond to decaying functions h.
Inthesequelweshowhowdifferentmodelsforsmooth
L∈L W ∈W w∈W
m v
signals in the literature can be written as smoothing
2tr(cid:0)X(cid:62)LX(cid:1) (cid:107)W ◦Z(cid:107)1,1 2w(cid:62)z problems of an initial non-smooth signal. We give an
tr(L) (cid:107)W(cid:107) 2w(cid:62)1=2(cid:107)w(cid:107)
1,1 1 example of three different filters applied on the same
– (cid:107)W(cid:107)2 2(cid:107)w(cid:107)2
F 2 signal in Figure 4 (Appendix).
diag(L) W1 Sw
1(cid:62)log(diag(L)) 1(cid:62)log(W1) 1(cid:62)log(Sw) Smooth signals by Tikhonov regularization.
(cid:107)L(cid:107)2 (cid:107)W(cid:107)2 +(cid:107)W1(cid:107)2 2(cid:107)w(cid:107)2+(cid:107)Sw(cid:107)2
F F 2 2 2 Solving problem (1) leads to smooth signals. By set-
tingg(x)= 1(cid:107)x−x (cid:107)2 wehaveaTikhonovregulariza-
where (cid:107)A(cid:107) is the elementwise norm-1 of A and ◦ is α 0
1,1 tionproblem,thatgivenanarbitraryx asinputgives
0
the Hadamard product (see Appendix). In words, the
its graph-smooth version x = (αL+I)−1x . Equiva-
0
smoothness term is a weighted (cid:96)-1 norm of W, encod-
lently, we can see this as filtering x by
0
ing weighted sparsity, that penalizes edges connecting
1
distantrowsofX. Theinterpretationisthatwhenthe h(λ)= , (6)
givendistancescomefromasmoothmanifold,thecor- 1+αλ
responding graph has a sparse set of edges, preferring where big α values result in smoother signals.
only the ones associated to small distances in Z.
Smooth signals from a probabilistic generative
Explicitly adding a sparsity term γ(cid:107)W(cid:107) to the ob- model.
1,1
jective function is a common tactic for inverse covari- Dong et al. (2015a) proposed that smooth signals can
ance estimation. However, it brings little to our prob- be generated from a colored Gaussian distribution as
(cid:16) (cid:17)
lem, as here it can be translated as merely adding a x=x¯+(cid:80) u xˆ , where xˆ ∼N 0,λ† and † denotes
i i i i i
constant to the squared distances in Z:
thepseudoinverse. Thereforexfollowsthedistribution
tr(cid:0)X(cid:62)LX(cid:1)+γ(cid:107)W(cid:107) = 1(cid:107)W ◦(2γ+Z)(cid:107) . (4) x∼N (cid:0)x¯,L†(cid:1).
1,1 2 1,1
Tosamplefromtheabove,itsufficestodrawaninitial
Note that all information of X conveyed by the trace non-smooth signal x ∼N (0,I) and then compute
0
term is contained in the pairwise distances matrix Z,
x=x¯+h(L)x , (7)
sothattheoriginalcouldbeomitted. Moreover,using 0
√
the last term of eq. (3) instead of the trace enables us
with h(L)= L†, or equivalently filter it by
todefineotherkindsofdistancesinsteadofEuclidean.
(cid:40)√
λ−1 ,λ>0
Note finally that the separate rows of X do not have h(λ)= (8)
0 ,λ=0
to be smooth signals in some sense. Two non-smooth
signalsx ,x canhaveasmalldistancebetweenthem,
i j andaddthemeanx¯. Wepointoutherethatusingeq.
and therefore a small entry Z .
i,j (7) on any x ∼N (0,I) and for any filter h(λ) would
0
yield samples from
3 WHAT IS A SMOOTH SIGNAL?
x∼N (cid:0)x¯,h(L)2(cid:1),
Givenagraph,differentdefinitionsofwhatisasmooth
therefore the probabilistic generative model can be
signal have been used in different contexts. In this
used for any filter h. However, it does not cover cases
section we unify these different definitions using the
where the initial x is not white Gaussian.
0
notion of filtering on graphs. For more information
aboutsignalprocessingongraphswerefertothework Smooth signals by heat diffusion on graphs
of Shuman et al. (2013). Filtering of a graph signal Anothertypeofsmoothsignalsintheliteratureresults
x∈Rm by a filter h(λ) is defined as the operation1 from the process of heat diffusion on graphs. See for
example the work by Zhang and Hancock (2008) for
(cid:88) (cid:88)
y =h(L)x= uih(λi)u(cid:62)i x= uih(λi)xˆi, (5) an application on image denoising by heat diffusion
i i smoothing on the pixels graph. Given an initial signal
where {ui,λi} are eigenvector-eigenvalue pairs of L, x0, the result of the heat diffusion on a graph after
and xˆ ∈ Rm is the graph Fourier representation of x time t is x=exp(−Lt)x0, therefore the corresponding
filter is
1We denote by h both the function h : R → R and
h(λ)=exp(−tλ), (9)
its matrix counterpart h : Rm×m → Rm×m acting on the
matrix’s eigenvalues. where bigger values of t result in smoother signals.
Manuscript under review by AISTATS 2016
4 LEARNING A GRAPH FROM graphs, we want to make sure that each node has at
SMOOTH SIGNALS least one edge with another node. It is also desirable
to have control of how sparse is the resulting graph.
Inordertolearnagraphfromsmoothsignals,wepro- To meet these expectations, we propose the following
pose, asexplainedinSection2, torewriteproblem(2) model with parameters α > 0 and β ≥ 0 controlling
using the weighted adjacency matrix W and the pair- the shape of the edges:
wise distance matrix Z instead of X:
minimize (cid:107)W ◦Z(cid:107) − α1(cid:62)log(W1) + β(cid:107)W(cid:107)2.
minimize (cid:107)W ◦Z(cid:107)1,1 + f(W). (10) W∈Wm 1,1 F
W∈Wm (12)
Since W is positive we could replace the first term The logarithmic barrier acts on the node degree vec-
by tr(WZ), but we prefer this notation to keep in tor W1, unlike the model of Proposition 1 that has a
mind that our problem already has a sparsity term on similar barrier on the edges. This means that it forces
W. This means that f(W) has to play two important the degrees to be positive, but does not prevent edges
roles: (1) prevent W from going to the trivial solution from becoming zero. This improves the overall con-
W = 0 and (2) impose further structure using prior nectivityofthegraph,withoutcompromisingsparsity.
information on W. This said, depending on f the Note however, that adding solely a logarithmic term
solutionisexpectedtobesparse,thatisimportantfor (β = 0) leads to very sparse graphs, and changing α
large scale applications. onlychangesthescaleofthesolutionandnotthespar-
sitypattern(Proposition2forβ =0). Forthisreason,
Inordertomotivatethisgeneralgraphlearningframe-
we add the third term.
work, we show that the most standard weight con-
struction,aswellasthestateoftheartgraphlearning We showed in eq. (4) that adding an (cid:96)-1 norm to
model are special cases thereof. controlsparsityisnotveryuseful. Ontheotherhand,
addingaFrobeniusnormwepenalizetheformationof
4.1 Classic Laplacian computations big edges but do not penalize smaller ones. This leads
tomoredenseedgepatternsforbiggervaluesofβ. An
In the literature one of the most common practices is
interesting property of our model is that even if it has
to construct edge weights given X from the Gaussian
two terms shaping the weights, if we fix the scale we
function
(cid:18) (cid:107)x −x (cid:107)2(cid:19) then need to search for only one parameter:
w =exp − i j 2 . (11)
ij 2σ2 Proposition 2. Let F(Z,α,β) denote the solution of
our model (12) for input distances Z and parameters
It turns out that this choice of weights can be seen as
α, β. Then the following property holds for any γ >0:
the result of solving problem (10) with a specific prior
on the weights W: (cid:18) (cid:19)
α
F (Z,α,β)=γF Z, ,βγ =αF (Z,1,αβ). (13)
Proposition 1. The solution of the problem γ
(cid:88)
minimize (cid:107)W ◦Z(cid:107) + 2σ2 W (log(W )−1)
1,1 ij ij
W∈Wm ij
Proof. See appendix.
is given by eq. (11).
Proof. The problem is edge separable and the objec-
tivecanbewrittenas(cid:80) W Z +2σ2W (log(W − This means that for example if we want to obtain a
i,j ij ij ij ij
W with a fixed scale (cid:107)W(cid:107) = s (for any norm), we
1)). Deriving w.r.t. W we obtain the optimal-
ij
ity condition Z + 2σ2log(W ) = 0, or W = can solve the problem with α = 1, search only for
ij ij ij
exp(−Z /(2σ2)), that proves the proposition. a parameter β that gives the desired edge density and
ij
thennormalizethegraphbythenormwehavechosen.
Notethathere,thelogarithminf preventstheweights
The main advantage of our model over the method by
from going to 0, leading to full matrices, and sparsifi-
Dong et al. (2015a), is that it promotes connectivity
cation has to be imposed explicitly afterwards.
by putting a log barrier directly on the node degrees.
Even for β = 0, we obtain the sparsest solution pos-
4.2 Our proposed model
sible, that assigns at least one edge to each node. In
Basedonourframework(10)ourgoalistogiveagen- this case, the distant nodes will have smaller degrees
eral purpose model for learning graphs, when no prior (because of the first term), but still be connected to
informationisavailable. Inordertoobtainmeaningful their closest neighbour similarly to a 1-NN graph.
Manuscript under review by AISTATS 2016
4.3 Fitting the state of the art in our choices of f(W). We use primal dual techniques that
framework scale, like the ones reviewed by Komodakis and Pes-
quet (2014) to solve the two state of the art models:
Dong et al. (2015a) proposed the following model for
theoneweproposeandtheonebyDongetal.(2015b).
learning a graph:
Using these as examples, it is easy to solve many in-
minimize tr(cid:0)X(cid:62)LX(cid:1)+α(cid:107)L(cid:107)2, teresting models from the general framework (10).
F
L∈L
Inordertomakeoptimizationeasier,weusethevector
s.t., tr(L)=s.
form representation from space W (see Table 1), so
v
that the symmetricity does not have to be imposed as
Parameter s > 0 controls the scale (Dong et al. set it
a constraint. We write the problem as a sum of three
tom),andparameterα≥0controlsthedensityofthe
functions in order to fit it to primal dual algorithms
solution. This formulation has two weaknesses. First,
reviewed by Komodakis and Pesquet (2014). The gen-
usingaFrobeniusnormontheLaplacianhasareduced
eral form of our objective is
interpretability: the elements of L are not only of dif-
ferent scales, but also linearly dependent. Secondly,
minimize f (w)+f (Kw)+f (w), (17)
optimizing it is difficult as it has 4 constraints on L: 1 2 3
w∈Wv
3 in order to constrain L in space L, and one to keep
the trace constant. We propose to solve their model where f1 and f2 are functions for which we can ef-
using our framework: Using transformations of Table ficiently compute proximal operators, and f3 is dif-
1, we obtain the equivalent simplified model ferentiable with gradient that has Lipschitz constant
ζ ∈ (0,∞). K is a linear operator, so f is defined
2
minimize (cid:107)W ◦Z(cid:107)1,1+α(cid:107)W1(cid:107)2+α(cid:107)W(cid:107)2F, on the dual variable Kw. In the sequel we explain
W∈Wm how this general optimization framework can be ap-
s.t., (cid:107)W(cid:107) =s. (14)
1,1 plied to the two models of interest, leaving the details
intheAppendix. Forabetterunderstandingofprimal
Using this parametrization, solving the problem be-
dualoptimizationorproximalsplittingmethodswere-
comes much simpler, as we show in Section 5. Note
fer the reader to the works of Combettes and Pesquet
that for α = 0 we have a linear program that assigns
(2011); Komodakis and Pesquet (2014).
weightstotheedgecorrespondingtothesmallestpair-
wise distance in Z, and zero everywhere else. On the In our model, the second term acts on the degrees
other hand, setting α to big values, we penalize big of the nodes, that are a linear function of the edge
degrees (through the second term), and in the limit weights. Therefore we use K = S, where S is the
α → ∞ we obtain a dense graph with constant de- linear operator that satisfies W1 = Sw if w is the
greesacrossnodes. Wecanalsoprovesomeinteresting vectorform of W. In the first term we group the posi-
properties of (14): tivity constraint of W and the weighted (cid:96)-1, and the
v
Proposition 3. Let H(Z,α,s) denote the solution of second and third terms are the priors for the degrees
model (14) for input distances Z and parameters α andtheedgesrespectively. Inordertosolveourmodel
and s. Then for γ >0 the following properties hold: we define
f (w)=1{w ≥0}+2w(cid:62)z,
1
H(Z+γ,α,s)=H(Z,α,s) (15)
f (d)=−α1(cid:62)log(d),
(cid:18) (cid:19) 2
s
H(Z,α,s)=γH Z,αγ, =sH(Z,αs,1) (16) f (w)=β(cid:107)w(cid:107)2, with ζ =2β,
γ 3
where 1{} is the indicator function that becomes zero
Proof. See appendix.
when the condition in the brackets is satisfied, infinite
In other words, model (14) is invariant to adding any otherwise. Note that the second function f is defined
2
constant to the squared distances. The second prop- on the dual variable d = Sw ∈ Rm, that here is very
ertymeansthatsimilarlytoourmodel,thescaleofthe conveniently the vector of the node degrees.
solutiondoesnotchangetheshapeoftheconnectivity.
For model (14) we can define in a similar way
If we fix the scale to s, we obtain the whole range of
edge shapes given by H only by changing parameter
f (w)=1{w ≥0}+2w(cid:62)z,
α. 1
f (c)=1{c=s},
2
5 OPTIMIZATION f (w)=α(cid:0)2(cid:107)w(cid:107)2+(cid:107)Sw(cid:107)2(cid:1), with ζ =2α(m+1),
3
Anadvantageofusingtheformulationofproblem(10) and use K = 21(cid:62) so that the dual variable is c =
is that it can be solved efficiently for a wide range of Kw =(cid:107)W(cid:107) , constrained by f to be equal to s.
1,1 2
Manuscript under review by AISTATS 2016
Algorithm 1 Primal dual algorithm for model (12). 1. Random Geometric Graph (RGG): We sample x
uniformlyfrom[0,1]2 andconnectnodesusingeq.
1: Input: z,α,β, w0 ∈W , d0 ∈Rm, γ, tolerance (cid:15)
v + (11) with σ =0.2, then threshold weights <0.6.
2: for i=1,...,i do
max
3: yi =wi−γ(2βwi+S(cid:62)di) 2. Non-uniform: We sample x in [0,1]×[0,5] from
4: y¯i =di+γ(Swi) a non-uniform distribution p ∝ 1/(1+αx )
x1,x2 2
5: pi =max(0,yi−2γz) and connect nodes using eq. (11) with σ = 0.2.
6: p¯i =(y¯i−(cid:112)(y¯i)2+4αγ)/2 (cid:46) elementwise We threshold weights smaller than the best con-
7: qi =pi−γ(2βpi+S(cid:62)pi) nection of the most distant node (≈0.01).
8: q¯i =p¯i+γ(Spi) 3. Erd˝os R´enyi: Random graph as proposed by
9: wi =wi−yi+pi; Gilbert (1959) (p=3/m).
10: di =di−y¯i+q¯i;
4. Barab´asi-Albert: Random scale-free graph with
11: if (cid:107)wi−wi−1(cid:107)/(cid:107)wi−1(cid:107)<(cid:15) and
12: (cid:107)di−di−1(cid:107)/(cid:107)di−1(cid:107)<(cid:15) then preferential attachment as proposed by Barab´asi
13: break and Albert (1999) (m0=1, m=2).
14: end if
Signal Types.To create a smooth signal we filter a
15: end for
Gaussian i.i.d. x by eq. (5), using one of the three
0
filter types of Section 3. We normalize the Laplacian
Using these functions, the final algorithm for our
((cid:107)L(cid:107) =1) so that the filters g(λ) are defined for λ∈
2
model is given as Algorithm 1, and for the model by
[0,1]. See Table 3 (Appendix) for a summary.
Dong et al. as Algorithm 2 in the Appendix. Vector
z ∈Rm×(m−1)/2 is the vector form of Z, and parame- 1. Tikhonov: g(λ)= 1 as in eq. (6).
+ 1+10λ √
ter γ ∈(0, 1+ζ+(cid:107)K(cid:107)) is the stepsize. 2. Generative Model: g(λ) = 1/ λ if λ > 0, g(0) =
0 from model of eq. (8) (x¯=0).
5.1 Complexity and Convergence 3. Heat Diffusion: g(λ)=exp(−10λ) as eq. (9).
For all cases we use m=100 nodes, smooth signals of
Both algorithms that we propose have a complex-
lengthn=1000, andadd10%((cid:96)-2sense)noisebefore
ity of O(m2) per iteration, for m nodes graphs, and
computingpairwisedistances. Weperformgridsearch
they can easily be parallelized. As the objective func-
tofindthebestparametersforeachmodel. Werepeat
tions of both models are proper, convex, and lower-
the experiment 20 times for each case and report the
semicontinuous,ouralgorithmsareguaranteedtocon-
average result of the parameter value that performs
vergetotheminimum(KomodakisandPesquet2014).
best for each of the different metrics.
Metrics.Since we have the ground truth graphs for
6 EXPERIMENTS
each case, we can measure directly the relative edge
error in the (cid:96)-1 and (cid:96)-2 sense. We also report the rel-
We compare our model against the state of the art
(cid:80)
ative error of the weighted degrees d = W . This
modelbyDongetal.(2015a)solvedbyourAlgorithm i j ij
is important because both models are based on priors
2 for both artificial and real data. Comparing to the
onthedegreesasweshowinsection4. Wealsoreport
modelbyLakeandTenenbaum(2010)wasnotpossible
the F-measure (harmonic mean of edge precision and
evenforthesmallgraphsofourartificialexperiments,
recall),thatonlytakesintoaccountthebinarypattern
as there is no scalable algorithm in the literature and
of existing edges and not the weights.
the use of CVX with the log-determinant term is pro-
hibitive. Other models based on the log-det term, for Baselines. The baseline for the relative errors is a
which scalable algorithms exist, are irrelevant to our classic graph construction using equation (11) with a
problem as a sparse inverse covariance is not a valid grid search for the best σ. Note that this exact equa-
Laplacian and are known to not perform well for our tionwasusedtocreatethetwofirstartificialdatasets.
setting (see Dong et al. (2015a) for a comparison). However, using a fully connected graph with the F-
measure does not make sense. For this metric the
6.1 Artificial data baselineissettothebestedgepatternfoundbythresh-
olding (11) with different thresholds.
Thedifficultyofsolvingproblem(10)dependsbothon
Table 2 summarizes all the results for different combi-
the quality of the graph behind the data and on the
nationsofgraphs/signals. Inmostofthem, ourmodel
type of smoothness of the signals. We test 4 different
performs better for all metrics. We can see that the
types of graphs using 3 different types of signals.
signals constructed following the generative model (7)
Graph Types. We use two 2-D manifold based do not yield better results in terms of graph recon-
graphs,oneuniformlyandonenon-uniformlysampled, struction. Using smoother “Tikhonov” signals from
and two graphs that are not manifold structured: eq. (6) or “Heat Diffusion” signals from (9) by set-
Manuscript under review by AISTATS 2016
Table 2: Performance of Different Models on Artificial Data.
Tikhonov Generative Model Heat Diffusion
base Dong etal Ours base Dong etal Ours base Dong etal Ours
Rand. Geometric
F-measure 0.685 0.885 0.913 0.686 0.877 0.909 0.758 0.837 0.849
edge (cid:96)-1 0.866 0.357 0.298 0.798 0.371 0.348 0.609 0.524 0.447
edge (cid:96)-2 0.676 0.376 0.336 0.658 0.397 0.390 0.576 0.531 0.468
degree (cid:96)-1 0.142 0.146 0.065 0.261 0.147 0.112 0.209 0.227 0.142
degree (cid:96)-2 0.708 0.172 0.079 0.689 0.174 0.128 0.474 0.264 0.176
Non Uniform
F-measure 0.686 0.863 0.858 0.633 0.840 0.832 0.766 0.839 0.830
edge (cid:96)-1 0.821 0.423 0.349 0.864 0.487 0.472 0.594 0.565 0.473
edge (cid:96)-2 0.706 0.434 0.344 0.735 0.480 0.474 0.550 0.587 0.451
degree (cid:96)-1 0.160 0.184 0.055 0.235 0.185 0.100 0.233 0.255 0.128
degree (cid:96)-2 0.612 0.209 0.073 0.632 0.215 0.161 0.427 0.324 0.157
Erdo˝s R´enyi
F-measure 0.288 0.766 0.893 0.199 0.755 0.896 0.377 0.629 0.655
edge (cid:96)-1 1.465 0.448 0.391 1.566 0.478 0.427 1.379 0.832 0.841
edge (cid:96)-2 1.060 0.442 0.402 1.105 0.457 0.440 1.033 0.735 0.726
degree (cid:96)-1 0.094 0.107 0.046 0.099 0.105 0.066 0.182 0.179 0.183
degree (cid:96)-2 0.986 0.161 0.066 1.312 0.181 0.151 0.892 0.236 0.273
Barab´asi-Albert
F-measure 0.345 0.710 0.868 0.382 0.739 0.838 0.352 0.690 0.765
edge (cid:96)-1 1.531 0.614 0.533 1.496 0.652 0.624 1.468 0.740 0.675
edge (cid:96)-2 1.061 0.568 0.506 1.036 0.611 0.571 1.041 0.662 0.590
degree (cid:96)-1 0.175 0.264 0.111 0.199 0.264 0.207 0.254 0.317 0.148
degree (cid:96)-2 0.554 0.340 0.201 0.556 0.333 0.287 0.568 0.414 0.283
tingλ=20yieldedslightlyworseresultsinbothcases For each of the graphs, we run standard spectral clus-
(not reported here). It also seems that the results are tering(asintheworkofNg,Jordan,Weiss,etal.2002
slightlybetterforthemanifoldrelatedgraphsthanfor but without normalizing the Laplacian) with k-means
theErd˝osR´enyiandBarab´asi-Albertmodels,aneffect 100 times. We also run label propagation we choose
that is more prevalent when we use signals of length 100 times a different subset of 10% known labels.
n = 100 smooth signals instead of 1000 (c.f. Table 4
In Fig. 1 we plot the behavior of different models
ofAppendix). Thiswouldbeinterestingtoinvestigate
for different density levels. The horizontal axis is the
theoretically.
average number of non-zero edges per node. In the
left plot we see the clustering quality. Even though
6.2 Real data
the best result of both algorithms is almost the same
Wealsoevaluatetheperformanceofourmodelonreal (0.24vs0.25),ourmodelismorerobustintermsofthe
data. Inthiscase,theactualgroundtruthgraphisnot graph density choice. A similar behavior is exhibited
known. We therefore measure the performance of dif- forlabelpropagationplottedinthemiddle. Theclas-
ferentmodelsonspectralclusteringandlabelpropaga- sificationqualityisbetterforourmodelinthesparser
tion, two algorithms that depend solely on the graph. graph density levels.
Note that an explicit Laplacian normalization is not
The robustness of our model for small graph densi-
needed for the learned models (it is even harmful as
ties can be explained by the connectivity quality plot-
foundexperimentally),sincethisroleisalreadyplayed
ted in the right. The continuous lines are the num-
by the regularization.
ber of different connected components in the learned
Learning the graph of USPS digits graphs,thatisameasureofconnectivity: thelesscom-
We first learn the graph connecting 1001 different im- ponents there are, the better connected is the graph.
ages of the USPS dataset, that are images of digits The dashed blue line is the number of disconnected
from 0 to 9 (10 classes). We follow Zhu, Ghahra- nodes of model Dong et al. 2015a. The latter fails to
mani, Lafferty, et al. (2003) and sample the class sizes assign connections to the most distant nodes, unless
non-uniformly. For each class i ∈ {1...10} we take the density of the graph reaches a fairly high level. If
round(2.6i2) images, resulting to classes with sizes we want a graph with 6 edges per node, our model re-
from 3 to 260 images each. We learn graphs of dif- turnsagraphwith3componentsandnodisconnected
ferentdensitiesusingbothmodels. Asbaselineweuse nodes. ThemodelbyDongetal.returnsagraphwith
a k-Nearest Neighbors (k-NN) graph for different k. 35componentsoutofwhich22aredisconnectednodes.
Manuscript under review by AISTATS 2016
Spectral Clustering Label Propagation Connectivity
0.45 0.35 70
[Dong etal] miss rate (Dong etal) # disconnected (Dong etal)
Ours miss rate (Ours) 60 # components (Dong etal)
0.4 k-NN 0.3 miss rate (k-NN) # components (Ours)
rorre gniretsulC0.03.53 rorre noitacifissa0.02.52 23450000 #n ucmombepro onfe cnlatss s(ke-sNN)
lC
0.25 0.15
10
0.2 0.1 0
5 10 15 20 5 10 15 20 5 10 15 20
graph density graph density graph density
Figure1: Graphlearnedfrom1001USPSimages. Left: Clusteringquality. Middle: Labelpropagationquality.
Right: Number of completely disconnected components (continuous lines) and number of disconnected nodes
for model by Dong etal. (blue dashed line). Our model and k-NN have no disconnected nodes.
0.4 0.4 0.4
miss rate (Dong etal) miss rate (Dong etal) miss rate (Dong etal)
0.35 miss rate (Ours) 0.35 miss rate (Ours) 0.35 miss rate (Ours)
miss rate (k-NN) miss rate (k-NN) miss rate (k-NN)
ro 0.3 unclassifiable (Dong etal) ro 0.3 unclassifiable (Dong etal) ro 0.3 unclassifiable (Dong etal)
rre n0.25 uunnccllaassssiiffiiaabbllee ((Ok-uNrNs)) rre n0.25 uunnccllaassssiiffiiaabbllee ((Ok-uNrNs)) rre n0.25 uunnccllaassssiiffiiaabbllee ((Ok-uNrNs))
o o o
ita 0.2 ita 0.2 ita 0.2
c c c
ifis0.15 ifis0.15 ifis0.15
s s s
a a a
lC 0.1 lC 0.1 lC 0.1
0.05 0.05 0.05
0 0 0
2 4 6 8 10 12 14 16 18 2 4 6 8 10 12 14 16 18 2 4 6 8 10 12 14 16 18
graph density graph density graph density
Figure 2: Label propagation for the problem “1” vs. “2” of MNIST with different class size proportions: 1 to 4
(left), 1 to 1 (middle) or 4 to 1 (right). Missclassification rate for different number of edges per node.
Note that in real applications where the best density etal.2015afailstorecoveredgesbetweendifferentdig-
level is not known a priori, it is important for a graph its“2”unlessthereturnedgraphisfairlydense,unlike
learning model to perform well for sparse levels. This ourmodelthatevenforverysparsegraphlevelstreats
isespeciallythecaseforlargescaleapplications,where the different classes more fairly. The effect is stronger
more edges mean more computations. when the set of 2’s is also the smallest of the two.
Time:Algorithm1implementedinMatlab2 learneda
7 CONCLUSION
10-edge/nodegraphof1001USPSimagesin5seconds
(218 iterations) and Algorithm 2 in 1 minute (2043
We introduce a new way of addressing the prob-
iterations) on a standard PC for tolerance (cid:15)=1e-4.
lem of learning a graph under the assumption that
Learning the graph of MNIST 1 vs 2 tr(cid:0)X(cid:62)LX(cid:1) is small. We show how the problem can
be simplified into a weighted sparsity problem, that
To demonstrate the different behaviour of the two
implies a general framework for learning a graph. We
models for non-uniform sampling cases, we use the
show how the standard Gaussian weight construction
problem of classification between digits 1 and 2 of the
from distances is a special case of this framework. We
MNIST dataset. This problem is particular because
proposeanewmodelforlearningagraph,andprovide
digits “1” are close to each other (average square dis-
ananalysisofthestateoftheartmodelofDongetal.
tance of 45), while digits “2” differ more from each
(2015a) that also fits our framework. The new formu-
other (average square distance of 102). In Figure 2
lationenablesustoproposeafastandscalableprimal
we report the average miss-classification rate for dif-
dual algorithm for our model, but also for the one of
ferent class size proportions, with 40 1’s and 160 2’s
Dongetal.2015athatwasmissingfromtheliterature.
(left), 100 1’s and 100 2’s (middle) or 160 1’s and
Our experiments suggest that when sparse graphs are
40 2’s (right). Results are averaged over 40 random
tobelearned, butconnectivityiscrucial, ourmodelis
draws. The dashed lines denote the number of nodes
expected to outperform the current state of the art.
contained in components without labeled nodes, that
can not be classified. In this case, the model of Dong We hope not only that our solution will be used for
many applications that require good quality graphs,
2Code for both models is available as part of the open-
source toolbox GSPBox by Perraudin et al. (2014b) using but also that our framework will trigger defining new
code from UNLocBoX, Perraudin et al. (2014a). graph learning models targeting specific applications.
Manuscript under review by AISTATS 2016
Acknowledgements Kalofolias, Vassilis et al. (2014). “Matrix completion
on graphs”. In: arXiv preprint arXiv:1408.1717.
The author would like to especially thank Pierre Van-
Komodakis, Nikos and Jean-Christophe Pesquet
dergheynstandNikolaosArvanitopoulosfortheircon-
(2014). “Playing with duality: An overview of
structive comments on the organization of the paper
recent primal-dual approaches for solving large-
and the experimental evaluation. He is also grateful
scale optimization problems”. In: arXiv preprint
to the authors of Dong et al. (2015a) for sharing their
arXiv:1406.5429.
code,toNathanaelPerraudinandNaumanShahidfor
Lake, Brenden and Joshua Tenenbaum (2010). “Dis-
discussions when developing the initial idea, and to
covering structure by learning sparse graph”. In: Pro-
AndreasLoukasforhiscommentsonthefinalversion.
ceedings of the 33rd Annual Cognitive Science Confer-
ence. Citeseer.
References Ng, Andrew Y, Michael I Jordan, Yair Weiss, et al.
(2002). “On spectral clustering: Analysis and an algo-
Banerjee, Onureena, Laurent El Ghaoui, and Alexan- rithm”. In: Advances in neural information processing
dre d’Aspremont (2008). “Model selection through systems 2, pp. 849–856.
sparsemaximumlikelihoodestimationformultivariate Perraudin,N.etal.(Feb.2014a).“UNLocBoXAmat-
gaussian or binary data”. In: The Journal of Machine lab convex optimization toolbox using proximal split-
Learning Research 9, pp. 485–516. ting methods”. In: ArXiv e-prints. arXiv:1402.0779.
Barab´asi, Albert-L´aszl´o and R´eka Albert (1999). Perraudin, Nathana¨el et al. (Aug. 2014b). “GSPBOX:
“Emergence of scaling in random networks”. In: Sci- A toolbox for signal processing on graphs”. In: ArXiv
ence 286.5439, pp. 509–512. e-prints. arXiv:1408.5781 [cs.IT].
Belkin, Mikhail and Partha Niyogi (2001). “Laplacian Shuman, David et al. (2013). “The emerging field
Eigenmaps and Spectral Techniques for Embedding of signal processing on graphs: Extending high-
and Clustering.” In: NIPS. Vol. 14, pp. 585–591. dimensionaldataanalysistonetworksandotherirreg-
Belkin, Mikhail, Partha Niyogi, and Vikas Sindhwani ular domains”. In: Signal Processing Magazine, IEEE
(2006). “Manifold regularization: A geometric frame- 30.3, pp. 83–98.
work for learning from labeled and unlabeled exam- Wang, Fei and Changshui Zhang (2008). “Label prop-
ples”. In: The Journal of Machine Learning Research agationthroughlinearneighborhoods”.In:Knowledge
7, pp. 2399–2434. and Data Engineering, IEEE Transactions on 20.1,
Cai, Deng et al. (2011). “Graph regularized nonneg- pp. 55–67.
ative matrix factorization for data representation”. Zhang, Fan and Edwin R Hancock (2008). “Graph
In: Pattern Analysis and Machine Intelligence, IEEE spectral image smoothing using the heat kernel”. In:
Transactions on 33.8, pp. 1548–1560. Pattern Recognition 41.11, pp. 3328–3342.
Combettes, Patrick L and Jean-Christophe Pesquet Zhang, Tong, Alexandrin Popescul, and Byron Dom
(2011). “Proximal splitting methods in signal process- (2006). “Linear prediction models with graph regular-
ing”. In: Fixed-point algorithms for inverse problems ization for web-page categorization”. In: Proceedings
in science and engineering. Springer, pp. 185–212. ofthe12thACMSIGKDDinternationalconferenceon
Daitch, Samuel I, Jonathan A Kelner, and Daniel A Knowledge discovery and data mining.ACM,pp.821–
Spielman(2009). “Fitting agraph tovector data”. In: 826.
Proceedings of the 26th Annual International Confer- Zhang, Yan-Ming et al. (2010). “Transductive Learn-
ence on Machine Learning. ACM, pp. 201–208. ing on Adaptive Graphs.” In: AAAI.
Dempster, Arthur P (1972). “Covariance selection”. Zheng, Miao et al. (2011). “Graph regularized sparse
In: Biometrics, pp. 157–175. coding for image representation”. In: Image Process-
Dong, Xiaowen et al. (2015a). “Laplacian Matrix ing, IEEE Transactions on 20.5, pp. 1327–1336.
Learning for Smooth Graph Signal Representation”. Zhu,Xiaojin,ZoubinGhahramani,JohnLafferty,etal.
In: Proceedings of IEEE ICASSP. (2003).“Semi-supervisedlearningusinggaussianfields
— (2015b). “Learning Laplacian Matrix in Smooth and harmonic functions”. In: ICML. Vol. 3, pp. 912–
Graph Signal Representations”. In: arXiv preprint 919.
arXiv:1406.7842v2.
Gilbert, Edgar N (1959). “Random graphs”. In: The
Annals of Mathematical Statistics, pp. 1141–1144.
Jiang, Bo et al. (2013). “Graph-Laplacian PCA:
Closed-form solution and robustness”. In: Computer
Vision and Pattern Recognition (CVPR), 2013 IEEE
Conference on. IEEE, pp. 3492–3498.
Manuscript under review by AISTATS 2016
A Derivations and proofs
0.4
Tikhonov: h(6) 1
0.35 /1+106
Gen.Model: h(6)/ 61 A.1 Detailed explanation of eq. (3)
0.3 HeatDi,usion: h(6) eqxp(!106)
/
0.25
6)(h 0.2 m m
0.15 (cid:107)W ◦Z(cid:107) =(cid:88)(cid:88)W (cid:107)x −x (cid:107)2
1 ij i j 2
0.1
i=1j=1
0.05
m m
(cid:88)(cid:88)
0 0.2 0.4 0.6 0.8 1 = (xi−xj)(cid:62)Wij(xi−xj)
6
i=1j=1
Figure 3: The filters of Table 3 for α=10. m m m m
(cid:88)(cid:88) (cid:88)(cid:88)
=2 x(cid:62)W x −2 x(cid:62)W x
i ij i i ij j
i=1j=1 i=1j=1
m m
=2(cid:88)x(cid:62)x (cid:88)W −2tr(cid:0)X(cid:62)WX(cid:1)
i i ij
i=1 j=1
=2tr(cid:0)X(cid:62)DX(cid:1)−2tr(cid:0)X(cid:62)WX(cid:1)
=2tr(cid:0)X(cid:62)LX(cid:1),
where D is the diagonal matrix with elements D =
ii
(cid:80)
W .
i ij
Tikhonov: h(6) 1
/ 1+106 A.2 Proof of proposition 2
1
0 Proof. We change variable W˜ =W/γ to obtain
-1 F(Z,α,β)=
=γargmin(cid:107)γW˜ ◦Z(cid:107) −α1(cid:62)log(γW˜1)+β(cid:107)γW˜(cid:107)2
Generative Model: h(6) / 61 W˜ 1,1 F
q 1 =γargminγ(cid:107)W˜ ◦Z(cid:107)1,1−α1(cid:62)log(W˜1)+βγ2(cid:107)W˜(cid:107)2F
W˜
0 α
=γargmin(cid:107)W˜ ◦Z(cid:107) − 1(cid:62)log(W˜1)+βγ(cid:107)W˜(cid:107)2
1,1 γ F
-1 W˜
(cid:18) (cid:19)
α
=γF Z, ,βγ ,
Heat Di,usion: h(6) exp(!106) γ
/ 1
where we used the fact that log(γW˜1) = log(W˜1)+
0 const.(W). The second equality is obtained from the
first one for γ =α.
-1
Figure 4: Different smooth signals on the Non Uni- A.3 Proof of proposition 3
form graph used for our artificial data experiments.
Proof. For equation (15)
All signals are obtained by smoothing the same initial
x ∼N(0,I)withthreedifferentfilters. Thisinstance
0 H(Z+γ,α,s)=
of the graph is disconnected with 2 components.
=argmin(cid:107)W ◦Z+γW(cid:107) +α(cid:107)W(cid:107)2 +α(cid:107)W1(cid:107)2
1,1 F
W∈Wm
s.t., (cid:107)W(cid:107) =s
1,1
=argmin(cid:107)W ◦Z(cid:107) +γ(cid:107)W(cid:107) +α(cid:107)W(cid:107)2 +α(cid:107)W1(cid:107)2
1,1 1,1 F
W∈Wm
s.t., (cid:107)W(cid:107) =s
1,1
=argmin(cid:107)W ◦Z(cid:107) +γs+α(cid:107)W(cid:107)2 +α(cid:107)W1(cid:107)2
1,1 F
W∈Wm
s.t., (cid:107)W(cid:107) =s
1,1
=H(Z,α,s),