Based Graph Matching Algorithm
Zhi-Yong Liu1, Hong Qiao1, and Lei Xu2
1 State Key Laboratory of Management and Control for Complex Systems, Institute of Automation, Chinese Academy of Sciences, Beijing, 100190, China
{zhiyong.liu,hong.qiao}@ia.ac.cn
2 Department of Computer Science and Engineering, Chinese University of Hong Kong, Shatin, Hong Kong
Abstract. In this paper we propose a regularized relaxation based graph matching algorithm. The graph matching problem is formulated as a con- strained convex quadratic program, by relaxing the permutation matrix to a doubly stochastic one. To gradually push the doubly stochastic ma- trix back to a permutation one, a simple weighted concave regular term is added to the convex objective function. The concave regular function is not a concave relaxation of the original matching problem. However, it is shown that such a simple concave regular term has a comparative performance as the concave relaxation of the PATH following algorithm, which works only on undirected graphs. A concave-convex procedure (CCCP) together with the Frank-Wolfe algorithm is adopted to solve the matching problem, and some experimental results witness the state- of-art performance of the proposed algorithm.
Keywords: graph matching, CCCP, Frank-Wolfe algorithm, convex relaxation.
1 Introduction
Graph matching plays a central role in many graph based techniques. For in- stance, graph is frequently used as the structural representation of objects in computer vision and pattern recognition, and consequently the graph matching algorithm is usually used to solve the object matching problem [1]. Graph match- ing involves identifying each vertices pair between two graphs in some optimal way, or inherently finding a good permutation matrix between the two adjacency matrices of both graphs. The problem is in nature a NP-hard combinatorial opti- mization problem with factorial complexity, except for some graphs with special structure, such as the planar graphs which has shown to be of polynomial com- plexity [2]. Consequently, an exhaustive enumeration based matching algorithm is computationally prohibited, except for some small-scale problems.
To make the problem computationally tractable, many approximate approaches have been proposed, to seek an acceptable trade-off between the complexity and
Y. Zhang et al. (Eds.): IScIDE 2011, LNCS 7202, pp. 9–16, 2012.
Springer-Verlag Berlin Heidelberg 2012c
matching accuracy. As summarized in [3], approximate matching algorithms can be roughly categorized into three groups, tree search based methods, spectral meth- ods and continuous optimization (relaxation techniques). Here, we concentrate on the relaxation techniques, which involve relaxing the combinatorial matching problem to a continuous one [4–6]. The key point lies in the fact that optimiza- tion over a continuous set is usually easier to be approximated than its discrete counterpart. Two well-known relaxation techniques in literature include gradu- ated assignment [5] and PATH following algorithm [6], both of which work directly on adjacency matrices.
Given two graphs GD= (VD, ED) and GM = (VM, EM) where V and E re- spectively denote the sets of vertices and edges and by formulating the matching problem as
P = arg min
P ∈Pf(AD, AP (M)), AP (M)= P AMPT, (1) where A denotes adjacency matrix, P is a permutation matrix, AP (M) denotes the adjacency matrix of the permutated graph of GM by P ,P denotes the set of permutation matrix, and f (·) is a cost function which usually takes a convex quadratic form, both gradual assignment and PATH following algorithms relax domain of the objective from set of permutation matrix P to it convex hull, i.e., set of doubly stochastic matrixD. In this way, eqn. (1) becomes a convex quadratic program. A soft assignment schema controlled by a parameter was introduced by the graduated assignment algorithm to control the non-convexity of the problem, and a double normalization is used to constrain the matrix a doubly stochastic one. As the parameter increases to large enough, a permutation matrix is expected to be obtained though usually a clean-up step is further needed.
Different from the gradual assignment algorithm, the PATH following algo- rithm introduces a weighted concave relaxation to control the non-convexity of the objective. Specifically, by adopting the square of Frobenius matrix norm as the cost, the objective in eqn. (1) is firstly relaxed as follows
fv(P ) = vec(P )TQvec(P ), P ∈ D, (2) where vec(P ) is the column-wise vector representation of P , and Q = (I⊗ AD− ATM⊗ I)T(I⊗ AD− ATM⊗ I) ∈ RN2×N2 is a symmetric definite positive matrix.
The concave relaxation introduced by PATH following algorithm is given by fc(P ) =−tr(P ) − 2vec(P )T(LTM ⊗ LTD)vec(P ), P ∈ D, (3) where ij = (DM(i, i)− DD(j, j))2, with D and L denoting the degree and Laplacian matrices of the graph, respectively. The concave relaxation holds the same global minimum as the original matching problem, that is, minP ∈Pfc(P ) = minP ∈Pfv(P ). Based on the convex and concave terms above, the objective function of the PATH following algorithm is given by
fpath(γ, P ) = γfv(P ) + (1− γ)fc(P ), P ∈ D, (4)
where γ∈ [0, 1] controls the non-convexity of the objective: A large γ means that fpath(P, γ) tends to be convex; by contrast, a small γ makes fpath(P, γ) concave.
Thus, by gradually decreasing γ for 1 to 0, the objective becomes finally a concave one, and its minimization results in a permutation matrix.
However, the PATH following algorithm cannot be used to solve the matching problem between directed graphs because the term in eqn. (3) cannot guarantee to be concave for directed graphs. In this paper we introduce a much simpler con- cave term which can be applied on both directed and undirected graphs. Though the simple concave term is not a relaxation of the original matching problem, it is shown that it has a comparative performance as eqn. (3) on matching accuracy.
Section 2 is devoted to the proposed method, some experimental illustrations and discussions are given in section 3, and finally section 4 concludes the paper.
2 Proposed Method
The objective function for the graph matching problem is firstly proposed, and then the CCCP together with Frank-Wolfe algorithm is proposed to minimize the objective.
2.1 Objective Function
The proposed objective function takes a similar form as eqn. (4), with the same convex relaxation but with a different concave term. To make the algorithm applicable also for matching problems on directed graphs, we propose to use the following simple concave term,
fc(P ) =−vec(P )Tvec(P ), P ∈ D. (5) Then, the graph matching problem is formulated as follows,
min fγ(P ) = γvec(P )TQvec(P ) − (1 − γ)vec(P )Tvec(P )
s.t. P ∈ D. (6)
It is obvious that minimization of the concave term given by eqn. (5) results in an extreme point ofD, i.e., a permutation matrix. Thus, by gradually decreasing γ from 1 to 0, minimization of the objective will make P gradually converge to a permutation matrix. At the beginning when γ = 1, the objective degenerates to the convex relaxation, whose global minimization denoted by Pvcan be obtained by such as the Frank-Wolfe algorithm. Actually, a permutation matrix can be directly obtained by an optimal assignment procedure which casts the doubly stochastic matrix Pv to a permutation matrix Ppvia
Pp= arg max
P ∈PPv, P , (7)
where· denotes an inner product. The assignment can be solved by the Hungar- ian algorithm [7], with the computational complexity O(N3). Such a hard-cut
operation based graph matching algorithm, named QCV (quadratic convex), may however bring a big error in the result, as to be witnessed by the experi- mental results in section 3. By contrast, as γ gradually decreases, P is gradually pushed away from Pv in such a way that update of P is guided to approach a permutation matrix with smaller matching error. This point can be intuitively understood in the following way. During the convergence process the update di- rection of P comprises two parts, gv(P ) and gc(P ), the directions provided by the convex and concave terms, respectively. Guidance from gv(P ) is to minimize the increase of the convex term, which, if can be globally minimized during the whole process, is equal to the difference between the best matching error and the global minimization of the convex relaxation gotten by Pv. On the other hand, gc(P ) provides no informative search direction since any permutation matrix gives the same global minimum for the concave term. Thus, in the global mini- mization sense it is under the guidance of gc(P ) that P is expected to approach a permutation matrix with a relatively small matching error.
In contrast to the above simple concave term, the concave relaxation given by eqn. (3) also provides a search direction gc(P ) which also provides a meaningful guidance for the update of P in the global optimization sense. However, starting from Pv, the search direction gc(P ) provided by the concave relaxation is the same as gv(P ), i.e., the direction from Pvto the global optimal point. Therefore, gc(P ) is somewhat redundant to gv(P ), and gc(P ) just strengthens gv(P ) but does not provide additional useful guidance. This is to some extent confirmed by the experimental comparisons in section 3 which witnesses that, on matching undirected graphs, the two concave terms have a comparative performance on accuracy (see Tab. 1 and Figs. 1 and 2 for details).
2.2 Algorithm
For each fixed γ, the objective function given by eqn. (6) is a constrained quadratic program which is generally neither convex nor concave, to solve which some local search technique seems to be unavoidable. Here, we firstly utilize the concave-convex procedure (CCCP) to decompose the objective function into a sequential constrained convex quadratic program, which is then solved by the Frank-Wolfe algorithm.
The CCCP algorithm consists of sequentially minimizing the following con- strained convex function,
fk+1(Pk+1) = fv(Pk+1) + vec(Pk+1)T∇fc(Pk), Pk+1∈ D, (8) where Pk+1 denotes the P to be found in step k, and fv, fc the convex and concave terms respectively. Since ∇fc(Pk) = −2(1 − γ)vec(Pk) is a constant with respective to Pk+1, eqn. (8) is formulated as the following constrained convex quadratic program:
minfcccp(P ) = γvec(P )TQvec(P ) − 2(1 − γ)vec(Pk)Tvec(P ), s.t. P ∈ D. (9)
The Frank-Wolfe algorithm, as a reduced gradient method, is adopted to solve the above constrained convex quadratic program. Specifically, it comprises the following four steps.
step 1: Initialize P0= P∗ and let t = 0, where P∗ denotes the result obtained by the previous CCCP loop.
step 2: Find an extreme point Xt(a permutation matrix) ofD by solving the linear program
min∇fcccp(Pt), Xt, s.t. Xt∈ D, (10) where∇fcccp(P ) is given by
∇fcccp(P ) = 2γ(ADADP − ADP AM − ADP AM + P AMAM)− 2(1 − γ)P∗. (11) step 3: Find a step size α ∈ [0, 1] to minimize fcccp(Pt+ α(Xt− Pt)), and
update Pt+1= Pt+ α(Xt− Pt).
step 4: If|∇fcccp(Pt), Xt− Pt| < ε|fcccp(Pt) +∇fcccp(Pt), Xt− Pt| where ε is a small positive constant, return Pt+1. Otherwise, let t = t + 1 and go back to step 2.
In the algorithm, the linear program in step 2 can be solved by the Hungarian algorithm with a complexity O(N3), and the line search can be efficiently im- plemented by for instance a backtracking algorithm [8]. The stopping criterion in step 4 is applicable, thanks to the convexity of the objective function fcccp.
Finally, the graph matching algorithm is summarized by Algorithm 2.1.
Algorithm 2.1. GraphMatching(AD, AM)
Pn← 1N×N/N, n ← 0, γ ← 1 whileγ ≥ 0&P /∈ P
do Pn+1← CCCP(Pn, γ), γ ← γ − δγ, n ← n + 1 return (Pn)
Storage complexity of the algorithm is O(N2) and since the computational complexities of the Hungarian algorithm and matrix multiplication involved in the algorithm are both O(N3), the computational complexity of the algorithm is roughly O(N3).
3 Experimental Illustrations
Three experiments on synthetic data are conducted to evaluate the proposed algorithm. Five algorithms, including Umeyama’s algorithm (U for short) [9], PATH following algorithm [6] (on undirected graphs only), QCV algorithm given by eqn. (7), graduated assignment algorithm [5] (GA for short), and the proposed algorithm are experimentally compared on both undirected and directed graphs.
For problems with small scale (N = 8 for instance), an exhaustive search is used to get the ground true solution. All of the algorithms were implemented by Matlab 2009b, with a MEX function to implement the Hungarian algorithm.
The first experiment is to simulate the scenario of graph matching without any prior. In the experiment, 100 pairs of graphs with size N = 8 are randomly generated by the following procedure: For each entry Aij(Aji= Aijin the case of undirected graph) randomly generate a uniformly distributed number r∈ [0, 1]; if r > 0.5 (meaning that sparsity of the graph is around 0.5), randomly generate its weight Aij = w∈ [0, 1], or otherwise Aij = 0. The first experimental results are listed in Tab. 1, from which it is witnessed that the PATH (on undirected graphs only) and our algorithms have a much better performance on accuracy than the U, QCV and GA algorithms. It is also observed that PATH and our algorithm have a comparative performance on accuracy, as echoed by the discussions in section 2.1.
Table 1. Comparative experimental results on four types of graphs with N = 8, summarized from 100 random runs
graph types error OPT U QCV GA PATH Ours mean 3.5473 8.7955 6.3547 6.3483 4.1323 4.0591 undirected graph models
std 1.0678 3.3016 2.5465 2.1100 1.4218 1.3254 mean 5.5269 11.0624 7.9136 9.6330 - 6.3928 directed graph models
std 0.8921 2.2747 1.8446 2.1368 - 1.2297
The second experiment is to evaluate the noise robustness of the algorithms.
In the experiment, the second graph in a graph pair was generated based on the first one by adding some noises which is controlled by a noise level. Specifically, given a noise level δ ∈ [0, 1] and a randomly generated adjacency matrix AD, AM is gotten by the following steps.
1. Set AM ← AD, and for each (AM)ij, randomly generate two variables r1 and r2∈ [0, 1].
2. If (AM)ij> 0: if r1< δ, (AM)ij← 0; or otherwise, (AM)ij ← (AM)ij+ δr2. 3. If (AM)ij = 0: if r1< δ, (AM)ij← r2.
4. Randomly generate a permutation matrix P , and set AM ← P AMPT. The noise level increases from 0 to 1 by a step size 0.1. On each noise level, 100 graph pairs with N = 8 are generated according to the above process. The experimental results on the two types of graphs are shown in Fig. 1. Similar to the first experimental result, the PATH and our algorithm outperform significantly the other three algorithms.
The third experiment is to evaluate the scalability of the five algorithms with respective to graph size, on both accuracy and complexity. In the experiment, 10 groups of graph pairs with different sizes are included for comparison, with the size increasing from 5 to 50 by a step size 5. In each group, 50 graph pairs are generated in the same way as the second experiment with a noise level 0.1.
0 0.2 0.4 0.6 0.8 1 0
2 4 6 8 10 12 14
noise level
matching error
on undirected graph models OPT
U QCV GA PATH Ours
0 0.2 0.4 0.6 0.8 1
0 2 4 6 8 10 12 14 16 18
noise level
matching error
on directed graph models OPT
U QCV GA Ours
Fig. 1. Changes of the matching error with respective to noise levels, summarized from 100 random runs
On accuracy, similar experimental results as the above two are obtained on all of the 20 groups of graph pairs, as witnessed by Figs. 2. The time-cost of the five algorithms are shown in Fig. 3, in which the slopes of the five curves corresponding to U, QCV, GA, PATH, and our algorithm are around 2.412, 2, 633, 3.135, 2.762, and 2.737, respectively, which imply that the GA suf- fers the biggest complexity and U is the simplest one.
10 20 30 40 50
0 50 100 150 200 250
graph size
matching error
on undirected graph models U
QCV GA PATH Ours
10 20 30 40 50
0 50 100 150 200 250
graph size
matching error
on directed graph models U QCV GA Ours
Fig. 2. Changes of the matching error with respective to graph sizes, summa- rized from 50 random runs on each size
0.5 1 1.5 2
−4
−3.5
−3
−2.5
−2
−1.5
−1
−0.5 0 0.5 1
graph size (log10N) time−cost (log10Sec.)
on undirected graph models U
QCV GA PATH Ours
0.5 1 1.5 2
−4
−3.5
−3
−2.5
−2
−1.5
−1
−0.5 0 0.5 1
graph size (log10N) time−cost (log10Sec.)
on undirected graph models U
QCV GA Ours
Fig. 3. Time-cost of the five algorithms with respective to the graph size, summa- rized from 50 random runs on each size
4 Conclusions
In this paper we showed that a simple weighted concave function has a compar- ative performance to the concave relaxation for the graph matching problem.
The point is that the simple concave function can be utilized on matching differ- ent types of graphs, but by contrast, the concave relaxation is usually difficult to find. Actually, to the best of our knowledge, only the concave relaxation on undirected graphs without self-loops has been reported in literature [6]1.
1 Just before finishing the final version of this paper, we reported two type of concave relaxations for directed graph models, see [10, 11].
The CCCP together with the Frank-Wolfe algorithm is then proposed to solve the matching problem. On four different types of graphs, our algorithm showed a state-of-art performance on accuracy.
Acknowledgement. The work described in the paper was supported by a grand from NSFC (grant 60975002), and the National Basic Research Program of China (973 Program) (grant 2009CB825404).
References
1. Eshera, M.A., Fu, K.S.: An image understanding system using attributed sym- bolic representation and inexact graph-matching. IEEE Transactions on Pattern Analysis and Machine Intelligence 8(5), 604–618 (1986)
2. Hopcroft, J.E., Wong, J.K.: Linear time algorithm for isomorphism of planar graphs (preliminary report). In: Proceedings of the Sixth Annual ACM Symposium on Theory of Computing, STOC 1974, pp. 172–184. ACM, New York (1974)
3. Conte, D., Foggia, P., Sansone, C., Vento, M.: Thirty years of graph matching in pattern recognition. International Journal of Pattern Recognition and Artificial Intelligence 18(3), 265–298 (2004)
4. Xu, L., Oja, E.: Improved Simulated Annealing, Boltzmann Machine and At- tributed Graph Matching. In: Almeida, L.B., Wellekens, C.J. (eds.) EURASIP 1990. LNCS, vol. 412, pp. 151–160. Springer, Heidelberg (1990)
5. Gold, S., Rangarajan, A.: A graduated assignment algorithm for graph matching.
IEEE Transactions on Pattern Analysis and Machine Intelligence 18(4), 377–388 (1996)
6. Zaslavskiy, M., Bach, F., Vert, J.P.: A path following algorithm for the graph matching problem. IEEE Transactions on Pattern Analysis and Machine Intelli- gence 31(12), 2227–2242 (2009)
7. Kuhn, H.W.: The hungarian method for the assignment problem. Naval Research Logistics Quarterly 2(1-2), 83–97 (1955)
8. Boyd, S., Vandenberghe, L.: Convex Optimization. Cambridge University Press (2004)
9. Umeyama, S.: An eigendecomposition approach to weighted graph matching prob- lems. IEEE Transactions on Pattern Analysis and Machine Intelligence 10(5), 695–
703 (1988)
10. Liu, Z.Y., Qiao, H., Xu, L.: An extended path following algorithm for graph match- ing problem. IEEE Transactions on Pattern Analysis and Machine Intelligence (to appear, 2012)
11. Liu, Z.Y.: Graph matching: a new concave relaxation fuction and algorithm. Au- tomatica Sinica (to appear, 2012)
