Consistency-driven Alternating Optimization for Multi-Graph Matching

16 Pages·2014·1.84 MB·English

JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 1
Consistency-driven Alternating Optimization for
Multi-Graph Matching: a Unified Approach
Junchi Yan, Jun Wang Member, IEEE, Hongyuan Zha, Xiaokang Yang Senior Member, IEEE and Stephen Chu
Abstract—The problem of graph matching in general is NP- GRAPH matching (GM) plays a critical role in both
complete and many approximate pairwise matching techniques theoryandapplicationforcomputersciencecommunity.
have been proposed. For a general setting in real applications,
Briefly speaking, graph matching aims to find correspondence
it typically requires to find the consistent matchings across a
between two feature sets, with a wide spectrum of applica-
batch of graphs. Sequentially performing pairwise matching is
pronetoerrorpropagationalongthepairwisematchingsequence, tions that require feature matching in image processing and
and the sequences generated in different pairwise matching computer vision, as diverse as image registration [2], object
orderscanleadtocontradictorysolutions.Motivatedbydevising recognition [3], shape matching [4], object tracking [5], and
a robust and consistent multiple-graph matching model, we
action recognition [6], among others. Different from the point
propose a unified alternating optimization framework for multi-
based matching methods such as RANSAC [7] and Iterative
graph matching. In addition, we define and use two metrics
related to graph-wise and pairwise consistencies. The former Closet Point (ICP) [8], GM methods incorporate both the
is used to find an appropriate reference graph which induces unary node-to-node, as well as the second-order edge-to-edge
a set of basis variables and launches the iteration procedure. similarity as structural information. As discussed in [9], by
The latter defines the order in which the considered graphs
encoding the additional second-order edge similarity in the
in the iterations are manipulated. We show two embodiments
graph representation and matching process, GM methods can
under the proposed framework that can cope with the non-
factorizedandfactorizedaffinitymatrix,respectively.Ourmulti- usuallyleadtobetternodecorrespondencesolutions.Although
graph matching model has two major characters: i) the affinity extensive research has been done for decades [10], deriving
informationacrossmultiplegraphsareexploredineachiteration optimal GM is still a challenging problem in both theory
by fixing part of the matching variables via a consistency-
and practice since the GM problem can be formulated as a
driven mechanism; ii) the framework is flexible to incorporate
quadratic assignment problem (QAP) [11], being well-known
various existing pairwise graph matching solvers in an “out-
of-box” fashion, and also can proceed with the output of other NP-complete [12][13].
multi-graphmatchingmethods.Theexperimentalresultsonboth Graph matching is mostly considered under the two-graph
syntheticdataandrealimagesempiricallyshowthattheproposed scenario[9],[14],[15],[16],[17],[18]etc.However,inmany
framework performs competitively with the state-of-the-art.
realapplications,givenabatchofgraphsreferringtoidentical
or related structures, it is required to find the global node
mappings among all graphs, which is essentially related to
I. INTRODUCTION
graphclustering[19],[20],[21],[22],classification[23],[24],
[25]andindexing[26].Duetotheadvanceofmodernimaging
and scanning technologies, there is an increasing need for
Apreliminaryversionofthepresentedpaperappearedin[1].Thecurrent
multi-graph matching in various applications. In [27], objects
papermakesseveralfurtherextensionsandimprovements:i)incorporatethe
factorizedgraphmatchinginourformulation,andsolveitviaapath-following arescannedusinginfrared,optical,cartographicandSynthetic
algorithm; ii) propose two metrics for graph-wise and pairwise consistency, Aperture Radar (SAR) that are desired to be combined to
whichareusedfortheeffectivefindingofbasisvariablesandupdatingorder
improve representations. In addition, 3D shape analysis often
forthealternatingoptimizationprocedure.Inparticular,thepreliminarywork
doesnotaddressthesetwoproblemsasonlyafewgraphsaretestedin[1]. requires to model objects by multi-view assembly [28].
However,theybecomemorecriticalwhenarelativelylargenumberofgraphs Fromthemethodologyperspective,pairwisegraphmatching
are considered; iii) perform extensive evaluations that involves more graph
methods aim to find the correspondences that preserve the
samples,additionaldatasets,andmorestate-of-the-artsforcomparison.
J.Yan(correspondenceauthor)iswiththeDepartmentofElectronicEngi- structural invariance between two graphs. Though being suc-
neering,ShanghaiJiaoTongUniversity,Shanghai,200240China.Meanwhile, cessfulwhenthestructuralinvariancepropertyisroughlyhold,
heisalsosecondarilyaffiliatedwithIBMResearch–China,Shanghai,201203
pairwise matching methods tend to fail when the differences
China.e-mail:[email protected]
J.WangiswithInstituteofDataScienceandTechnology,AlibabaGroup, between the input graphs are significant or the graph noises
Seattle,WA,USA.e-mail:[email protected] are nontrivial. This is because the optimization solver may be
H.ZhaiswithSoftwareEngineeringInstitute,EastChinaNormalUniver-
trapped in local optima and/or far from the semantic ground
sity,Shanghai,200062China,andtheSchoolofComputationalScienceand
Engineering,CollegeofComputing,GeorgiaInstituteofTechnology,Atlanta, truth especially when large noises exist. However, a set of
Georgia,30332,USA.e-mail:[email protected] graphs open the possibility to explore the global information
X. Yang is with Shanghai Key Laboratory of Media Processing and
andmatchingconsistencyregularizationthathelpimprovethe
Transmission, Department of Electronic Engineering, Shanghai Jiao Tong
University,Shanghai,200240China.e-mail:[email protected] matching accuracy. For instance, directly producing node cor-
S. Chu is with IBM Research – China, Shanghai, 201203 China. e-mail: respondences between two significantly deformed graphs via
[email protected]
pairwise matching techniques is extremely error-prone; while
The work is supported by NSF IIS-1116886, NSF IIS-1049694, NSFC
61129001,NSFC61025005andthe111Program(B07022). ideally, an indirect matching with a sequence of interpolating
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 2
graphsyieldsmorerobustresults.Forarbitrarygraphs,simply extensive empirical tests.
applying the existing pairwise matching algorithms on each
pairofgraphisfarfromsatisfaction.Suchapairwisematching II. RELATEDWORK
based scheme encodes no global reasoning and regularization As a general problem for matching structural data, graph
due to the separated and independent matching process, re- matchinghasbeenextensivelystudiedfordecadesnotonlyin
sulting in information-loss and vulnerability to local noises. image processing and computer vision, but also in computer
Moreover, different orders of sequential pairwise matching scienceandmathematics[10],[34].Hereweviewtheproblem
would usually lead to inconsistent matching results, which from several key aspects that account for the main threads of
means XTiaXaj (cid:54)= XTibXbj for matching graph i and j, and the related work, mostly in the area of image processing and
the anchor graph is a and b respectively. Here X is the node- pattern recognition. We will discuss pairwise matching and
mappingassignmentmatrixwhoseformaldefinitionwouldbe multiple graph matching separately.
given later in the paper. This fact has been illustrated and
discussed in our conference paper [1].
A. Advances on pairwise graph matching
Compared with the vast existing work on the pairwise
Machine Learning Methodologies: Conventional graph
graph matching problem [9], [14], [15], [16], [17], [18],
matching methods first compute an affinity matrix and keep
[29], [30] etc., the topic of multi-graph matching is in fact
the affinity metrics unchanged during the entire matching
a relatively less-investigated one in the image processing
process. Recent work leverage various leaning algorithms for
and computer vision community [27], [31], [32]. Motivated
estimating the optimal affinity matrix [9], [17], [35], [36],
by the limitations of the existing approaches, we aim to
[37], and the methods can fall into either supervised [17],
design a robust method to obtain consistent node-to-node
unsupervised [9], or semi-supervised [35] learning paradigm.
correspondenceacrossacollectionoflooselyrelatedobjectsor
The problem of estimating the optimal setting of the affinity
weightedgraphs.Despitethefactthatnotheoreticalguarantee
metricisoutofthescopeofthispaper.Weassumetheaffinity
isprovidedforobtainingagloballyoptimalsolutions(notethe
matrix is prior known which is the same with most of related
problem is NP-complete), in this paper, we do take account
work [4], [14], [15], [16], [29], [38], [39] etc.
the global information of the entire set of graphs to generate
Higher-Order Affinity Modeling: Combing the unary and
approximated solutions that tend to satisfy the structural
second-order edge information has been heavily investigated
information and matching consistency across all the graphs.
since such a type of matching schemes plays a good tradeoff
To evaluate the effectiveness of the proposed multi-graph
between computational complexity and representation capa-
matching approach, we perform extensive empirical studies
bility [4], [14], [15], [38], [39]. Recently, the higher-order
on both synthetic datasets that are built using the popular
information has been encoded to achieve more robust match-
protocols [15], [29], as well as the well-known benchmark
ing paradigms. Several representative hyper-graph matching
datasets including the CMU motion sequence, the POSE-
methodshavebeenproposedforpatternrecognitionandimage
sequence, and the WILLOW-ObjectClass data.
processing applications [35], [40], [41], [42], [43] that encode
The major contributions of our work are:
higher-order information to enhance the matching distinctive-
First, we propose a unified framework, which has two
ness. However, a slight increase in the order would lead to a
alternating optimization algorithm variants to cope with non-
combinatorialexplosionofthestatespace.Asaresult,mostof
factorized and factorized affinity matrix in the objective func-
higher-order methods, such as the related work we mentioned
tion, respectively. Our method bears several merits: i) the
here,aretypicallyappliedonverysparsegraphswithnomore
pairwiseaffinitiesovermultiplegraphsarejointlyexplored,in
than third-order features. In addition, it is worth noting that
the hope of being more robust against local noises, which is
the term pairwise used in this paper refers to matching two
empirically observed in our extensive tests; ii) the framework
graphs,whileinthecontextofhyper-graphmatching,pairwise
is flexible to incorporate various existing pairwise matching
sometimes refers to the second-order edge affinity.
techniques in an “out-of-box” fashion, in two stages for both
Optimization Methods: Most approaches first formulate
before and during the iteration; iii) our method does not
an objective function, and then employ certain optimization
impose any consistency requirement on the initial matching
methods to derive optimal matching results [4], [14], [38],
inputs.Bycontrast,othermethods[32],[33]requiretheinitial
which can vary among a wide spectrum of optimization
solutions must be inconsistent, which allows them to improve
techniques and strategies. Several recent work first relax the
the overall accuracy by enforcing consistency. Therefore, our
objectivefunctiontoaconvex-concaveformulation[29],[44].
methodcanstart with,andfurtherimprove theresultsof[32],
Then the optimal solutions are obtained using the so-called
[33] also in an “out-of-box” fashion.
path following strategy and a modified version of the Frank-
Second, we define and use the metrics related to graph-
Wolfe algorithm [45]. Probabilistic matching paradigms are
wise/pairwise consistency to address two key challenges in
also developed, which have shown power in interpreting and
our framework: i) appropriate selection of the reference graph
addressing hyper-graph matching problems [15], [42], [43].
which induces a set of basis variables to proceed the iterative
optimization procedure effectively; ii) adaptive setting of the
B. Advances on multiple graph matching
iterative updating order, which improves the convergence
speed by a notable margin. Compared with the preliminary Themulti-graphmatchingproblemisalsoreferredasCom-
version [1], the consistency-driven framework outperforms in mon Labeling [31] or Multiple Isomorphism [46]. In [31], a
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 3
common labeling is defined as a bijective mapping between
all graph nodes in the considered graphs to a virtual node set.
The common labeling is constructed by a consistent multiple
isomorphism [46], where an isomorphism involves assigning
each node from one graph to one of the other graphs.
Comparedwiththepairwisematchingproblem,themultiple
graph matching has not been extensively studied, and only
a few address this problem using principled formulations
(a) twographsandtheassociatedaffinitymatrix (b) thestartreestructure
or techniques. A loosely related early work [47] aims at
finding a maximum spanning tree on a super graph whose
Fig.1:Illustrationsfor:a)thepairwisematchingaffinitymatrixinducedbytwosample
vertices represent graphs, and its edges represent matching graphs. Darker cells in the matrix denote smaller value within the normalized value
range[0,1];b)thestartreestructureforthereferencegraphtogetherwiththespecified
correspondencesbetweentwographs.Theedgesonthesuper-
“updated” graph and the resting “fixed” graphs in one iteration of the proposed two
graphareweightedbytheir“quality”suchasmatchingaffinity methods.Fivepairwisetermsareconsidered:Kru,Kf1u,Kf2u,Kf3u,Kf4u.
score as we will describe later in this paper. Very recently,
Sole-Ribalta and Serratosa [48] (and its journal version [31])
Specifically, his paper intensively uses X to denote the
ij
generalize the classical Graduated Assignment Graph Match- pairwise matching matrix between graph i, j, and X =
ing (GAGM) algorithm [14], [16] from two-graph to multi-
{X ,i,j =1,...,N}forthesetofpairwisematchingmatrix
ij
graphcase,whichgeneratesthecommonlabelingbymatching over a graph set denoted by G = {G ,...,G }. We call X
1 N
all graph nodes to a virtual node set. This method is assumed thematchingconfigurationw.r.t.G.Ifnotexplicitlystated,we
toinheritthesoundperformanceoftheoriginalmethod,while follow the convention that using N for the cardinality of G,
being more time-consuming as it repeatedly applies GAGM
n for the number of nodes in a graph, and m for edges.
across graph pairs iteratively. Other two recent work [33],
Now we introduce two definitions regarding with matching
[32]starttheirprocedureonaninconsistentpairwisematching
consistency induced by the initial pairwise matching results
configuration that covers all graph pairs. Given N graphs
fromanypairwisematchingsolver.Aswouldbeshownlaterin
with n nodes per each, the former method first builds an
thispaper,thesetwodefinitionsplayakeyroleintheproposed
N-dimensional probabilistic hyper-cube with size n in each
alternating optimization framework. In general, they are used
dimension, whereby each cell indicates the probability for the
to solve two common problems associated with an alternating
N-nodecorrespondencefromN graphsrespectively.Thenthe
optimization method, respectively: i) how to choose a set of
consistency and binarization are fulfilled in a post-processing
base pairwise matching variables to proceed the iteration; ii)
step.Thelatteremploysspectralapproximationtoeigenvector
how to set an optimal rotating order for alternating updating.
decomposition on the matching configuration matrix stacked
by all raw pairwise matching solutions (assignment matrix), Definition 1. Given N graphs G={Gk,k = 1,...,N}
and recover the consistent matching solutions. There are also and the pairwise matching configuration X = {Xij,i,j =
severalothermorerecentworkonmatchingabatchofgraphs 1,...,N}, the graph-wise consistency of graph Gk is defined
by self-boosting [49] and multi-view point registration [50]. asC (G ,X)=1−(cid:80)Ni=−11(cid:80)Nj=i+1(cid:107)Xij−XikXkj(cid:107)F/2 ∈[0,1]where
g k nN(N−1)/2
X is the pairwise permutation matrix over the graph set G.
ij
III. PRELIMINARIESFORGRAPHMATCHING Definition 2. Given a set of graphs G={G ,k = 1,...,N},
k
In this section, several notations used in this paper and two matching configuration X, and a reference graph G , for any
r
definitions related to the proposed framework will be intro- othergraphG ,thepairwiseconsistencybetweenG andG is
u u r
dpuaicrewdi,sfeogllroawphedmbaytcahrientgrousnpdeecrtiotwnoonpaerxaidsitginmgsf:onrmonu-lfaatciotonrsizfeodr defined as Cp(Gu,Gr,X)=1− (cid:80)Ni=1(cid:107)Xri−nNXruXui(cid:107)F/2 ∈[0,1]
and factorized representations of the affinity matrix. We leave the details of using the two definitions to Section
IV in the context of the proposed multi-graph matching
framework. In the rest of this section, we depict the two
A. Notations and definitions
existing pairwise graph matching formulations, which are the
Throughout the paper, R denotes the real number domain.
origins of the proposed multi-graph matching algorithms.
Bold capital letters denote for a matrix X, bold lower-case
letters for a column vector x, and hollow bold letters for a set
B. Non-factorized pairwise graph matching
X. All non-bold letters represent scalars. XT is the transpose
GiventwographsG (V ,E ,A )andG (V ,E ,A ),where
of X; diag(x) is a diagonal matrix whose diagonal elements 1 1 1 1 2 2 2 2
V denotes nodes, E, edges and A, attributes. There is
are x; vec(X) is the column-wise vectorized matrix X. tr(X) an affinity matrix K ∈ Rn1n2×n1n2 defined such that
is the trace of matrix X. In ∈Rn×n is an identity matrix and Kia;jb,(i,j = 1,...,n1),(a,b = 1,...,n2) measures the
the subscript n will be omitted when it can be inferred from edge pair affinity (vi,vj) vs. (va,vb) from two graphs. The
context. 1m×n,0m×n ∈ Rm×n are matrices of all elements diagonal term Kia;ia describes the unary affinity of a node
match (v ,v ). Rigorously, the affinity value K for the
being ones and zeros. 1 is the abbreviation for 1 . X◦Y i a ia;jb
n n×1 edge pair (v ,v ) vs. (v ,v ) is located at the ((a-1)n +i)-th
and X⊗Y denote the Hadamard and Kronecker products of i j a b 1
row and ((b-1)n +j)-th column of K. One example of K is
matrices. |X| is the cardinality of the set X, (cid:107)X(cid:107)F =tr(XTX) illustrated in Fig2.1(a). By introducing an assignment matrix
for the Frobenious norm and (cid:107)·(cid:107)p is the p-norm. X∈{0,1}n1×n2 for n1 =|V1|,n2 =|V2|, we set Xia =1 if
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 4
node vi matches node va (0 otherwise). The problem of GM matrix,andKq12 ∈Rm1×m2 fortheedgeaffinitymatrix,which
involves finding the optimal correspondence X, such that the is used to measure the similarity of each pair of nodes and
sumofthenodeandedgecompatibilitybetweentwographsis edges, respectively. As shown in [29], then the affinity matrix
maximized.Withoutlossofgenerality,weassumen ≥n for for pairwise matching can be factorized to:
1 2
different sizes of graphs. This leads to the following two-way
constrained quadratic assignment problem (QAP): K12 =(H2⊗H1)diag(vec(L12))(H2⊗H1)T
X∗ =argmaxvec(X)TKvec(X) where Hi =[Gi,Ini]∈{0,1}ni×(mi+ni), i=1,2
X (cid:20) Kq −Kq GT (cid:21)
s.t. X1n2 ≤1n1 1Tn1X=1Tn2 X∈{0,1}n1×n2 L12 = −1G21Kq12 G1K12q12G2T2 +Kp12
The constraints refer to the two-way one-to-one node map- BasedontheabovefactorizationofK ,theoriginalpairwise
12
ping:anodefromgraphG canmatchatmostonenodeinG graph matching objective function can be written as [29]:
1 2
andeverynodeinG2iscorrespondingtoonenodeinG1.There J12=vec(X12)T(H2⊗H1)diag(vec(L12))(H2⊗H1)Tvec(X12)
is no (one/many)-to-many matchings between two graphs. (cid:16) (cid:17)
The above formulation is general as it allows two graphs =tr L12(HT1X12H2)◦(HT1X12H2
having unequal number of nodes (n1 (cid:54)= n2). Similar to By factorizing L to L = U VT = (cid:80)c ui viT,
the previous work [12], [29], we follow the protocol that 12 12 12 12 i=1 12 12
cpoenrmveurttastiXonfmroamtrixan(Xa1snsi2g=nm1enn1t)mbyataridxdi(nXg1dnu2m≤my1nn1o)detsotoa [wvh112e,rev212U,1·2·· ,=vc12[]u∈112R,u(n2122+,·m·2·),×ucc1,2J]12∈canRb(en1d+emriv1)e×dca,sV[29]=:
oanndegaruagpmhe(ni.te.thaeddaifnfignsitlyacmkvatarriixabblyeszteorothse)ainssicgansmeenn1tm(cid:54)=atnri2x. J12 =(cid:88)tr(cid:16)Ai12X12Bi12XT12(cid:17) (2)
This is a standard technique from linear programming [51] i
and is adopted by [12], [29], [52] etc., being able to handle where Ai =H diag(ui )HT, Bi =H diag(vi )HT
12 1 12 1 12 2 12 2
superfluous nodes in a statistically robust manner. By taking
this step, the graphs are of equal sizes. This preprocessing Notetheabovefactorizationformulationdoesnotrequiretwo
opens up the applicability of existing multi-graph methods graphs being of equal sizes. In the rest of paper, similar to
[31],[32],[46]astheyareallbasedontheassumptionthatall
the non-factorization case, we will add dummy nodes in case
graphsareofequalsizes.Therefore,thefollowingformulation
of unequal sizes of graph, to make the two-way constraint
is used throughout the paper which assumes n =n =n.
1 2
become a bijection. This step will transform X to a strict
12
X∗12 =argmaxvec(X12)TK12vec(X12) (1) permutation matrix, which leads to the following convex-
X12
s.t. X 1 =1 1T X =1T X ∈{0,1}n1×n2 concave relaxation formulation as discussed in [29].
12 n2 n1 n1 12 n2 12
Here X is a permutation matrix which is constructed by
12 D. Convex-concaverelaxationsforfactorizedgraphmatching
augmenting the assignment matrix with slack columns in case
By ensuring X a permutation matrix such that XT X =
n >n ,suchthatthereisalwaysabijectionbetweengraphs. 12 12 12
1 2 I , the objective can be relaxed to the convex formulation,
The above QAP formulation is used by most existing bnynaddingaconstanttermC asusedin[29].Thuswehave:
12
GM methods [4], [14], [15], [16], [38], which directly deals
1
with the large pairwise affinity matrix K. The affinity matrix J1v2ex(X12)=J12(X12)− 2C12(X12) (3)
polradyesr aancdenstercaolnrdo-loerdinerGrMelatbieocnasusbeetiwteeenncogdreasphasl.l Tthheefitrwsto- where C12=(cid:88)(cid:16)tr(Ai12Ai12X12XT12)+tr(Bi12Bi12XT12X12)(cid:17)
i
common properties of K, full-rankness and indefiniteness,
On the other hand, as shown in [29], the concave relaxation
pose key challenges to conventional GM methods such as Jcav can be written as:
12
spectralmethods[4],[38]andgradient-basedapproaches[14],
(cid:16) (cid:17)
[15], [16]. In addition, its relatively large size also impedes J1c2av =tr Kq1T2(GT1X12G2◦GT1X12G2) (4)
the applicability when large-size graphs are considered. Thus (cid:16) (cid:17) (cid:16) (cid:17)
−tr (G Kq GT)TX +tr KpTX
the factorized graph matching formulation is studied and 1 12 2 12 12 12
introduced by [29] in parallel which avoids involving the
whole K for optimization.
IV. MULTIPLEGRAPHMATCHINGBY
C. Factorized pairwise graph matching CONSISTENCY-DRIVENALTERNATINGOPTIMIZATION
The authors in [29] show that the affinity matrix can First of all, as the same with the preprocessing used in
be factorized as a Kronecker product of smaller matrices,
pairwise graph matching that makes graphs of equal size, we
which decouples the graph structure from the affinity. Specif-
add dummy nodes to make all input graphs of the same size.
ically a graph G can be denoted by {P,Q,G}, where P =
[p ,...,p ] ∈ Rdp×n, Q = [q ,...,q ] ∈ Rdq×m are the Note this protocol is also used by [31], [32], [46] for multi-
1 n 1 m
feature matrices computed for n nodes and m edges. Here d graph matching. Based on this consensus, in what follows we
p
and d are the number of dimensions for unary features and will present our algorithm with two variants, by factorizing
q
edge features. The topology of the graph is specified by the the affinity matrix and not, respectively. The non-factorized
node-edge incidence matrix G ∈ Rn×m such that the non-
formulation is a more dominant representation used in most
zero elements in each column of G indicate the starting and
existingwork[4],[39],[15],[16],whiletheotheravoidsusing
ending nodes in the corresponding edge. Refer to [29] for
moredetails.Morespecifically,giventwographs{P ,Q ,G } the whole affinity matrix in optimization [29], [30]. We will
1 1 1
and {P2,Q2,G2}, let Kp12 ∈Rn1×n2 denote the node affinity provide a unified approach to cope with these two cases.
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 5
A. Non-factorized multiple graph matching Algorithm 1 Consistency-driven Non-factorized Alternating
Given N graphs and the pairwise affinity matrix K for Optimization for Multi-Graph Matching
ij
eachpairofgraphsGi,Gj,withoutlossofgenerality,themulti- Input:
graph matching objective function can be written as follows: 1: N graphswithnnodesofeachgraph;
N 2: PairwiseaffinitymatrixKij(i,j=1,...,N);
X∗ =argmax (cid:88) vec(X )TK vec(X ) (5) 3: Maximumiterationcount:Tmax;
ij ij ij Output:
X
i,j=1,i(cid:54)=j 4: ConsistentassignmentmatrixXij (i,j=1,...,N);
s.t. X 1 =1 1TX =1T X =XT ∈{0,1}n×n Procedure:
ij n∀i,j =n1,.n..,iNj ; |nX|=iNj (N −ji1) 5: pOebrtfaoirnmitnhgepraaiwrwmiseatcghrainpghmcoantcfihgiunrgatoiovnerXN=grap{hXs;ij} by exhaustively
(awughmereentXed a=ssig{nXmije,nit,)jma=trix1,s.e.t.s,uNch}thiast Xthiej =peXrmTjiu.taNtiootne 768::: SSIneeitttiuraelpifdzeearetitnnhcgeeblgiesrstatpOshouGldutrtioi=nnsamsXcae∗kxnrGdk=inCgXgok(rrGd,ekkr,wX=.)r,.1t{.,k.C.=p.(,GN1u,,.,k.G.(cid:54)=r,,NXr;});
that by the definition of affinity matrix as introduced earlier 9: fort=1:Tmax do
in the paper, in general Kij (cid:54)= Kji, Kij (cid:54)= KTji. Though 1101:: forFuixinNO-2uXdttufd,oforf =1,...,N,f (cid:54)=r,u,updateXturbysolving
they encode redundant information, we include all Kij in the thetwo-graphmatchingproblem(7);
objective function, which would be used in the optimization 12: Compute the updated objective score J(Xur) using the updated
Xt fortheobjective(7);
procedure as we will show in the following. ur
One obvious observation on the above function is that 13: Computetheso-farbestobjectivescoreJ∗(Xur)usingtheso-far-
bestX∗ fortheobjective(7);
the pairwise matching variables X are redundant and can ur
be determined by a compact basisijvariable set {X ,k = 14: ifJur >Ju∗r,X∗ur =Xtur;
rk 15: endfor
1,...,N,k (cid:54)= r} (a star tree rooted at the reference graph) 16: endfor
such that any X can be computed by XT X . Consequently,
ij ki kj
we design our non-factorized multi-graph matching algorithm
as follows. First, in order to induce a basis matching set to Algorithm 2 Consistency-driven Factorized Alternating Opti-
proceedtheoptimizationprocedure,areferencegraphGr isset mization for Multi-Graph Matching
by a certain means. Then in each iteration, one can choose a
certaingraphG inrotationbyacertainorderO ,andupdate Input:
u udt
its mapping with G i.e. X , by fixing the other N −2 X 1: N graphswithnnodesofeachgraph;
r ur rf
regardinggraphsG (f =1,...,N,f (cid:54)=r,u)againstG .This 2: Graph-wisenode-edgeincidencematrixGi,nodeandedgeaffinitymatrix
f r Kp,Kq (i,j=1,...,N);
idea is illustrated in Fig.1(b). Now, we reach the following ij ij
objective function by dropping the constant terms1 for X : 3: Maximumiterationcount:Tmax,path-followingstepsizeδ
rf Output:
N 4: ConsistentassignmentmatrixXij(i,j=1,...,N);
vec(Xur)TKurvec(Xur)+f=1(cid:88),f(cid:54)=r,uvec(Xuf)TKufvec(Xuf) P5r:ocOebdtuarine:the raw matching configuration X = {Xij} by exhaustively
performingtwo-graphmatchingoverN graphs;
tNhiosteleiandsthteoptehremfuotlaltoiwoninmgaetrqiuxaftoiornmiwnethheavveecXtourfiz=edXfuorrmXr:f, 67:: SSeettruepfdearetinncgelgisrtapOhuGdtri=namscaexnGdkinCggo(rGdekr,wX.)r,.t{.kC=p(G1u,.,.G.r,,NX});
vec(X )=(X ⊗I)vec(X ) (6) 8: InitializethebestsolutionsX∗rk=Xrk,k=1,...,N,k(cid:54)=r;
uf fr ur 9: fort=1:Tmax do
FromnowonwewilluseF ∈Rn2×1 todenoteX ⊗I,by 10: foruinOudt do
fr fr 11: FixN-2Xt ,forf =1,...,N,f (cid:54)=r,u,updateXt byusing
replacing vec(X ) via vec(X ) = F vec(X ) we rewrite fu ru
uf uf fr ur thefollowingpath-followingprocedure:
the objective function more compactly as follows: 12: forα=0:δ:1do
13: update Xt via applying MFW on the convex-concave relax-
N ru
J(X )=vec(X )T(K + (cid:88) FT K F )vec(X ) (7) ation:J =(1−α)Jvex+αJcav byformula(10).
ur ur ur fr uf fr ur 14: endfor
f=1,f(cid:54)=r,u 15: Compute the updated objective score J(Xru) using the updated
It becomes clear that in each iteration the sub-problem (7) Xtru fortheobjective(9);
is a standard pairwise graph matching problem, which can be 16: Computetheso-farbestobjectivescoreJ∗(Xru)usingtheso-far-
bestX∗ fortheobjective(9);
ru
solved by various pairwise graph matching techniques based 17: ifJru>Jr∗u,X∗ru=Xtru;
on the QAP formulation in an “out-of-the-box” manner2 18: endfor
19: endfor
So far, we have proposed our alternating optimization
framework that directly deals with the affinity matrix. Note
thatintheabovediscussion,weassumethereferencegraphG
r
andtheupdatingordero arebothpre-given.Wewouldstill B. Factorized multiple graph matching
udt
holdthisassumptionforpresentingthefactorizedvariantinthe
Recent work [29] exploits the underlying structure of affin-
next paragraph. In the end of this section, we would address
ity matrix K, which is not necessarily negative definite, in the
these two common problems using a principled paradigm.
hope of designing better optimization scheme for addressing
As we will show later, it is closely related to the graph-
thenon-convexissue.Thepresentedalgorithmintheprevious
wise/pairwise consistency metrics as defined in Section III.
section cannot be directly applied to the original factorized
graphmatchingformulationin[29]asthere-weightedaffinity
1Forefficiency,weomitthetermsregardingKfuastheyhavebeenencoded matrix in (7) lacks a convenient interpretation as the original
byKuf.Thispolicyisusedthoughtthepaperfortheproposedalgorithms. affinitymatrixthatcanbefactorized(decoupled)intopairwise
2Several matching methods are not or not explicitely based on QAP: e.g.
similarity and self-structure. In this section, we present an
thegrapheditdistancebasedmethods[53]andtherecentadvances[54],[55]
canenjoybettererror-toleranceagainstmissingnodesandedges. alternating optimization approach to solve the multiple graph
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 6
matching problem based on the factorized formulation. As where the gradient can be obtained as follows:
pointed out by [29], the main advantage of factorized model
N
is avoiding computing and storing the large-sized affinity ∇ Jα =(1−α)∇ (Jvex+ (cid:88) Jvex) (12)
Xru Xru ru fu
matrix by decoupling it into smaller edge-wise and node-wise
f=1,f(cid:54)=r,u
matricesaswehaveshowninthepairwisegraphmatchingfor- N
mulation. Another notable difference from our non-factorized +α∇Xru(Jrcuav+ (cid:88) Jfcauv)
multi-graph matching algorithm, is that it derives a specific f=1,f(cid:54)=r,u
convex-concave path-following algorithm as we will present Having obtained the optimal Y, the optimal step size λ can
inthissection,whiletheformercanre-useanynon-factorized be found at the optimal point of the following parabola:
graph matching solvers [15], [16], [38], [39].
Jα =(1−α)Jvex(X +λY)+αJcav(X +λY) (13)
In general, the factorized multiple graph matching problem ru ru
can be expressed by maximizing the following objective =aλ2+bλ+const
function, where J is in the form of formula (2):
ij
Next, we would focus on solving the two common issues
N
X∗ =argmax (cid:88) J (X ) (8) for both non-factorized and factorized variants: i) finding the
X ij ij optimal reference graph G that induces the basis variable
i,j=1,i(cid:54)=j r
set {X ,k = 1,...,N,k (cid:54)= r} to initialize the iterative
s.t. X 1 =1 1TX =1T X =XT ∈{0,1}n×n rk
ij n n n ij n ij ji optimization procedure; ii) setting the optimal updating order
∀i,j =1,...,N; |X|=N(N −1)
O to speed up convergence and to suppress possible error
udt
For the factorized graph matching formulation, we will propagation over alternating optimization. This is because
directly use the matrix form X instead of the vectorized
although in each iteration the objective (7) or (9) contains
form x = vec(X) as used in the non-factorized case. As the
affinities over multiple graphs to explore jointly, nevertheless,
samewiththenon-factorizedmulti-graphmatchingproblem,a
reference graph Gr is first selected, followed by an alternating this new objective is coupled with the estimated Xrf from the
optimization procedure with a certain sequential updating or- previous iteration thus cannot fully avoid error propagation.
der.Similartothenon-factorizedcase,wereachthefollowing
objective function without using the affinity matrix:
C. Consistency-driven alternating optimization
N
(cid:88)
J(X )=J (X )+ J (X ) (9) Firstweaddresstheproblemrelatedtofindingthereference
ru ru ru fu fu
f=1,f(cid:54)=r,u graph. Our iterative alternating optimization involves a set of
Note J is a function of X as X = X X for fixed basis mappings between a reference graph and others in an
ru fu fr ru
X . Its convex-concave relaxation form, which is weighted expectation maximization manner. It is critical to generate
fr
by α∈[0,1], can be written as follows: reasonably accurate initial solutions to avoid trapping into
an unsatisfactory local optimum. Specifically, one is given a
Jα(X )=(1−α)(cid:16)Jvex(X )+ (cid:88)N Jvex(X )(cid:17) (10) raw matching configuration Xraw obtained from independent
ru ru ru fu fu
f=1,f(cid:54)=r,u pairwise matching without knowing the ground truth. Under
+α(cid:16)Jrcuav(Xru)+ (cid:88)N Jfcauv(Xfu)(cid:17) tGhrisancodnidtsitaiossno,coinaeteidsbsaeseiksisnegtathna“tasppparnoaprnieawte”mraetfcehreinngcecognrafipgh-
f=1,f(cid:54)=r,u uration Xspan by Xspan = XrawTXraw. One intuitive metric
ij ri rj
where Jvex is written by formula (3), and Jcav by formula to approximate the accuracy is the graph-wise affinity score
(4). As α increases from 0 to 1, this path-following strategy by(cid:80)N vec(X )TK vec(X )givenG .However,due
k=1,k(cid:54)=r kr kr kr r
enjoys the global optimal solution when α = 0 due to the to outliers and local deformation, as well as the difficulty in
convex form, and finally leads to an integer solution when setting up the perfect affinity matrix in a parametric manner3,
α = 1, for the concave form. Thus it is insensitive to the the ground truth may not correspond to the highest score
initial point, and usually no post-binarization is needed [56], modeled by the parametric objective function. Thus graph-
[57] which otherwise may cause additional performance loss. wise score is not a robust indicator to accuracy as shown in
For each α ∈ [0,1], the sub-problem can be solved using
Fig.2, which will be discussed in details in our experiments.
the Frank Wolfe’s algorithm (FW) [45], which is widely
Our embodiment is setting G = max C (G ,X) by
used constrained nonlinear programming. In implementation, r Gk g k
we use the modified Frank Wolfe’s algorithm (MFW) [58] Definition (1). The rationale is that the erroneous correspon-
to speed up its convergence which probably finds a better dences are at random yet correct correspondences concur
searching direction Y by a convex combination of previously thus consistent across pairwise matchings. This observation
obtained solutions. In the experiments, we used the directions
motivates us choosing G as the one maximizing the graph-
computed in 2 previous steps. The FW (MFW) algorithm first r
wise consistency – Definition (1), in the hope that it would
computes the optimal direction Y and then determines the
optimal step size λ ∈ [0,1], with respect to the objective inducethemostaccuratebasissetforinitialization,orequally:
function(10).TheoptimalYwithrespecttoX initerationk min (cid:80)N (cid:107)Xspan − XrawTXraw(cid:107) , where (cid:107)(cid:107) is the p-
ru r i,j=1 ij ri rj p p
canbecalculatedbysolvingthefollowinglinearprogramming norm for a matrix. Since the resultant matching matrices are
problem by the Hungarian method [59]:
(cid:16) (cid:17) 3Currentlytheaffinityfunctionismostlymodeledbyparametricfunctions,
max tr ∇ Jα(Xk )T(Y−Xk ) (11)
Y Xru ru ru thefixedparametersbywhatevermanualsetting[16],orlearnedfromtraining
samples[17],[36]etc.arestillunabletoperfectfitthescorewithaccuracy,
s.t. Y1 =1 ,YT1 =1 ,Y≥0
n n n n n×n leadingtotheexistenceofdiscrepancybetweenscoreandaccuracy.
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 7
TABLEII:ParametersettingsforsolutionpathillustrationasshowninFig.2.
noisetype parametersettings results
deform N=20,nin=12,nout=0,ε=.15,ρ=1,σ2=.05 Fig.2(a)
outlier N=20,nin=8,nout=4,ε=.05,ρ=1,σ2=.05 Fig.2(b)
density N=20,nin=12,nout=0,ε=.05,ρ=.5,σ2=.05 Fig.2(c)
D. Implementation details and convergence discussion
(a) Deform (b) Outlier (c) Density Another consideration is due to the sub-optimality nature
of existing pairwise matching techniques for solving the NP-
hard sub-problem of Eq.(7) or (9) per iteration. Therefore the
newly obtained solution cannot guarantee to achieve higher
Fig. 2: Comparison of four strategies for setting reference graph, updating order and
affinityscorethanthepreviousone,althoughoftenempirically
enforcingscore-non-descendingpathselectionconstraint,refertoTableIandTableII
for structured descriptions. Accuracy over iteration path by four solvers in different observed. To avoid the score degenerating case, we suggest
colors:RRWM(blue),GAGM(green),FGM(black),IPFP(red).Thecorrespondingfour
enforce a score-non-descending path selection strategy by
horizontalsolidlinesineachplotindicatetheperformanceofabaselinebyrandomly
choosingagraphasthereferenceonetoderiveconsistentmatchingconfiguration.The comparing the objective score calculated by the new solution
curvesaregeneratedbyaveraging50randomtestsoverTmax =2iterations,i.e.the Xt and the one from the currently maintained solution X∗ .
numberofinneriterationroundsis(N−1)∗Tmax=38.Zoominforbetterdisplay. ur ur
The new solution would be dropped if it decreases the score.
As such, we obtain a score-non-descending solution path.
TABLE I: Settings of reference graph, alternating updating order, and path selection
strategyforthesolutionpathsasillustratedinFig.2. When this strategy is not imposed, X would always be
ur
referencegraph updatingorder pathselection curve overwritten. We would compare both the “path selection” and
graphwiseconsistency pairwiseconsistency enforced dashed “non-path selection” strategies in our experiments.
graphwiseconsistency randomorder enforced solid
graphwiseconsistency pairwiseconsistency non-enforced dash-dot Now, we summarize the above discussion and analysis into
graphwisescore pairwisescore enforced dotted twovariantsofouralternatingoptimizationalgorithms,which
are depicted in Alg.1 and Alg.2 for the non-factorized and
factorized model respectively. The score-non-descending path
all binary ones, thus the setting of the reference graph is selection step is plugged in the two algorithm charts.
insensitive to p. Without loss of generality, we set p=2. Fig.2depictsacomparisonunderdifferentsettingsbyusing
Given the basis variable set for rotating updating, the next the proposed framework, regarding with the way of setting i)
problem is finding an “optimal” updating order. The pairwise the reference graph, ii) the updating order and iii) the path
consistency C (G ,G ,X) – Definition (2) is used to define selection strategy. Four pairwise matching solvers are used to
p u r
the order in which the sample graphs in the iterations are study if the solution path is impacted by the specific pairwise
considered, so that the least accurate ones are updated first, matching technique or if a more general pattern exists. We
thus avoiding error propagation and speeding up the conver- will give a detailed analysis to Fig.2 later in our experiments.
gence. Motivated by the similar observation in the case of In the last, we provide a general discussion about the
finding the reference graph, the pairwise consistency measure convergencebehavioroftheproposedmethods.Bothproposed
is considered to approximate the pairwise matching accuracy. algorithmsinvolveupdatingasetofbasisassignmentsolutions
As a result, the updating variables are ranked in ascending of length N: {Xt ,k = 1,...,N,k (cid:54)= r} w.r.t. the reference
rk
order with respect to C (G ,G ,X). Similarly, C (G ,G ,X) graphG initerationt.NowweconcatenateallindividualXt
p u r p u r r rk
is also insensitive to which matrix norm (cid:107)(cid:107) (we set p = 2 into a single matrix Xt = [Xt ,Xt ,...,Xt ]. Since each
p r r1 r2 rN
here) is used as only binary matching matrices are involved. Xt is a discrete binary permutation matrix thus Xt is also
rk r
Note the authors in [31], [33] also define an inconsistency exhaustively enumerable. As a result, as iteration continues,
index for the overall pairwise matching configuration X, or, there exist t and s (t<s), such that the t-th iteration solution
multipleisomorphismastermedintheirpapers.However,that Xrt wouldequaltoapreviousvalueXrs (Otherwiseitimplies
index is not related to the fine-grained graph-wise/pairwise the solution space is innumerable thus contradicts with the
consistency, and unable to define the reference graph and fact). Furthermore, because both algorithms are deterministic,
updating order. Perhaps more importantly, the two proposed thustheiterativesolutionpathwouldrestartthesolutionpath:
definitions in general provide the sketch for deriving other Xtrk →Xsrk =Xtrk →Xsrk →··· inaloopingmanner.There-
consistency metrics, which can be tailored in particular appli- fore, both algorithms would converge to a looping sequences
cations. For instance, the graph-wise consistency can be cal- when t<s or to a fixed point when t=s.
culated by weighting each term (cid:107)X −X X (cid:107) in Definition
ij ir rj
(1)byaparameterw ifthepriorestimationonthequalityof
ij V. EXPERIMENTSANDDISCUSSION
theinitialpairwisematchingsolutionisgiven.Insummary,the
presented two measures themselves already bear some gener- The experiments are performed on both synthetic and real-
alities since the binary matrices are insensitive to the choice image datasets. The synthetic test is controlled by quanti-
of p-norm and they can allow weighted variants as discussed tatively varying the disturbance of deformation, outlier and
above. Moreover, our general consistency-driven mechanism edge density. The real-image datasets are tested with varying
can also benefit from other newly designed metrics. viewing angles, scales, shapes, and spurious outliers.
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 8
The matching accuracy over all graphs, is calculated by the protocol of [29], [36] that set the final affinity matrix
averaging all pairwise matching accuracy (cid:80)iN=−11(cid:80)Nj=i+1Accij. re-weighted by the edge length affinity and angle affinity:
Each Accij computes the matches between thNe(Nc−o1r)r/e2spon- Kij,ab = βKilje,nab +(1−β)Kiajn,agb, where β ∈ [0,1] is the
dence matrix Xalg given by the ground truth Xtru: Acc = weighting parameter. This is because this dataset is more
ij ij ij
tr(XalgXtru) geometrically ambiguous than the sequence data if only edge
ij ij . Note we only calculate the accuracy for com-
tr(1nj×niXtijru) length is used. The angle for each edge is computed by the
mon inliers and ignore the matching results over outliers. absoluteanglebetweentheedgeandthehorizontallineasused
Theabovetestingmethodologiesfollowastandardprotocol in [29]. The edge affinity and angle affinity are calculated in
widelyadoptedbymanyrelatedworkssuchas[15],[16],[29]. the same way in the CMU-POSE test. For node-wise affinity,
we use the 128-dim SIFT feature [61] s ∈ R128 associated
with each single landmark by K =exp(−(cid:107)s −s (cid:107) ).
A. Dataset description and affinity setting ia;ia i a 2
Synthetic data The synthetic test follows the widely used
protocol of [1], [14], [15], [16], [29], [38]. For each trial, B. Comparing methods
a reference graph with nin nodes is created by assigning First we give a short description for the used pairwise
a random attribute to its edge, which is uniformly sampled matching solvers, which serve as building-blocks for our and
from the interval [0,1]. Then the “perturbed” graphs are other multi-graph matching algorithms [32], [33]. Other state-
created by adding a Gaussian deformation disturbance to the of-the-art multi-graph methods [31], [32], [33] are briefly
edge attribute qirj, which is sampled from N(0,ε) i.e. qipj described in the sequel. The tests run on a laptop with 2.9G
= qirj+N(0,ε) where the superscript ‘p’ and ‘r’ denotes for Intel Core I7 and 8G memory with a single thread. The code
“perturb”and“reference”respectively.Each“perturbed”graph of the comparing methods are all from original authors, apart
is further added by nout outliers, which can also be helpful from [14], which is implemented by the author of [29].
to make the graphs of equal sizes when the input graphs are We select the following widely used pairwise graph match-
different sizes. Its edge density is controlled by the density ing methods as the pairwise solvers in our methods. In
parameter ρ ∈ [0,1] via random sampling. The edge affinity particular, the parameter settings of these methods as reported
is computed by Kij,ab = exp(−(qij−σ2qab)2) where σ2 is the belowmostlyfollowtheoriginalsettingsfromtheauthors.We
edge similarity sensitivity parameter. No single-node feature further slightly tune the parameters to balance efficiency and
is used and the unary affinity Kii,aa is set to zero. efficacyasmatchingmultiplegraphsismoretime-consuming.
CMU-POSE Sequence This data contains three se- Graduated assignment (GAGM) GAGM [14] performs
quences. The first two are from the CMU house (30 gradient ascent on a relaxed QAP objective driven by the
marks, 101 frames), hotel (30 marks, 111 frames) sequence deterministic annealing procedure. The relaxation level is
(http://vasc.ri.cmu.edu//idb/html/motion/) which are widely controlledbyacontinuationparameterβ,whichisupdatedby
used in [1], [15], [16], [17], [29], [36]. The third sequence αβ →β ≤β . We set α=1.1, β = 0.5 and β =50.
t t+1 max 0 max
is sampled from the sedan (VolvoC70) sequence (19 marks, In the CMU-POSE data test, we also replace GAGM with
225 frames) which is viewed from various angles covering a its variant, i.e. Soft Constrained Graduated Assignment [16]
range of 70 degrees from the POSE dataset [60]. For each which is a recent improvement based on GAGM.
test, an image sequence is sampled which is spaced evenly Integer projected fixed point method (IPFP) IPFP [39]
by three frames. We utilize this dataset for the outlier test. approximates the original objective function via transforming
We select nin=10 landmarks out of all nant annotated points it into a series of linear assignments over iterations, and in
(nant=30 for the two CMU sequences, nant=19 for Pose eachiterationfindstheoptimalsolutioninthediscretedomain
sedan sequence), and randomly chose nout=4 nodes from viaHungarianmethod.Thenthemethodfindsthere-weighted
the rest nant−nin nodes as outliers. The graphs are fully solutioninthecontinuousdomainalongthedirectionfromthe
connected to encode all information to suppress outliers. current solution to the optimal discrete solution. We set the
The affinity matrix is constructed by the edge length maximum round of iterations as 10 in all experiments, which
similarity Klen = exp(−(dij−dab)2) where d , d are is in accordance with the empirical observation by [39].
ij,ab σ2 ij ab
the Euclidean distance between two points that are further Re-weighted random walk matching (RRWM) RRWM
normalized to [0,1] by dividing the largest edge length. The [15] introduces a random walk view on the problem and to
unary affinity is also set to zero in line with [15], [16], [29]. some extent can be regarded as a re-weighted version for
WILLOW-ObjectClassTheobjectclassdatasetreleasedin GAGM [14] and Spectral Matching [4] methods. We fix its
[36]isconstructedbyimagesfromCaltech-256andPASCAL parameters α=0.2 and β=30 in all experiments.
VOC2007. Each object category contains different number of Fast Bipartite graph matching (FBP) FBP [54] is a re-
images:109Face,50Duck,66Winebottle,40Motorbike,and cently proposed graph edit distance based method. It approxi-
40 Car images. For each image, 10 feature points were man- matesthesecond-ordergraphmatchingproblemtoafirst-order
ually labeled on the target object. The edge sampling for the assignment problem by considering the local edge structure
affinity matrix follows the same way as [29] by constructing via fast bipartite matching. It requires to build the node-to-
the sparse delaunay triangulation among the marks, since the node cost matrix C ∈ Rn×n which encodes both node-wise
triangulation can efficiently encode the object structure when and local structure information around each node. However,
no outlier exists. For defining the edge affinity, we follow in each iteration of our framework, only the (re-weighted)
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 9
affinitymatrixKisgivenwhiletheattributematrix(thispaper configuration4 (when N≤n) to obtain consistent matchings.
is dealing with attributed graphs) for each graph is unknown Permutation synchronization (SYNC) SYNC [32] em-
and cannot be recovered from K. Thus the attribute/adjacency ploys spectral analysis and approximation to eigenvector de-
matrix based steps for computing C as used in [54], [55] are composition on the matching configuration matrix comprised
inapplicable in our framework. To compute the element C ofallinitialpairwisematchingsolutions,andrecoverthecon-
ia
when only K is given, we employ the Hungarian method to sistent matching solutions. Note that the SYNC method only
compute the assignment cost for the rest of nodes, by fixing can work when the number of nodes is less than the number
the mapping i → a. Specifically, the input to the Hungarian of graphs, thus similar to the sub-step in the MHCL method,
methodistherow/columnofK(byinvertingthevaluesoasto Maximum Spanning Tree (MST) is used when n≤N.
change it from affinity to cost) that can be re-shapen into the
cost matrix w.r.t. i → a i.e. vec2mat(max(K)−K ) – see C. Incorporating both pairwise and multi-matching solvers
ia,:
Fig.1(a) for illustration. Note this adaption may weaken the Havingdescribedtheabovemulti-graphmatchingmethods,
performanceofFBPwhichisoriginallydesignedforeffective we provide a further comparison for the advantage of our
exploring the individual graph attribute/adjancey information. framework,whichcoverstwoaspects:i)betterreusingvarious
Factorized graph matching (FGM) FGM [29] factorizes pairwise matching techniques in an “out-of-the-box” fashion;
the affinity matrix which allows for relaxing the objective to ii) proceeding with the output of other multi-graph matching
a convex function and a concave one respectively. We set the methods to improve the solution.
increment step δ = 0.2 for the path-following algorithm, and For the first aspect, our framework not only leverages the
the iteration count in the Modified Frank-Wolfe method as 5. pairwise matching solver to compute the putative matchings
Ontheotherhand,thefollowingpeermulti-graphmatching in the first stage for initialization and consistency estimation,
methods are evaluated and compared with our methods: but also reuses the solver during the iterative alternating
Graduated assignment based common labeling (GACL) optimization. In this second stage, the affinity matrix is re-
GACL [48], [31] is a more recent work that extends the newed, containing the global information beyond two graphs.
graduated assignment method to solve the multi-graph match- In contrast, the peer methods MHCL/SYNC use the pairwise
ing problem, by matching all graph nodes to a virtual node matching technique in one-shot: only for the purpose of
set. It continuously updates the set of probabilities for the obtainingtheinitialpairwisematchingconfiguration,whilethe
mapping between the nodes in each graph to the virtual node post-procedure is immune from any pairwise matching solver.
set by optimizing the objective function until converge. The Fortheabovereason,inthefollowingplotswhenourmethods
final solution is consistent and discrete. Here we use the are tested, we use the naming convention “X-Y” to term the
exponentialobjectivefunctionasthesameformwithours,but different pairwise matching solvers used in two stages. While
inverseaffinitytocostKgacl =max(Kours)−Kours element- the compared other multi-graph matching methods are set to
wise, as GACL concerns with minimizing the cost function. always use the first stage solver as termed by “X”.
We set the parameter in GACL (refer to the original paper) For the second aspect, the proposed framework can either
β = .5, β = 20, α = 1.1 for the annealing outer loop use the proposed reference graph driven mechanism, or reuse
0 max
αβ → β ≤ β , and set the maximum iteration round any other multi-graph matching methods to generate an initial
t t+1 max
of the inner loop T =5 for cost efficiency. matching configuration X. As we will show later, these two
in
Modified hype-cube based common labeling (MHCL) approaches perform similarly in accuracy, while using the
the original hypercube framework is used in [33], [46] while reference graph based initialization only involves one-shot
specifically tested in the three-graph matching case, due to O(N3) times of multiplication of permutation matrix which
exponential memory overhead. The algorithm has two main can be very fast, compared with the spectral analysis based
steps: first, it performs a pairwise matching via existing method SYNC related to SVD computing with the complexity
pair solvers and averages either discrete assignment solu- ofO(N2n3)andthetediousthree-layerloopingprocedurefor
tions or probabilistic ones into an N-dimensional assign- GACLwhosetimecomplexitywouldbefurtherdiscussedlater
ment hypercube, as an extension to the permutation matrix in this section. Moreover, SYNC/MHCL are both unable to
for two graphs. Each element of the hypercube represents proceed their processing given the output of our algorithms,
the matching probability for the tuple of N nodes from because they require the input solutions are inconsistent and
N graphs respectively. Second, a post-binarization step is then transform the inconsistent solution to a consistent one.
applied to to obtain the consistent common labeling. While To make the plots in the rest of the paper more digestible,
the original hypercube data structure can only work with a themethodswiththeiracronymsarelistedforcross-reference.
very small number of graphs, not scalable to large N due • PAIRThebaselinewhichperformsexhaustivePAIRwise
to its memory space overhead is nN. Thus we modify the graph matching over the whole graph set. Note the
original method to a more light-weighted one including the consistency over pairwise matchings is not guaranteed;
following four steps: i) adopt a pairwise graph matching • AREF Abbreviated for Adaptive REFerence Graph Se-
solver to obtain all possible pair assignment solutions X0 ; lection, which consists of three steps: i) perform exhaus-
ij
ii) re-calculate the assignment matrix X by averaging the tive pairwise matching; ii) select the most “consistent”
ij
indirect mapping X =(cid:80)N X X /N; iii) use Hungarian
ij k=1 ik kj 4Therawconfigurationinducesafullyconnectedsupergraphforeachnode
method to binarize X ; iv) use the SYNC method (when
ij isagraph,andtheedgeisthepairwisematchingscore.AMaximumSpanning
N >n)orMaximumSpanningTree(MST)[62]overtheraw Tree(MST)canbefoundwhichspansanewconsistentconfiguration[47].
JUNCHIYANETAL.CONSISTENCY-DRIVENALTERNATINGOPTIMIZATIONFORMULTI-GRAPHMATCHING:AUNIFIEDAPPROACH 10
graph as the reference G according to the metric of then the cost of setting the reference graph and updating
r
graph-wise consistency; iii) use its pairwise matching order is O(N3n3) where O(n3) refers to the permutation
solutionsX withothergraphstospanthenewpairwise matrix multiplication that in fact can be signification speeded
rk
matching solutions by X =XTX ; up at the code level typically in a linear time w.r.t. n. In
ij ri rj
• FREF Abbreviated for Fixed (random) REFerence each iteration, as shown in formula (7), it adds up N − 2
Graph Selection, which replaces the first two steps of indirect pairwise affinities in addition with the direct affinity
AREF by randomly selecting a reference graph while objective. Building this whole objective function as the input
keeping the third step the same. It is very efficient since to RRWM costs O(Nn4) per iteration. The whole overhead
onlyalinearnumberofpairwisematchingsareperformed is N(τ +O(Nn4))+O(N3n3)+N2τ as our method
rrwm rrwm
w.r.t. the number of graphs N; stopswhenallassignmentmatrixareupdatedforonetime,i.e.
• GACL Abbreviated for Graduated Assignment based afterN−1roundofiterations.Thetimecomplexityestimation
Common Labeling as proposed in [31], [48]; for applying other pairwise graph matching solvers including
• SYNC Abbreviated for SYNChronization in [32]; IPFP/GAGM can also be derived in a similar manner.
• MHCL Abbreviated for Modified Hype-Cube based We further study the time complexity of other peer multi-
Common Labeling, which is introduced earlier in this graph matching methods. i) SYNC: it involves a one-shot
paper as an modified version of [33], [46] to handle the SVD to find the n leading eigenvectors on the grouped
memory bottleneck when more graphs are involved; assignment matrix of size nN × nN, whose complexity is
• FREF+/GACL+/SYNC+/MHCL+ Perform Alg.1/2 us- O(N2n3), and N iterations of Hungarian method with cost
ing FREF/GACL/SYNC/MHCL’s output as initial input. O(Nτ ) = O(Nn3). Thus the total cost is τ =
Hun sync
ThedifferencetoAlg.1/2issettingX∗ instep(8)asthe O(N2n3)+N2τ ; ii) MHCL: it performs the Hungarian
kr rrwm
solutions from other methods instead of pairwise match- method to binarize each matrix which is averaged by the pu-
ings.All“plus-sign”casesalsousetheconsistency-driven tative assignment matrices from N(N−1) pairwise matchings.
2
method to set the reference graph and updating order, ThecostisO(N2n3)+O(N3n3),wherethefirsttermiscon-
and this can be done in a free-rider fashion since other cernedwiththeHungarianmethod,andthesecondisrelatedto
methods also require computing pairwise matchings. It the averaging step involving N3 times of permutation matrix
is worth noting that FREF+ in fact corresponds to the multiplication, which can usually be optimized to O(N3n) at
method in our conference version [1] since either the thecodelevel.ThenSYNCisperformedtoobtainconsistency.
reference graph or the updating order is randomly set; Thus its time complexity is O(N3n3)+N2τ ; iii) GACL:
rrwm
• ALG12 One of our two methods, being Alg.1 when non- accordingtotheanalysisbytheauthorsin[31],itscomplexity
factorized graph matching solvers are adopted including is O(N2m2)+τ per iteration which is further surrounded
Sh
RRWM/IPFP/GAGM, or Alg.2 when factorized solvers by two layers of loops. Specifically, the outer loop concerns
are used like FGM. Note our method in fact can also be a continuation method which gradually increases the control
termedasAREF+basedontheabovenamingconvention. parameter β. The inner loop concerns iteratively updating the
This is because the result of AREF in fact initialize both assignmentmatrixtowardsthevirtualnodesetbygivenafixed
Alg.1andAlg.2instep(8).Forthisreason,intherelated β.ThusthetotalcomplexityisT O(N2m2)whereT isthe
ga ga
plots, AREF and ALG1 are in the same red color; total number of iterations over two layers of loops, which in
2
• ALG12* same as ALG12, but does not enforce the score- fact is significantly larger than N or n.
non-descending path selection strategy. Factorized multi-graph matchingDuetocurrentlyweare
unable to identify any other factorized multi-graph matching
method to our best knowledge, thus we only present the
D. Time complexity analysis
result for our method. Specifically, Besides the step for initial
Since the existing multi-graph matching methods SYNC pairwisematching,referencegraphandupdatingordersetting,
[32], MHCL [33], [46] start with the pairwise matching Alg.2 further involves i) factorizing the N(N − 1) pairs
solutions,andourmethodsalsousepairwisematchingsolvers of graphs that involves SVD decomposition at the cost of
for initialization and optimization. Thus we will consider O((n+m)3) per each; and ii) using the MFW method used
a specific pairwise matching solver and estimate the time in each iteration that includes calculating the gradient of the
complexity when applying it for the different methods. objective function ∇2 J at the cost of O((n+m)2) for
Xur ij
Non-factorized multi-graph matching Without loss of each pair (i,j), which is repeatedly performed over N −1
generality, we choose RRWM [15] as the pairwise matching terms (cid:80)N J + J in the objective function; iii)
f=1,(cid:54)=r,u fu ru
solver, and study the overall time complexity when applying performing the Hungarian method at the cost of O(n3) to
it in our non-factorized method for Alg.1. RRWM involves obtaintheoptimallinesearchdirectionYgiventhecalculated
an iterative procedure, in each iteration, the power iteration gradient ∇2 J; iv) the step-size search given the optimal
Xur
method [63] and Sinkhorn method [64] are employed. The searchdirectionthatisatthesamecostofcomputinggradient.
costoftheformerisO(m2)andthelatter,whichisdenotedas Thus the total complexity is O(N2(n+m)3)+T O(N(n+
fw
τ = O(n2) largely depends on the convergence speed, thus m)2 +n3)+N3n3 +N2τ where T is the number of
Sh rrwm fw
wehaveτ =T (O(m2)+τ )whereT referstothe iterations for the MFW method.
rrwm rrwm Sh rrwm
number of iteration rounds. For our method, first RRWM is Table III summarizes the time complexity by the general
performed on each pair of graphs whose cost is O(N2τ ), theoreticalanalysis.Ingeneral,Alg.1wouldberelativelymore
rrwm