Dynamic Programming Bipartite Belief Propagation For Hyper Graph Matching

1 Introduction

The feature matching problem aims to establish consistent correspondences between two feature sets, and is a fundamental subroutine in computer vision and pattern recognition [Leordeanu and Hebert, 2005; Torresani et al., 2008; Li et al., 2010]. When establishing correspondences between features, it is important to consider not only local similarity between features, but also structured similarity between sets of features. Graphical models provide an elegant framework to encode feature matching problems, as they are naturally able to express structured relationships between sets of variables.

However, graphical model based methods [Zhang et al., 2016; Torresani et al., 2008; Ahn et al., 2015] are typically limited to pairwise matching, which can be used to express (for example) matches between rigid structures. Unfortunately, pairwise relations are insufficient to encode rich geometrical structures [Lee et al., 2011]; thus higher-order information should be included and the problem formulated as one of hyper graph matching [Yan et al., 2015; Nguyen et al., 2015; Lee et al., 2011; Duchenne et al., 2009].

Such a formulation can be more robust to noise and outliers, though hyper graph matching is less well studied than pairwise matching, due to the increased difficulty of optimization.

Among hyper graph matching methods, most recent ones fall into one of two categories. In the first category, the correspondence vector (which encodes the matching) is relaxed to be continuous, and then post-processing (such as the Munkres/Hungarian algorithm) is used to discretize the correspondence vector; Hyper Reweighted Random Walks Methods (HRRWM) [Lee et al., 2011] belong to this category. The second category, which includes Discrete Hyper Graph Matching [Yan et al., 2015] and Tensor Block Coordinate Ascent Graph Matching (TBCAGM) [Nguyen et al., 2015], relaxes the matching objective as a multi-linear or multi-quadratic problem, and then alternating optimization methods are used to optimize over the relaxed objective. Due to the relaxation such alternating optimization methods may return multiple inconsistent correspondence vectors, though in practice their performance is state-of-the-art. Regularizers can be added to the objective to encourage the multiple correspondence vectors to be consistent, however this may result in a looser relaxation. Another issue for multi-linear or multi-quadratic relaxations is their non-convexity; due to non-convexity, alternating optimization methods may easily get stuck in local optima.

In this paper, we formulate the hyper graph matching problem as a constrained higher order MAP inference problem in a graphical model. To solve this we propose a convex linear programming (LP) relaxation. In the primal form of the relaxation, rather than introduce several global correspondence vectors, we introduce only one global correspondence, but several local correspondence vectors, which represent correspondences only within a small sub-graph (e.g. an edge, a triplet or other small sub-graph). In the dual formulation, the hyper graph matching problem is decomposed as a linear bipartite matching problem (which can be efficiently solved by the Munkres/Hungarian algorithm [Munkres, 1957]), and several non-linear matching problems, which we term belief propagation sub-problems, with relatively small scale. These sub-problems are solved iteratively, and during optimization the global correspondence vector is encouraged to be consistent with the local correspondences. If the global correspondence vector is consistent with all local correspondences, our algorithm is guaranteed to provide an exact solution.

Although some previous work has used similar formulations for pairwise matching problems [Torresani et al., 2008; Zhang et al., 2016], the exhaustive search procedures involved when solving belief propagation sub-problems are computationally expensive when applied to hyper graph matching. As a response, by exploiting the internal structure of hyper graph matching problems, we further decompose the belief propagation sub-problems as several overlapping subproblems with optimal substructure. This allows efficient dynamic programming methods to be developed to solve these problems. We title the proposed methods Dynamic Programming Bipartite Belief Propagation. Experiments show that the proposed methods outperform the current state-of-the-art in hyper graph matching problems.

2 Notation and ProblemFormulation

A hyper graph H = (V, C) consists of nodes v ∈ V and hyper arcs c ∈ C, where a hyper arc c is a subset of V. The order of a hyper arc is its cardinality |c| and the order of a hyper graph is the largest order among all hyper arcs.

We seek a correspondence between two hyper graphs: the model graph HP = (VP , CP ) and the data graph HQ = (VQ, CQ), and their associated attributes AP and AQ. For any hyper arc c from CP or CQ, there is a corresponding attribute, denoted by aPc or aQc . Attributes could encode the orientation of a keypoint, the length of an edge, the angle between two edges, etc..

The goal of hyper graph matching is to find the optimal correspondence between VP and VQ. Here to simplify the derivation, we assume that V = VP = VQ = {1, 2, · · · , n}, though this can easily be relaxed (e.g. to handle outliers) by adding 'dummy' nodes in either VP or VQ. Given a correspondence vector yi = j (meaning that the ith node in VP corresponds to the j th node in VQ), finding an optimal oneto-one matching between VP and VQ can be formulated as a constrained MAP inference problem in a graphical model:

where 1(S) is an indicator function (i.e. 1 if S is true and 0 otherwise). θi(yi) and θc(yc) are value functions (or 'potentials') of the form

where ϕ are real-valued functions, which measure the similarity between unary or higher order attributes. In practice, the similarity function are often defined as Duchenne et al.; Nguyen et al. [2009; 2015]

where γ1,|c| ⩾ 0 and γ2,|c| ⩾ 0 are problem-specific parameters for order |c| hyper arcs.

3 LP Relaxation and Its Dual

The problem (1) is NP-hard in general, thus relaxations are required. The approaches of Zhang et al. [2016] and Torresani et al. [2008] are restricted to pairwise problems. Thus we proposed a linear programming (LP) relaxation for both pairwise and hyper graph matching problems as follows:

where 〈f, g〉 = ∑x f(x)g(x) denotes the inner product of two functions. In the LP relaxation, µ = [µi(yi)]i∈V,yi∈[n] is the global correspondence vector, and µi(yi = j) = 1 means that the ith node in VP corresponds to the j th node in VQ. Each µc(yc) is a local correspondence vector, and µc(yc = [ai]i∈c) = 1 means that the ith node in VP corresponds to the (ai)th node in VQ for all i ∈ c.

It is well-known that 'off-the-shelf' solvers such as CPLEX are slow for LP problems like (4) [Werner, 2007]. Thus we use dual block coordinate descent based belief propagation to solve such problems efficiently. First we introduce dual variables (shown in Table 1), then by standard Lagrangian duality the dual objective is as follows,

In order to satisfy the one-to-one matching constraint, the feasible set Yc is

The dual variables λc→i(yi) are referred to as messages, and u, v are referred to as matching variables as in previous work [Zhang et al., 2016]. Similar to Zhang et al. [2016] for pairwise problems, we alternately optimize over u, v and λ. When optimizing u and v, we fix all λ's. When optimizing λ, in each step we pick a particular c from CP , and allow those λc→i(yi), i ∈ c to vary. Then the dual problems (5) can be decomposed as the following two groups of sub-problems:

where λc denotes [λc→i(yi)]i∈c, and the augmented potential ϑi(yi) is defined as

The matching sub-problem (7) can be solved by the Hungarian algorithm [Zhang et al., 2016]. However, the belief propagation sub-problems are not easy to solve in general. For general higher order potentials, the time complexity for solving (8) is O(n|c| ), which increases very quickly as the order of hyper arcs grows. However, for hyper-graph matching problems the higher order potentials are usually sparse [Yan et al., 2015; Nguyen et al., 2015; Lee et al., 2011]. Thus we will use the sparsity to derive an efficient solver.

Remarks For sparse potentials whose non-zero entries are positive, Rother et al. [2009]'s and Potetz and Lee [2008]'s approaches can be applied to derive an efficient message passing procedure. However, these methods are not compatible with one-to-one matching constraints which must be satisfied in matching problems.

4 Dynamic Programming Approaches for the Belief Propagation Sub-problem

It is intractable to compute affinities between two hyper graphs by considering all possible mappings between hyper arcs. Even for pairwise matching problems, the number of terms required to evaluate all possible pairs of connections is already n4 , which is quite expensive for large-scale problems. Thus we seek methods that can be used to sparsify the potentials. Here we derive dynamic programming approaches for the belief propagation sub-problem that exploit the sparsity of higher-order potentials.

4.1 Decomposition of the Belief Propagation Sub-problem

To obtain our dynamic programming belief propagation scheme, we start at a closed-form solution of the belief propagation sub-problem. To simplify our derivation we define

then we can give the closed-form solution of (8) as

where the domain Yyi c is defined as

The closed-form solution (11) can be computed by exhaustive search. However, the time complexity required would be O(n|c| ), which is not scalable for hyper graph matching problems. Thus we decompose the most expensive part (the max-marginal) into several overlapping problems with optimal substructure.

Proposition 1. The max-marginal procedure in (11) can be be performed as follows

The second sub-problem (last row of (13)) only considers non-zero entries of θc(yc), which can be done efficiently via exhaustive search. The solution of the first term (second last row of (13)) can be reformulated as a linear bipartite matching problem as follows.

Proposition 2. Let y(l) i′ denote the entries corresponding to the lth largest ζi′(yi′), and let

then we have that

By Proposition 2, the first term in (13) can be efficiently computed via |c| partial sort operations, and solving a bipartite matching problem whose cost matrix size is |c|2 . While the size of |c| is very small (usually 3 or 4) in most higher order matching problems [Yan et al., 2015; Nguyen et al., 2015] the bipartite matching problem can be solved via exhaustive search with complexity O(|c|!). For large cliques, we can use Hungarian methods to further accelerate the procedure.

4.2 The Partial Reparametrization Formulation

The terms bi(yi) and bc(yc) are in fact reparametrizations of the augmented potentials, that is ∀y

Thus in our implementation, we only store bi(yi), θc(yc), all messages λ, and matching variables u, v. Using this formulation, the augmented potentials ϑi(yi) and θi(yi) can be recovered easily. Another benefit of this formulation (rather

Algorithm 1: The DPMU Procedure

input : A clique c, potential bi(yi), i ∈ c, messages λc→i(yi). output: b∗i (yi), λ∗c→i(yi), i ∈ c 1 for i ∈ c do 2 ζi(yi) = bi(yi) - λc→i(yi); 3 Compute ξi(yi), i ∈ c by Proposition 2; 4 with ζi(yi) as input; 5 for i ∈ c do b∗i (yi) ← ξi(yi); 6 for yc ∈ non-zero entries of θc(yc ) do 7 fv = θc(yc ) + ∑i∈c ζi(yi); 8 for i ∈ c do 9 if fv > b∗i (yi) then b∗i (yi) ← fv; 10 for i ∈ c do 11 b∗i (yi) ← 1 |c| b∗i (yi); 12 λ∗c→i(yi) ← λc→i(yi) - bi(yi) + b∗i (yi);

than storing all potentials) is that significant computation can be avoided. Recall the definition of bi(yi) in (10); if we only store messages, then several messages have to be added together every time we compute bi(yi). However the optimal b∗i (yi) can be computed directly from max-marginals by the following proposition.

Proposition 3. The optimal b∗i (yi) can be computed by

By Proposition 3 and (11), the optimal messages can be computed as

The above results, Propositions 2 and 1 together result in our Dynamic Programming Message Updating (DPMU) algorithm in Algorithm 1. Putting DPMU into a coordinate descent framework results in our Dynamic Programming Bipartite Belief Propagation (DPBBP) algorithm in Algorithm 2.

Decoding Our Dynamic Programming Bipartite Belief Propagation methods are dual based methods, thus a feasible integer solution is required. We consider two strategies to decode an integer solution. In the first strategy, integer solutions are decoded by the bipartite matching. This strategy is natural since in each iteration of belief propagation, a bipartite matching procedure is used to optimize over the dual variables u and v, and they provide a feasible primal integer solution at the same time. In the second strategy, we use a multi-linear relaxation method to further refine our results. Although such methods can easily get stuck at non-optimal points, our algorithm usually provides near optimal solutions in practice, so that using multi-linear relaxations to refine the results usually yields satisfactory results. Multi-linear relaxations are more expensive than bipartite matching, thus we only perform this refinement procedure at the last step.

Algorithm 2: The Dynamic Programming Bipartite BP

input : Potentials θi(yi), i ∈ V, θc(yc ),c ∈ CP ; MaxIter; threshold ϵ1 and ϵ2. output: y⋆ . 1 fmax = -∞, u = 0, v = 0 ; 2 bi(yi) = θi(yi), ∀i ∈ V, yi ∈ [n]; 3 for k ∈ {1, 2, . . . ,MaxIter} do 4 for c ∈ CP do 5 Compute optimal λ∗c→i(yi), b∗i (yi), i ∈ c by Algorithm 1; 6 Updating λc→i(yi), bi(yi), i ∈ c to λ∗c→i(yi), b∗i (yi), i ∈ c; 7 Compute optimal u⋆ , v⋆ for sub-problem (7), and decoding a candidates ˆy by the Hungarian algorithm; 8 Update [u, v] to [u⋆ , v⋆ ]; 9 fk = ∑i∈V θi(ˆ yi) + ∑c∈CP θc(ˆyc ); 10 If fk > fmax then fmax = fk, y⋆ = ˆy; 11 gk ← current dual objective of (5); 12 if |fmax - gk| < ϵ1 or |gk - gk-1| < ϵ2 then break;

4.3 Analysis

In this section, we analyse the convergence, time complexity and accuracy of the proposed method. The convergence of the proposed methods is guaranteed by the following proposition.

Proposition 4. For arbitrary bounded input, the Dynamic Programming Bipartite Belief Propagation produces a monotonically decreasing dual objective sequence.

Since the MAP value is a natural lower bound of the dual objective, the proposed methods must converge to a fixed dual objective. At the fixed point of the proposed algorithm, the local and global correspondence vectors have the following weak consistency properties:

Proposition 5 (Weak Consistency). Assume that Dynamic Programming Bipartite Belief Propagation converges to a fixed point; let ˆy be the optimal assignment produced by linear bipartite matching in (7), and ∀c ∈ CP let

then we have that ∀c ∈ CP ,

By the weak consistency properties, we can see that the local and global correspondence vectors are encouraged to attain the same optimal objective value on unary terms. By this property, we immediately obtain the condition for global optimality.

Proposition 6. Assume that Dynamic Programming Bipartite Belief Propagation converges to a fixed point; if ∀i ∈ V, there some ˆ yi s.t. ∀¯ yi ̸= ˆ yi, b(ˆ yi) > b(¯ yi), then Dynamic Programming Bipartite Belief Propagation provides a globally optimal solution of the problem (1).

Time Complexity In each iteration, we must solve | CP | belief propagation sub-problems and one linear bipartite matching problem. The time complexity of solving a bipartite matching problem is O(n3 ). The time complexity for solving belief propagation sub-problems by Algorithm 1 is O(|c|n log |c| + |c|! + |c||NNZ(θc(yc))|), where we use NNZ(θc(yc)) to denote the set of non-zero entries of θc(yc). On the other hand the time complexity for solving belief propagation sub-problems by exhaustive search is O(n|c| ). Although both expressions grow exponentially in |c|, the proposed methods are still significantly faster than exhaustive search for sparse potentials. This is because in practice the value |c| is usually very small (often 2 or 3), so that |c| and even |c|! can be viewed as (small) constants. Treating |c| as a constant, the time complexity for each iteration of our algorithm is O(n3 + | CP |n+ ∑c∈CP |NNZ(θc(yc))|).

5 Experiments

In this section, we apply Dynamic Programming Bipartite Belief Propagation (DP-BBP) to several hyper graph matching tasks. Our method is applied with two decoding schemes. In the first decoding scheme, we simply use solutions for bipartite matching as candidates, while in the second one, we use multiple linear relaxation based methods to further refine the decoding results. The first algorithm is referred to as DP-BBP, and the second is referred to as DP-BBP-MLR. Our method is compared with several existing popular and state-of-the-art hyper graph matching algorithms, including:

• The Tensor Block Coordinate Descent Method for multilinear relaxation ("BCA") [Nguyen et al., 2015];
• The Tensor Block Coordinate Descent Method for multiquadratic relaxation with IPFP [Leordeanu et al., 2009] as the quadratic solver ("BCA-IPFP");
• The Tensor Block Coordinate Descent Method for multiquadratic relaxation with Max Pooling Matching [Cho et al., 2014] as the quadratic solver ("BCA-MP");
• The Probabilistic Hyper Graph Matching method ("HGM") [Zass and Shashua, 2008];
• The Tensor Matching method ("TM") [Duchenne et al., 2009];
• The hyper graph matching via reweighted random walking method ("HRRWM") [Lee et al., 2011].

Two state-of-the-art quadratic graph matching algorithms, Factorized Graph Matching (FGM) [Zhou and De la Torre, 2012] and Bipartite Belief Propagation (BBP) [Zhang et al., 2016], are also included as baselines. We mainly compare the accuracy and objective value for different algorithms.

Implementation details In our experiments, the proposed algorithms DP-BBP and DP-BBP-MLR are implemented in c++ and python. As in previous work [Zhang et al., 2016], when there is a gap between the dual and decoded primal, a most-fractional-first branch-and-bound strategy is used to tighten the gap. We run at most 100 branch-and-bound iterations, and in each branch-and-bound, we run at most 10 iterations of DP-BBP or DP-BBP-MLR. All other methods are based on publicly available implementations.

Higher Order Potential Computation In our experiments, we mainly consider third-order similarities, which are invariant to scale and rotation. Given two sets of points P = [pi]ni=1 and Q = [qi]ni=1 , we use Delaunay triangulation to construct the set of hyper arcs CP and CQ. The higher order potentials are computed as

where αic (αiyc ) is the ith angle of the triangle c(yc), and the distance between two angles d(αic , αiyc ) is computed as

The pairwise potentials for pairwise methods are as follows,

where dij (dyiyj ) are Euclidean distances between nodes i and j (yi and yj), and τij,yiyj is a user-specified parameter.

5.1 CMUHouse Dataset

The CMU house dataset has been widely used in previous work to evaluate matching algorithms [Nguyen et al., 2015]. The house sequence consists of 111 frames of a toy house captured from different view points. In each image there are 30 manually marked landmark points with known correspondences. We match all images spaced by 10, 20, . . . , 90 frames and compute the average performance per separation gap. In this experiment, unary terms are set to zero [Nguyen et al., 2015]. Again for each image we use Delaunay triangulation to connect landmark points. Each triplet in the connected graph becomes a hyper arc in the hyper graph (or each edge in the case of pairwise matching methods). Higher order potentials are computed as in (17) and pairwise potentials are computed as in (18). For hyper graph matching methods, we only use higher order potentials and we use pairwise potentials only for pairwise matching methods, where we set the parameter τij,yiyj = 1/2500 [Leordeanu et al., 2009].

We consider two different settings. In the first setting, all landmark points are used to construct the hyper graph. In the second setting, for each image pair, we randomly remove 5 landmarks from each image, which results in outliers (i.e. landmarks that do not have any correspondences). Then for landmarks in the second image in the pair, a scale transformation with scale parameter 4 and a rotation transformation with angle 3π/4 are applied.

Results are shown in Figure 1. In the first setting, for small image separation all methods find the exact matching, while the accuracy of BCA, BCA-MP, HGM and TM drops as the separation increases. In the second setting, pairwise matching methods fail to obtain satisfactory results as they are sensitive to scale variation and rotation. For hyper graph matching methods, for image separation less than 80 frames our algorithm attains similar (or slightly better) accuracy compared to the best results among BCA-MP, BCA-IPFP and HRRWM. For the highest separation our methods' accuracy is more than

10% higher compared to these methods.

5.2 The Cars andMotorbikes Dataset

The Cars and Motorbikes dataset consists of 30 pairs of images of cars and 20 pairs of images of motorbikes from the Pascal 2007 dataset [Everingham et al., 2009]. Each pair contains a number of ground-truth correspondences and several outliers. In the dataset, the scale and viewing angle have only small variations. To investigate the performance of matching methods under large scale and viewing angle variations, for each image pair, we add a rigid transformation to the keypoints in the second image with scale parameter 4 and rotation angle π/4. We set the parameter τij,yiyj = min(dij, dyiyj ) as in previous work [Zhou and De la Torre, 2012; Zhang et al., 2016]. We randomly added 1-20 outliers from the background to the matching problem, with the result shown in Figure 2. From the figure we can see that our methods provide better average accuracy in all cases.

6 Conclusion

In this paper, we formulate hyper graph matching problems as constrained MAP inference problems in graphical models. The proposed method explicitly models only one global correspondence vector, but several local correspondence vectors, allowing us to decompose difficult high order optimization problems into a linear bipartite matching problem and several belief propagation sub-problems with small scale. We further propose a dynamic programming procedure to efficiently solve the belief propagation sub-problems. Experiments demonstrate the proposed method outperforms stateof-the-art methods for hyper graph matching.

Acknowledgments

This work was supported by NSFC Project 61301193, 61231016 and 61671385, China 863 Project 2015AA016402, ARC Project DP140102270 and DP160100703.