1

An Efficient Multilinear Optimization Framework
for Hypergraph Matching

arXiv:1511.02667v2 [cs.CV] 24 May 2016

Quynh Nguyen, Francesco Tudisco, Antoine Gautier, and Matthias Hein
Abstract—Hypergraph matching has recently become a popular approach for solving correspondence problems in computer vision as
it allows the use of higher-order geometric information. Hypergraph matching can be formulated as a third-order optimization problem
subject to assignment constraints which turns out to be NP-hard. In recent work, we have proposed an algorithm for hypergraph
matching which first lifts the third-order problem to a fourth-order problem and then solves the fourth-order problem via optimization of
the corresponding multilinear form. This leads to a tensor block coordinate ascent scheme which has the guarantee of providing
monotonic ascent in the original matching score function and leads to state-of-the-art performance both in terms of achieved matching
score and accuracy. In this paper we show that the lifting step to a fourth-order problem can be avoided yielding a third-order scheme
with the same guarantees and performance but being two times faster. Moreover, we introduce a homotopy type method which further
improves the performance.
Index Terms—Hypergraph Matching, Tensor, Multilinear Form, Block Coordinate Ascent

F

1

I NTRODUCTION

G

RAPH matching is among the most challenging tasks
of graph processing and lies at the heart of various
fields in pattern recognition, machine learning and computer vision. In computer vision, it has been used for solving
several types of problems, for instance, object recognition
[32], feature correspondence [11], [42], texture regularity discovery [20], shape matching [16], [38], [46], object tracking
[2] and surface registration [48]. Graph matching also finds
its applications in document processing tasks like optical
character recognition [17], [26], image analysis (2D and 3D)
[6], [33], [37], [43], or bioinformatics [39], [41], [44].
In this paper, we focus on the application of graph
matching to the feature correspondence problem. Given sets
of points/features extracted from two images, the task is
to find for each point/feature in the first image a corresponding point/feature in the second image while maximizing as much as possible the similarity between them. The
most simple first-order approach is to match each feature
in the first set to its nearest neighbor in the second set
by comparing their local descriptors. This naive approach
fails in the presence of repeated patterns, textures or nondiscriminative local appearances. To overcome this problem,
second-order methods [9], [14], [24], [27], [28], [29], [42],
[49], [50] try to enforce geometric consistency between pairs
of correspondences, for instance using distances between
points. However, pairwise relations are often not enough
to capture the entire geometrical structure of point set.
Therefore, higher-order approaches [8], [15], [25], [47], [48]
that take into account higher-order relations between features/points have been proposed. This paper falls into this
line of research. In particular, we cast the correspondence
problem as a hypergraph matching problem using higher-

•

Q. Nguyen, F. Tudisco, A. Gautier and M. Hein are with the Department
of Mathematics and Computer Science, Saarland University, Germany.
Corresponding Author: Quynh Nguyen.
E-mail: quynh@cs.uni-saarland.de

order similarities instead of unary or pairwise ones used
by previous methods. Depending on the chosen similarities,
this formulation allows for both scaling and rotation invariance. However, instead of concentrating on how to build
these similarities, see [7], [15], [19], [35], the focus of this
paper is how to solve the resulting optimization problem.
Compared to graph matching, algorithms for hypergraph matching are less studied in the computer vision
literature. The reason might be the difficult optimization
problem which requires not only to deal with the combinatorial constraints but also to maximize a high degree polynomial function. Most prior work has relaxed the constraints
in order to use concepts from continuous optimization [8],
[15], [25], [47]. Our main idea is to reformulate the problem. Instead of maximizing a third-degree polynomial we
maximize the corresponding multilinear form. This allows
the direct integration of assignment constraints and yields a
simple yet efficient block-coordinate ascent framework. The
idea of this paper is based on our previous work [36], where
the third-order problem has been lifted to a fourth-order
problem in order to show the equivalence of maximizing
the score function and its associated multilinear form. In
this paper we show by a careful analysis that the lifting to
a fourth-order problem can be avoided. The resulting thirdorder algorithms based on maximization of the third-order
multilinear form have the same guarantees and properties
as in [36] and achieve similar state-of-the-art matching score
and matching accuracy while being two times faster. Moreover, we provide in this paper a quite generic framework,
which allows the adaptation of our whole approach to
matching problems with different assignment constraints, as
well as to other general third-order optimization problems
with combinatorial constraints.
Extensive experiments on both synthetic and realistic
datasets show that all our algorithms, including the thirdorder and fourth-order ones, significantly outperform the
current state-of-the-art in terms of both matching score

2

and matching accuracy, in particular for very challenging
settings where different kinds of noise are present in the
data. In terms of running time, our algorithms are also on
average significantly faster than previous approaches. All
proofs and additional experimental results can be found in
the appendix.
1.1

Related Work

The graph resp. hypergraph matching problem is known
to be NP-hard except for special cases where polynomialtime algorithms exist e.g. for planar graphs [21]. In order
to make the problem computationally tractable, a myriad of
approximate algorithms have been proposed over the last
three decades aiming at an acceptable trade-off between the
complexity of the algorithm and matching accuracy. They
can be categorized from different perspectives and we refer
to [13], [18] for an extensive review. In this section, we
review only those approximate algorithms that use Lawler’s
formulations for both graph and hypergraph matching as
they are closely related to our work. In particular, the
graph matching problem can be formulated as a quadratic
assignment problem (QAP)

max xT Ax
x∈M

while the hypergraph matching problem is formulated as a
higher-order assignment problem (HAP)

max
x∈M

n
X

Fijk xi xj xk

i,j,k=1

In these formulations, A and F refer to the affinity matrix
and affinity tensor respectively, and M is some matching
constraint set depending on specific applications. Also,
depending on whether the problem is graph matching or
hypergraph matching, the corresponding algorithms will be
called second-order methods or higher-order methods.
Among recent second-order methods, Leordeanu and
Hebert [27] introduce the Spectral Matching (SM) algorithm,
and Cour et al. [14] introduce Spectral Matching with Affine
Constraint (SMAC). Both algorithms are based on the best
rank-1 approximation of the affinity matrix. Torresani et
al. [42] design a complex objective function which can be
efficiently optimized by dual decomposition. Leordeanu et
al. [29] propose the Integer Projected Fixed Point (IPFP)
algorithm, which optimizes the quadratic objective function via a gradient-type algorithm interleaved with projection onto the assignment constraints using the Hungarian
method. Lee et al. [24] tackle the graph matching problem
using stochastic sampling, whereas Cho et al. [10] introduce
a reweighted random walk (RRWM) where the reweighting is done in order to enforce the matching constraints.
Zaslavskiy et al. [46] propose a path-following algorithm
based on convex-concave relaxations. Zhou and Torre [49]
factorize the affinity matrix and then compute a convexconcave relaxation of the objective function which is finally
solved by a path-following method. Along this line there is
also work of Liu et al. [31] and Zhou and Torre [50] extending previous methods [46], [49] to deal with directed graph
matching. Recently, Cho et al. [12] propose a novel max
pooling matching algorithm, in which they tweak the power

method to better cope with noise in the affinities. Although
the algorithm comes without theoretical guarantees, it turns
out to perform very well in practice, in particular, when one
has a large number of outliers.
The hypergraph matching problem (i.e. HAP) is much
less studied in the literature. Duchenne et al. [15] formulate
the hypergraph matching problem as a tensor eigenvalue
problem and propose a higher-order power method for solving it. Zass and Shashua [47] introduce a probabilistic view
on the problem with their Hypergraph Matching Method
(HGM). Their idea is to marginalize the tensor to a vector
and solve then the lower dimensional problem. Chertok and
Keller [8] extend this idea and marginalize the tensor to a
matrix, leading to a quadratic assignment problem which is
then tackled using spectral methods. Since both methods are
based on tensor marginalization, some part of the higherorder information is lost. Moreover, they cannot handle the
one-to-one matching constraint during the iterations since
it is only considered at the final discretization step. Lee et
al. [25] extend the reweighted random walk approach of [24]
to hypergraph matching. Their algorithm aims at enforcing
the matching constraint via a bi-stochastic normalization
scheme done at each iteration. In [48], Zeng et al. propose to
use pseudo-boolean optimization [4] for 3D surface matching, where higher-order terms are decomposed into secondorder terms and the quadratic pseudo-boolean optimization
(QPBO) algorithm [22] is employed to solve the problem.
We have proposed in [36] a tensor block-coordinate ascent
framework, see also [45], for hypergraph matching based
on a multilinear reformulation of the HAP where the thirdorder problem has been lifted to a fourth-order problem and
one directly optimizes over the set of assignment matrices.
1.2

Notation

Vectors are denoted by bold lowercase letters (e.g. x). The
elements of vectors are denoted by subscripts while superscripts denote iteration numbers. Tensors will be denoted by
calligraphic uppercase letters (e.g. F ), while their associated
multilinear forms will be denoted by the corresponding
standard uppercase letter (e.g. F). The symbol ⊗ is used to
denote the tensor product.

2

H YPERGRAPH M ATCHING F ORMULATION

Hypergraphs are powerful in matching problems because
they allow modeling of relations involving groups of more
than two vertices so that higher-order information can be
integrated. We concentrate our study on 3-uniform hypergraphs, i.e. each hyperedge describes a relation between
up to 3 vertices. Every 3-uniform hypergraph can be represented by a third-order tensor. Let V and V 0 be two
point sets with n1 = |V | ≤ n2 = |V 0 |. The matching
problem consists of finding a binary assignment matrix
X ∈ {0, 1}n1 ×n2 such that Xij = 1 if vi ∈ V matches
vj0 ∈ V 0 , and Xij = 0 otherwise. The set of all possible
assignment matrices that assign each vertex of V to exactly
one vertex of V 0 is given by
n1
X

M = X ∈ {0, 1}n1 ×n2
Xij ≤ 1,
i=1

n2
X
j=1

Xij = 1 .

3

k1

i2

k2

i1
j2

j1

Fig. 1: Illustration of constructing the affinity tensor. Each entry is
computed by comparing the angles of two triangles formed by three
candidate correspondences (i1 , i2 ), (j1 , j2 ) and (k1 , k2 ).

Suppose that we have a function F : (V × V 0 )3 →
R+ which maps each pair of triples {vi1 , vi2 , vi3 } ⊂
V and {vj0 1 , vj0 2 , vj0 3 } ⊂ V 0 to its similarity weight
F(i1 ,j1 ),(i2 ,j2 ),(i3 ,j3 ) . For example, F is calculated as the
similarity of the corresponding angles of the two triangles
4i1 j1 k1 and 4i2 j2 k2 as shown in Fig. 1. Then, the score
S(X) of a matching X ∈ M is defined as follows [15]

S(X) =

n1
X

n2
X

F(i1 ,j1 ),(i2 ,j2 ),(i3 ,j3 ) Xi1 j1 Xi2 j2 Xi3 j3 .

It turns out that the score function in Eq. (1) is a special
case of the multilinear form when all the arguments are the
same, i.e. S(x) = F(x, x, x). The multilinear form is called
symmetric if it is invariant under any permutation of its
arguments, i.e. F(x1 , x2 , x3 ) = F(xσ(1) , xσ(2) , xσ(3) ) for all
x1 , x2 , x3 ∈ Rn and σ a permutation of {1, 2, 3}. Note that
if the tensor F is symmetric then F is also symmetric. In
the following, we writeP
F( · , y, z) to denote a vector in Rn
n
such that F( · , y, z)i = j,k=1 Fijk yj zk for all 1 ≤ i ≤ n,
and write F( · , P
· , z) to denote a matrix in Rn×n such that
n
F( · , · , z)ij =
k=1 Fijk zk for all 1 ≤ i, j ≤ n. Note
that the positions of the dots do not matter in the case
of a symmetric multilinear form because the function is
invariant under any permutation of its arguments.
The hypergraph matching problem can be written as the
maximization of a third-order score function S subject to the
assignment constraints,

max S(x) = max F(x, x, x) = max
x∈M

x∈M

x∈M

i1 ,i2 ,i3 =1 j1 ,j2 ,j3 =1

For ease of notation, we introduce a linear ordering in V ×V 0
so that each X ∈ M can be rewritten as a vector x ∈ {0, 1}n
with n = n1 n2 and F as a tensor in Rn×n×n . With this
convention, the score function S : Rn → R can be computed
in the following way

S(x) =

n
X

Fijk xi xj xk .

(1)

i,j,k=1

We recall that a third-order tensor T ∈ Rn×n×n is called
symmetric if its entries Tijk are invariant under any permutation of {i, j, k}. In particular, the tensor Tijk := xi xj xk is
symmetric and thus the non-symmetric part of F is “averaged” out in the computation of S(x). Therefore, without
loss of generality, we can assume that F is symmetric1 . In
this paper, we consider only terms of order 3, i.e. Fijk = 0 if
i = j , i = k or j = k . However, note that it is possible
to integrate first-order terms on the main diagonal Fiii ,
and pairwise potentials Fijj for i 6= j . This is also the
major advantage of hypergraph matching over conventional
graph matching methods because one can combine local
appearance similarities, pair-wise similarities (e.g. distances)
and higher-order invariants (e.g. angles of triplets of points)
in a unified framework.

3

M ATHEMATICAL F OUNDATIONS FOR M ULTILIN EAR O PTIMIZATION F RAMEWORK

In this section, we derive the basis for our tensor block coordinate ascent framework for solving hypergraph matching
problems. The general idea is to optimize instead of the
score function S the associated third-order multilinear form
F : Rn × Rn × Rn → R which is defined as
F(x, y, z) =

n
X

Fijk xi yj zk .

(2)

i,j,k=1

1. Note that F can always be symmetrized
(without changing the
1 P
score function) as follows F̃ijk = 3!
σ∈S3 Fσ(i)σ(j)σ(k) , where
i, j, k = 1, . . . , n and S3 is the permutation group of three elements.

n
X

Fijk xi xj xk . (3)

i,j,k=1

In [36] we propose to solve instead the equivalent problem
defined as

max S 4 (x) = max F4 (x, x, x, x):=
x∈M

x∈M

n
X

4
Fijkl
xi xj xk xl (4)

i,j,k,l=1

where the new symmetric fourth-order tensor F 4 is obtained by lifting the third order tensor F via
4
Fijkl
= Fijk + Fijl + Fikl + Fjkl .

(5)

The reason for this lifting step is discussed below. The
main idea of [36] is then to optimize instead of the score
function the associated multilinear form. The following theorem establishes equivalence of these problems under the
assumption that the score function S 4 is convex on Rn .
Theorem 3.1 (Nguyen et al. [36]). Let F 4 be a symmetric
fourth-order tensor and suppose the corresponding score
function S 4 : Rn → R is convex. Then it holds for any
compact constraint set D ⊂ Rn ,

max F4 (x, x, x, x) = max F4 (x, x, y, y)
x∈D

x,y∈D

=

max

x,y,z,t∈D

F4 (x, y, z, t).

Note that a given score function S 4 need not be convex.
In [36] we show that every fourth-order score function can
be modified so that it becomes convex, while solving (4)
for the original and the modified score function turns out
to be equivalent. The maximization of the multilinear form
F4 (x, y, z, t) can then be done in two ways. Either one fixes
all but one or all but two arguments and maximizes the
remaining arguments over the set of assignment matrices
M . As shown in [36], this scheme boils down to solving
a sequence of Linear Assignment Problems (LAPs) or a sequence of (approximate) QAPs which can be done efficiently
by existing lower-order methods.
The convexity of the score function is crucial for the
proof of Theorem 3.1. The lifting step (5) has been proposed

4

in [36] as non-trivial third-order score functions are not
convex on Rn as the following lemma shows.
Lemma 3.2 (Nguyen et al. [36]). Let S be a third-order score
function defined as in Eq. (1). If S is not constant zero,
then S : Rn → R is not convex.
In the following we show that a statement similar to Theorem 3.1 holds also in the third-order case. This allows us to
propose algorithms directly for the third-order case without
the lifting step, leading to two times faster algorithms. The
new key insight to establish this result is that convexity
on whole Rn is not required. If one just requires that the
third-order score function is convex on the positive orthant,
then Lemma 3.2 is no longer true. In fact convexity is only
required in a “pointwise” fashion which we make precise in
Theorem 3.3. Similar to our previous Theorem 3.1, we state
the main Theorem 3.3 of this paper for general compact sets
D ⊂ Rn which includes as a special case the set of assignment matrices. Although this general form of the theorem
may allow generalizations of our proposed framework to
other applications, in this paper we only focus on the case
D = M for the application of hypergraph matching.
Theorem 3.3. Let D ⊂ Rn be a compact set, F : Rn × Rn ×
Rn → R a symmetric third-order multilinear form and
S : Rn → R the associated third-order score function
such that S(x) = F (x, x, x), then we have the following
chain of implications: 1) ⇒ 2) ⇒ 3) ⇒ 4).
1) ∇2 S(x) is positive semidefinite for all x ∈ D.
2) F(x, y − z, y − z) ≥ 0 for all x, y, z ∈ D.
3) It holds for all x, y, z ∈ D
F(x, y, z) ≤

max
u∈{x,y,z}

F(u, u, u).

(6)

4) The optimization of the multilinear form is equivalent
to the optimization of its associated score function

max F(x, y, z) = max F(x, x, y) = max F(x, x, x).

x,y,z∈D

x,y∈D

x∈D

As a twice continuously differentiable function is convex
if and only if its Hessian is positive semi-definite on its
domain, the first requirement is much weaker than requiring
convexity on whole Rn . Thus Theorem 3.3 allows us to
establish for D = M inequality (6) for a nontrivial class of
third-order score functions, which is the key result necessary
to establish the equivalence of the maximization of the score
function and the maximization of the multilinear form. This
is the key requirement to show monotonic ascent for our
block-coordinate optimization scheme. In Proposition 3.4 we
show how to modify an existing third-order score function
so that it satisfies condition 1) for D = M . Again this
modification turns out to be constant on M and thus leads to
equivalent optimization problems. This altogether leads to
two new algorithms for the hypergraph matching problem
(1) similar to our previous work [36]. The key idea of both
algorithms is to use block-coordinate ascent updates, where
all but one argument or all but two arguments of the multilinear form are fixed and one optimizes over the remaining
ones in either case. The inequality (6) then allows us to connect the solutions of the two optimization problems. In both
variants, we directly optimize over the discrete set M of
possible matchings, that is, there is no relaxation involved.

Moreover, we theoretically prove monotonic ascent for both
methods. In all our experiments, the proposed third-order
methods achieve competitive matching score and matching
accuracy to the fourth-order ones [36] while outperforming other state-of-the-art approaches. On the other hand,
we achieve a speed up of factor 2 over the fourth-order
algorithms due to the ability of working directly with the
original tensor.
Compared to [36] we avoid the lifting step while maintaining the same theoretical guarantees. Moreover, we use
the weaker conditions on the score function in Theorem
3.3 to propose two variants of our algorithms based on
homotopy methods as presented in Section 3.3.
In general, the multilinear form F might not fulfill the
conditions of Theorem 3.3. Therefore, we propose to modify
F by adding to it some multilinear form G so that the first
statement of Theorem 3.3 is satisfied for the new function. At
the same time the modification does not affect the solution
of our problem since the added term is constant on the set
of assignment matrices. This will be detailed later on in
Proposition 3.4. While finding the best possible modification
is difficult, one can consider some canonical choices. As the
first choice, one can implicitly modify F by considering a
2
new score function defined as Sα (x) = S(x) + α kxk2 , and
then show that there exists some α0 ∈ R such that ∇2 Sα
is positive semidefinite everywhere in M for all α ≥ α0 .
However, the modified score function Sα then becomes inhomogeneous and thus there does not exist any multilinear
form Fα such that Sα (x) = Fα (x, x, x).
A homogeneous modification would be Sα (x) = S(x) +
3
3
α kxk3 . Indeed, onePcan associate to kxk3 the multilinear
n
form G(x, y, z) =
i=1 hei , xi hei , yi hei , zi , where ei is
the i-th standard basis vector of Rn . Moreover, it holds that
3
G(x, x, x) = kxk3 . In this paper, we consider a family of
modifications which includes the one above by perturbing
the set of basis vectors as ēi = ε1+(1−ε)ei with 1 being the
3
all-ones vector. Note that kxk3 is included in this family as a
special case for ε = 0. The main idea is that one can later on
optimize ε to find the best lower bound on α. That is to say,
we aim at a minimal amount of S such that the conditions
from Theorem 3.3 are satisfied. As shown in the proof of
Proposition 3.4, ε = 1/3 is the optimal value for which the
best bound on α can be obtained. Thus, we use this value
in the rest of the paper. In particular, let ēi = 31 1 + 23 ei we
define the symmetric third-order tensor G as

G=

n
X

ēi ⊗ ēi ⊗ ēi ,

(7)

i=1

and the associated multilinear form as
n
X
G(x, y, z) =
hēi , xi hēi , yi hēi , zi .

(8)

i=1

Proposition 3.4. Let F : Rn × Rn × Rn → R be a symmetric
third-order multilinear form and G : Rn × Rn × Rn → R
as defined in Eq. (8). We consider the new multilinear
form Fα defined by
Fα (x, y, z) = F(x, y, z) + α G(x, y, z).
qP
n
2
For all α ≥ 27
max
j,k=1 (Fijk ) it holds
4
i=1,2,...,n

5

1) Fα is a symmetric third-order multilinear form.
2) ∇2 Sα (x) is positive semidefinite for all x ∈ M.
3) The new problem is equivalent to the original one

arg max Sα (x) = arg max S(x).
x∈M

x∈M

Discussion. In [34] a general convexification strategy for
arbitrary score functions has been proposed where, similar
to our modification, the added term is constant on the set
of assignment matrices M . However, as the added term
is in-homogeneous, it cannot be extended to a symmetric
multilinear form and thus does not fit to our framework.
In second-order graph matching also several methods use
convexified score functions in various ways [46], [49], [50].
However, for none of these methods it is obvious how to extend it to a third-order approach for hypergraph matching.
Once Fα is symmetric and ∇2 Sα (x) is positive semidefinite at every x in M , it follows from Theorem 3.3 that

max Fα (x, y, z) = max Sα (x) ≡ max S(x).

x,y,z∈M

x∈M

x∈M

Therefore, block-coordinate ascent schemes that optimize
Fα over assignment constraints can be derived in a similar
fashion to [36]. In particular, we propose below two variants,
one for solving
max Fα (x, y, z)

Algorithm 1 BCAGM3
1: Input: α

=

27
max
4 i=1,2,...,n

qP

n
0
0 0
2
j,k=1 (Fijk ) , (x , y , z )

∈ Rn ,

k = 0, m = 0
2: Output: x∗ ∈ M
3: loop
4:
x̃k+1 = arg max Fα (x, yk , zk )
x∈M

5:

ỹk+1 = arg max Fα (x̃k+1 , y, zk )

6:

z̃k+1 = arg max Fα (x̃k+1 , ỹk+1 , z)

7:
8:

if Fα (x̃k+1 , ỹk+1 , z̃k+1 ) = Fα (xk , yk , zk ) then
um+1 =
arg max
Fα (u, u, u)

y∈M

z∈M

u∈{x̃k+1 ,ỹk+1 ,z̃k+1 }

9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:

if Fα (x̃k+1 , ỹk+1 , z̃k+1 ) < Fα (um+1 , um+1 , um+1 )
then
xk+1 = yk+1 = zk+1 = um+1
m=m+1
else
return um+1
end if
else
xk+1 = x̃k+1 , yk+1 = ỹk+1 , zk+1 = z̃k+1
end if
k =k+1
end loop

x,y,z∈M

which leads to Algorithm 1 and the other for solving

max Fα (x, y, y)

x,y∈M

which leads to Algorithm 2. Both variants come along with
theoretical guarantees including strict monotonic ascent and
finite convergence.
3.1 Tensor Block Coordinate Ascent via a Sequence of
Linear Assignment Problems
The first algorithm uses a block coordinate ascent scheme
to optimize the multilinear function Fα (x, y, z) where two
arguments are fixed at each iteration and the function is
maximized over the remaining one. This boils down to
solving a sequence of LAPs, which can be solved to global
optimality by the Hungarian method [5], [23]. As we are
optimizing the multilinear form, we need to guarantee that
our algorithm produces finally a homogeneous solution,
i.e. x = y = z. Moreover, we require that our algorithms
achieve ascent not only in the multilinear form but also in
the original score function. These are the two main reasons
why we use the inequality (6) at step 7) of Algorithm 1
to get ascent in the score function. The following theorem
summarizes the properties of Algorithm 1.
Theorem 3.5. Let F be
a symmetric third-order tensor and
qP
n
2
α ≥ 27
max
j,k=1 (Fijk ) . Then the following
4
i=1,2,...,n

holds for Algorithm 1:
1) The sequence Fα (xk , yk , zk ) for k = 1, 2, . . . is strictly
monotonically increasing or terminates.
2) The sequence of scores S(um ) for m = 1, 2, . . . is
strictly monotonically increasing or terminates. For every m, um ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations.

We would like to note that all statements of Theorem
3.5 remain valid for α = 0 if one of the statements 1)-3) in
Theorem 3.3 holds for D = M . In practice, this condition
might be already satisfied for some constructed affinity
tensors. Thus we adopt the strategy in [36] by first running
Algorithm 1 with α = 0 until we q
get no further ascent and
Pn
27
2
only then we set α = 4 max
j,k=1 (Fijk ) . It turns
i=1,2,...,n

out that in our experiments often the first phase with α = 0
leads automatically to a homogeneous solution. In this case,
the algorithm can be stopped as the following lemma shows.
Lemma 3.6. Suppose Algorithm 1 runs with α = α0 for some
α0 ∈ R+ . If x̃k+1 = ỹk+1 = z̃k+1 =: u holds at some
iteration k then for all α ≥ α0 , it holds

u = arg max Fα (z, u, u).
z∈M

Lemma 3.6 shows that we can stop Algorithm 1 whenever
the first phase leads already to a homogeneous solution as
no further ascent is possible. In particular, there is no need to
go for a new phase with larger value of α since the current
iterate (u, u, u) is already optimal to all local updates of
Algorithm 1, i.e. steps (4)-(6), for all α ≥ α0 = 0. Also, the
fact that the iterate is already homogeneous implies that no
further possible ascent can be achieved at step (7).
3.2 Tensor Block Coordinate Ascent via Alternating between Quadratic and Linear Assignment Problems
The second algorithm uses a block coordinate ascent scheme
to optimize the multilinear form Fα where now one resp.
two arguments are fixed and we optimize over the remaining ones in an alternating manner. The resulting scheme
alternates between QAPs and LAPs. While the LAP can be

6

efficiently solved using the Hungarian algorithm as before,
the QAP is NP-Hard. Thus a globally optimal solution is
out of reach. However, our algorithm does not require the
globally optimal solution at each step in order to maintain
the theoretical properties. It is sufficient that the sub-routine
method, called Ψ yields monotonic ascent w.r.t. the current
variable, that is,
E
D
(9)
z = Ψ(A, yk ) s.t. hz, Azi ≥ yk , Ayk
where yk is the current iterate and the nonnegative symmetric A ∈ Rn×n is Fα (x̃k+1 , · , · ) in our algorithm. As
in Algorithm 1, we go back to the optimization of the
score function in step 6) using inequality (6). The following
theorem summarizes the properties of Algorithm 2.
Theorem 3.7. Let F q
be a symmetric third-order tensor and
Pn
2
α ≥ 27
max
j,k=1 (Fijk ) . Let Ψ be an algorithm
4
i=1,2,...,n

for the QAP which yields monotonic ascent w.r.t. the
current iterate. The following holds for Algorithm 2:
1) The sequence Fα (xk , yk , yk ) for k = 1, 2, . . . is strictly
monotonically increasing or terminates.
2) The sequence of scores S(um ) for m = 1, 2, . . . is
strictly monotonically increasing or terminates. For every m, um ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations.
Algorithm 2 BCAGM3+Ψ
1: Input:

2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:

α=

27
max
4 i=1,2,...,n

qP

n
0
0
2
j,k=1 (Fijk ) , (x , y )

∈ Rn ,

k = 0, m = 0,
z = Ψ(A, yk ) is a sub-routine method that provides
monotonic ascent for the QAP, i.e. hz, Azi ≥ yk , Ayk
Output: x∗ ∈ M
loop
x̃k+1 = arg max Fα (x, yk , yk )
x∈M

ỹk+1 = Ψ Fα (x̃k+1 , · , · ), yk
if Fα (x̃k+1 , ỹk+1 , ỹk+1 ) = Fα (xk , yk , yk ) then
um+1 = arg max Fα (u, u, u)
u∈{x̃k+1 ,ỹk+1 }
k+1
k+1
k+1

if Fα (x̃
, ỹ
, ỹ
) < Fα (um+1 , um+1 , um+1 )
then
xk+1 = yk+1 = um+1
m=m+1
else
return um+1
end if
else
xk+1 = x̃k+1 , yk+1 = ỹk+1
end if
k =k+1
end loop

In analogy to Theorem 3.5, all statements of Theorem
3.7 remain valid for α = 0 if one of the statements 1)-3)
in Theorem 3.3 holds for D = M .. Thus we use the same
initialization strategy with α = 0 as described above. There
are several methods available which we could use for the
sub-routine Ψ in the algorithm. We decided to use the recent
max pooling algorithm [12] and the IPFP algorithm [29], and

then use the Hungarian algorithm to turn their output into
a valid assignment matrix in M at each step. It turns out
that the combination of our tensor block coordinate ascent
scheme using their algorithms as a sub-routine yields very
good performance on all datasets. In case of Algorithm 2 a
statement as in Lemma 3.6 is not possible, as the subroutine
Ψ usually can only deliver a local solution to the subproblem and Lemma 3.6 relies on the fact that the subproblem
can be solved globally optimally. However, we observe in
our experiments that almost always Algorithm 2 does not
achieve further ascent when its iterate is homogeneous, thus
we recommend to stop Algorithm 2 in this case.
3.3 A Homotopy Tensor Block Coordinate Ascent
Scheme
Both algorithms consist of two phases – the first phase
uses α = 0 as q
initialization and the second phase uses
Pn
2
α = 27
max
j,k=1 (Fijk ) . However, if the value of
4
i=1,2,...,n

α is too large for the second phase, often no further improvement is achieved. This phenomenon can be explained by the
fact that a large modification term tends to homogenize the
variables more quickly, i.e. xk = yk = zk , which makes the
algorithm get stuck faster at a critical point.
To fix this, one first observes that the inequality (6) lies
at the heart of both algorithms. On one hand, it serves as a
sufficient condition for the equivalence of maximizing score
functions and multilinear form. On the other hand, it might
help the algorithms to jump to better solutions in case they
reach a stationary state. Both methods guarantee this inequality to be satisfied for all the tuples {x, y, z} ∈ M at the
same time by using a rather large value of α. However, this
might be too conservative as the algorithm itself converges
typically rather quickly and thus visits only a small number
of feasible tuples. Thus, we propose below to satisfy the
inequality (6) as the algorithm proceeds. This is done by
updating α accordingly during the algorithm which yields
a homotopy method.
Proposition 3.8. Let F : Rn × Rn × Rn → R be a symmetric
third-order multilinear form and G : Rn × Rn × Rn → R
as defined in Eq. (8). For each α ∈ R, let Fα be defined as
Fα (x, y, z) = F(x, y, z) + α G(x, y, z).
Then the following holds
1) Fα is a symmetric multilinear form.
2) For all non-homogeneous tuple (x, y, z) ∈ M (i.e. x 6=
y or y 6= z or z 6= x), the following inequality holds
Fα (x, y, z) ≤

max
u∈{x,y,z}

Fα (u, u, u)

(10)

if and only if
F(x, y, z) −

α ≥ Λ(x, y, z) :=

max
u∈{x,y,z}

F(u, u, u)

G(x, x, x) − G(x, y, z)

.

3) For every α ∈ R, it holds

arg max Sα (x) = arg max S(x).
x∈M

x∈M

Proposition 3.8 suggests that if α is chosen such that

α≥

max

x,y,z∈M
x6=y or y6=z or z6=x

Λ(x, y, z),

(11)

7

then the inequality (6) is satisfied for all x, y, z ∈ M , for
which, it follows from Theorem 3.3

max Fα (x, y, z) = max Sα (x) ≡ max S(x).

x,y,z∈M

x∈M

x∈M

Thus, Algorithm 1 or 2 can be applied again to solve the
problem. However, the computation of the optimal bound
given in (11) is not feasible as the number of feasible tuples
grows exponentially with the size of the problem. Therefore,
we adaptively update α as follows. First, we use block
coordinate ascent steps as in the previous methods with
α = 0. When such a step achieves no further ascent and
inequality (6) is violated, we slightly increase α so that
the inequality becomes satisfied for the current iterate, and
the algorithm is continued. The whole scheme is shown in
Algorithm 3. Note that a small ξ > 0 is used in step 13) to
prove the convergence, which is analyzed below.
Algorithm 3 Adapt-BCAGM3
1: Input: α = 0, (x0 , y0 , z0 ) ∈ Rn , k
2: Output: x∗ ∈ M
3: loop
4:
x̃k+1 = arg max Fα (x, yk , zk )

= 0, m = 0, ξ > 0

x∈M

5:

ỹk+1 = arg max Fα (x̃k+1 , y, zk )

6:

z̃k+1 = arg max Fα (x̃k+1 , ỹk+1 , z)

7:
8:

if Fα (x̃k+1 , ỹk+1 , z̃k+1 ) = Fα (xk , yk , zk ) then
um+1 =
arg max
Fα (u, u, u)

9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:

y∈M

z∈M

u∈{x̃k+1 ,ỹk+1 ,z̃k+1 }
k+1
k+1 k+1

if Fα (x̃
, ỹ
, z̃
) < Fα (um+1 , um+1 , um+1 )
then
xk+1 = yk+1 = zk+1 = um+1
m=m+1
else if α does not satisfy Eq. (11) then
α = Λ(x̃k+1 , ỹk+1 , z̃k+1 ) + ξ
xk+1 = x̃k+1 , yk+1 = ỹk+1 , zk+1 = z̃k+1
else
return um+1
end if
else
xk+1 = x̃k+1 , yk+1 = ỹk+1 , zk+1 = z̃k+1
end if
k =k+1
end loop

Note that α is always strictly increasing whenever Algorithm 3 enters step 12). This can be seen as follows.
If the current value of α is greater than or equal to
Λ(x̃k+1 , ỹk+1 , z̃k+1 ), then inequality (10) holds for the current tuple {x̃k+1 , ỹk+1 , z̃k+1 }, which implies the algorithm
should not have entered step 12), leading to a contradiction.
For any fixed value of α, Algorithm 3 works exactly
in the same fashion as Algorithm 1, thus, it also yields
strict monotonic ascent. The only difference is that the first
algorithm updates α only with a potentially large value,
whereas the new one splits this update into multiple phases.
In case α exceeds the bound given in (11) at some
iteration, the “strict” inequality will hold for all x, y, z ∈ M ,
in which case α will no longer be updated. This, combined
with the fact that the algorithm yields strict monotonic

ascent for fixed α and Fα is bounded from above, guarantees
the whole scheme converges in a finite number of iterations.
It should be emphasized that for any value of α, the
optimization of Sα (x) and S(x) on M is always equivalent,
but the equivalence of the optimization of the multilinear
form Fα (x, y, z) and that of the score function Sα (x) is only
guaranteed when α satisfies Eq. (11).
Finally, we adopt the same homotopy approach as described above for Algorithm 2 to obtain Algorithm 4. All
the properties of Algorithm 3 transfer to Algorithm 4. In
particular, the algorithm achieves strict monotonic ascent
within each phase and has finite termination.
Algorithm 4 Adapt-BCAGM3+Ψ
1: Input: α = 0, (x0 , y0 ) ∈ Rn , k =
2: Output: x∗ ∈ M
3: loop
4:
x̃k+1 = arg max Fα (x, yk , yk )
5:
6:
7:

0, m = 0, ξ > 0

x∈M

ỹk+1 = Ψ Fα (x̃k+1 , · , · ), yk
if Fα (x̃k+1 , ỹk+1 , ỹk+1 ) = Fα (xk , yk , yk ) then
um+1 = arg max Fα (u, u, u)
u∈{x̃k+1 ,ỹk+1 }

9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:

if Fα (x̃k+1 , ỹk+1 , ỹk+1 ) < Fα (um+1 , um+1 , um+1 )
then
xk+1 = yk+1 = um+1
m=m+1
else if α does not satisfy Eq. (11) then
α = Λ(x̃k+1 , ỹk+1 , ỹk+1 ) + ξ
xk+1 = x̃k+1 , yk+1 = ỹk+1
else
return um+1
end if
else
xk+1 = x̃k+1 , yk+1 = ỹk+1
end if
k =k+1
end loop

4

C OMPUTATIONAL A NALYSIS

8:

Even though all our algorithms are shown to converge in
a finite number of iterations, it is difficult to bound the
number of iterations. Therefore, we focus on analyzing the
cost of each iteration for both methods.
For Algorithm 1, each iteration requires solving three
linear assignment problems, in which the gradient of the
multilinear form needs to be computed, e.g. , the iterate
D
E
x̃k+1 = arg max Fα (x, yk , zk ) = arg max x, Fα ( · , yk , zk )
x∈M

x∈M

requires to compute the vector
Fα ( · , yk , zk )r =

n
X

Frst ysk zkt ,

∀ r = 1, 2, . . . , n.

s,t=1

Without any knowledge of yk and zk , this computation
can be done efficiently by passing through all the entries
of F once, i.e. in O(T ) where T is the number of nonzero
entries of the tensor [15]. However, this approach potentially
considers many unnecessary elements, for instance, those

8

entries indexed by (r, s, t) for which ysk = 0 or ytk = 0
do not contribute anything to the final result, thus, should
be avoided during computation. One special property of
our algorithm is that it always delivers binary assignment
matrices at each step, thus, the sparsity of the variables can
be taken into account at each step. As yk and zk are binary
assignment matrices, they have exactly n1 nonzero values.
Therefore, we first go over the positions of all nonzero components in each variable. In particular, for each pair (s, t)
such that ysk > 0 and zkt > 0, we can access all the related
entries Frst using a linked list or hashing method. With that,

the above computation can be done in O min{T, n21 K} ,
where K  n = n1 n2 is the average number of tensor
entries that are associated to each pair (s, t).
For solving the LAP, we use the Hungarian method2 .
The Hungarian method
for the LAP has cubic complexity

of each
O max{n1 , n2 }3 . Therefore, the total complexity

iteration is O max{n1 , n2 }3 + min{T, n21 K} . For K =
max{n1 , n2 } and n21 K < T , the total complexity becomes
O max{n1 , n2 }3 , which is roughly the same as the complexity of the Hungarian method.
For Algorithm 2, each iteration consists of one LAP and
one QAP. The complexity for solving LAP is the
 same as
above, that is, O max{n1 , n2 }3 + min{T, n21 K} . For the
QAP, the iterate

ỹk+1

=
=

arg max Fα (x̃k+1 , y, y)
y∈M
D
E
arg max y, Fα (x̃k+1 , · , · )y
y∈M

requires to compute the matrix
Fα (x̃k+1 , · , · )st =

n
X

Frst x̃k+1
,
r

∀ s, t = 1, 2, . . . , n.

r=1

Adopting the same technique as before, we first pass
through all the nonzero components of x̃k+1 , and at each
index r where x̃k+1
> 0, we access all the associated
r
tensor entries Frst via the use of linked list. This costs
O (min{T, n1 L}), where L  (n1 n2 )2 is the number of
nonzero tensor entries associated to each index r on average. Given the above matrix, let us assume our quadratic
solvers need t iterations to solve the problem, where each
iteration costs approximately the same as the computation of the matrix3 . Then the computation of ỹk+1 takes
O (t min{T, n1 L}). The total complexity of each iteration
is thus the sum of two complexities, one for solving LAP
and the other for solving QAP, that is, O(max{n1 , n2 }3 +
min{T, n21 K} + t min{T, n1 L}). For L = max{n1 , n2 }2
and n1 L < T , the whole complexity is simplified to
O(t max{n1 , n2 }3 ), which is equivalent to solving a sequence of LAPs using the Hungarian method. Therefore, we
observe in our experiments that often Algorithm 2 is quite
slower than Algorithm 1 since the later only needs to solve
three LAPs at each step.
Finally, the cost per iteration of Algorithm 3 and Algorithm 4 is the same as the cost per iteration of Algorithm 1
and Algorithm 2 respectively.
2. We employ the implementation of Cyrill Stachniss due to its
simplicity and efficiency: http://www2.informatik.uni-freiburg.de/
∼stachnis/misc.html
3. This assumption is reasonable for those quadratic solvers that
involve only matrix-vector multiplications at each step.

5

E XPERIMENTS

This section reports experimental results on standard benchmark datasets (one synthetic and three real) and compares
our algorithms with state-of-the-art hypergraph matching
as well as several graph matching methods.
Hypergraph Matching Methods Tensor Matching (TM)
[15], Hypergraph Matching via Reweighted Random Walks
(RRWHM) [25] and probabilistic hypergraph matching
method (HGM) [47].
Graph Matching Methods Max Pooling Matching
(MPM) [12] and Integer Projected Fixed Point (IPFP) [29],
Reweighted Random Walks for Graph Matching (RRWM)
[10] and Spectral Matching (SM) [27]. Since HGM, RRWM
and SM are outperformed by other methods [12], [25], [36],
we do not show their performance for saving space.
We denote our Algorithm 1 as BCAGM3 which uses
the Hungarian method at each iteration, and Algorithm 2
as BCAGM3+MP resp. BCAGM3+IPFP which uses MPM
[12] and IPFP [29] respectively as subroutine. MPM has
recently outperformed other graph matching algorithms in
the presence of a large number of outliers. Our homotopy
algorithms will be denoted as Adapt-BCAGM3 (Algorithm
3) and Adapt-BCAGM3+MP and Adapt-BCAGM3+IPFP
(Algorithm 4). We also compare our third-order methods
with the fourth-order counterparts from [36] that we denote
as BCAGM, BCAGM+MP and BCAGM+IPFP accordingly.
Note that we differentiate between the proposed third-order
algorithms and their fourth-order counterparts in [36] by
adding a 3 in their names. For all the algorithms, we use the
authors’ original implementation.
In the experiments below, we use the all-ones vector as a
starting point for all the algorithms as done in the previous
work [12], [15], [25] . Moreover, we use the Hungarian
method to turn the output of all algorithms into a proper
matching.
Generation of affinity tensor/matrix:
To build the affinity tensor for the third-order algorithms, we follow the approach of Duchenne et al. [15]. For
each triple of points we compute a feature vector f from
the angles of the triangle formed by those three points as
illustrated in Fig. 1. Let fi1 ,i2 ,i3 and fj1 ,j2 ,j3 denote the
feature vectors of two triples (i1 , i2 , i3 ) and (j1 , j2 , j3 ), we
compute the third-order affinity tensor as:
2

F(i1 ,j1 ),(i2 ,j2 ),(i3 ,j3 ) = exp(−γ kfi1 ,i2 ,i3 − fj1 ,j2 ,j3 k2 ) (12)
where γ is the inverse of the mean of all squared distances.
As shown by Duchenne et al. [15] this higher-order affinity is
more efficient than pairwise distances which have been used
in previous graph matching methods [14], [27], [29], [49]. In
particular, it can deal with scale transformations that is not
possible for affinities based on distances. Moreover, more
complicated higher-order affinities can also be derived to
achieve (partial) invariance with respect to other transformations such as affine and projective transformations [15].
This explains why higher-order approaches lead typically
to superior results compared to second-order approaches.
However, one should take into account the computational
cost of constructing the affinity tensor. Given two images
with n1 and n2 points respectively, then the number of
candidate correspondences is n = n1 n2 . Thus a third-order

400
300
200
100
0

0

200

400

600
500
400
300
200
100
0

600

0

BCAGM3

200

400

600

500

500

400

400
300
200
100
0

600

0

BCAGM3+IPFP

200

400

600

BCAGM3+MP

Matching Score

500

Matching Score

Adapt-BCAGM3

600

Adapt-BCAGM3+MP

Adapt-BCAGM3+IPFP

9

400
300
200
100

(a) Face

300
200
100
0

0

300
200
100
0

0

200

400

500
400
300
200
100
0

600

0

BCAGM3

200

400

100

# of outliers

150

0

200

400

200

1

1

0

0

200

400

600

0.8

0.8

BCAGM3+MP

0.6
0.4

200
100
0

0

200

400

600

BCAGM3

Adapt-BCAGM3+MP

Adapt-BCAGM3+IPFP

Adapt-BCAGM3

300

500
400
300
200
100
0

0

200

400

0.6
0.4
0.2

0.2

400

200

100

(b) Duck

500

150

n out

300

600

600

100

500

BCAGM3+IPFP

600

50

# of outliers

n out

Accuracy

400

50

600

Accuracy

Adapt-BCAGM3

500

600

Adapt-BCAGM3+MP

Adapt-BCAGM3+IPFP

0
600

600

0

0

500

0

400

50

200

0

150

0

BCAGM3+IPFP

200

400

600

0

200

50

100

# of outliers

n out

150

200

n out

(b) nin = 10, σ = 0.03, scale = 1.5

(a) nin = 10, σ = 0.01

100

600

100

# of outliers

300

Fig. 3: Matching point sets in R2 . The number of outliers is varied from
0 to 200. See Fig. 4 for the legend. (Best viewed in color.)

BCAGM3+MP

(c) Winebottle
Fig. 2: Scatter plots showing matching score of adaptive vs. nonadaptive methods over different datasets.

1

2000

2

Q
2
exp(−|dP
i1 i2 − dj1 j2 | /σs )

(13)

Q

where dP
i1 i2 and dj1 j2 are the Euclidean distances between
two points i1 , i2 ∈ P and j1 , j2 ∈ Q respectively, σs is a normalization term that is specified later in each experiment.
5.1

Advantages of adaptive methods

In this section, we compare our adaptive algorithms with
the non-adaptive versions. In particular, we conduct our
experiments on a realistic dataset [9] which consists of three
object classes: face, duck and wine bottle. Each class contains
at least 50 images with different instances.4 Every image
in the dataset is annotated with 10 ground truth points
for matching. We perform our experiments separately on
4. The dataset consists of 109 Face, 50 Duck, 66 Wine bottle and is
available from [1].

0.8
IPFP
MPM
TM
RRWHM
BCAGM
BCAGM+MP
BCAGM+IPFP
Adapt-BCAGM3
Adapt-BCAGM3+MP
Adapt-BCAGM3+IPFP

1500

Accuracy

affinity tensor consists of
entries, which is not only
infeasible to compute in reasonable time but also requires
huge memory resources even for moderate n1 and n2 .
Therefore we consider only sparse tensors as in [15], [25]. In
particular, we randomly sample n2 triples of points from the
first point set, and find for each of them k nearest neighbors
in the second point set using approximate nearest neighbor
search algorithm [3]. Each pair of tuples then determines
a non-zero entry of the tensor to be calculated while the
remaining entries are zero. In principle, a large value of k
would increase the chance that correct matches are captured
by these pairs of tuples which potentially improves the
result, however, this comes at additional searching cost.
Therefore, we empirically set k = 300 in all our experiments
to achieve reasonable running time for all the methods.
For second-order methods, we construct the affinity matrix by following previous work [10], [12], which estimates
the pairwise similarity as

Matching Score

n3 
3

1000

0.6
0.4
0.2

500

0

0

0.1

0.2

Deformation noise

0.3

0.4

σ

0

0.1

0.2

Deformation noise

0.3

0.4

σ

Fig. 4: Performance of all algorithms with varying deformation noise. The
number of inlier points is fixed as nin = 20. (Best viewed in color.)

each object class, for which we test all our algorithms on
2000 image pairs randomly selected from the same class.
In order to make the tasks more difficult, we randomly
add to one image several points that are extracted using a
SIFT detector. This simulates realistic outlier points resulting
from cluttered background in images. In particular, the
number of outliers is varied from 0 to 20 which yields 4200
test cases per object class in total. The experimental results
are shown in the scatter plot of Fig. 2 where the matching
score of the non-adaptive method (x-axis) is plotted against
the matching score of the adaptive version (y-axis). One can
clearly see that, in all the experiments, the adaptive methods
always achieve similar or better matching score than the
non-adaptive methods, which confirms the advantage of
gradually updating α.
To better understand our algorithms especially in terms
of matching accuracy, in the next sections we will compare
them with state-of-the-art methods over different datasets.
Since the proposed adaptive versions perform better than
the non-adaptive approaches, we show below only the
performance of the adaptive ones. The appendix provides
more detailed comparison between both approaches as well
as the fourth-order methods from [36].

10

5.2

Synthetic Dataset
1

Accuracy

0.8

0.8

IPFP
MPM
TM
RRWHM
BCAGM
BCAGM+MP
BCAGM+IPFP
Adapt-BCAGM3
Adapt-BCAGM3+MP
Adapt-BCAGM3+IPFP

0.6
0.4
0.2

0.6
0.4
0.2

0

0
20

40

60

80

100

20

Baseline

0.8

400

Matching Score

500

0.6
0.4
0.2
0
40

60

60

80

100

(b) 20 pts vs. 30 pts

1

20

40

Baseline

(a) 30 pts vs. 30 pts

Accuracy

In this section we are trying to match point sets in R which
is a standard benchmark for testing graph and hypergraph
matching methods [12], [14], [25], [28], [50]. We follow a
similar procedure as in [12] for the creation of two point
sets P and Q. In order to create P we sample nin inlier
points from a Gaussian distribution N (0, 1). These points
are then copied to Q and subsequently modified by several
operations in the following order: 1) the coordinates of
every point are multiplied with a scale factor; 2) additive
Gaussian noise N (0, σ 2 ) is added to every point where σ is
hereinafter referred to as deformation; 3) on top of that nout
random points, that we refer to as outliers, are added to Q
by sampling from N (0, 1). Depending on different scenarios
described below, some operations in the above chain might
be absent but their order is the same. Based on that, we
test all matching algorithms under different changes in the
data: outliers, deformation and scaling. In this section, we
consider two main settings: outlier and deformation setting.
In all plots in this section, each quantitative result was
obtained by averaging over 100 random trials. The accuracy
is computed as the ratio between the number of correct
matches and the number of inlier points. We use σs = 0.5
in the affinity construction for the second-order methods as
suggested by [12].
In the outlier setting, we perform two experiments. In the
first experiment, we vary the number of outliers from 0 to
very challenging 200, which is 20-times more than the number of inliers nin , while σ is set to 0.01 to simulate noise in
realistic settings and no scaling is used here (i.e. scale = 1.0).
In the second case, we slightly increase σ to 0.03 and scale
all the points in Q by a factor of 1.5, making the problem
more realistic and also much more difficult. This experiment shows the robustness of higher-order methods over
second-order ones w.r.t. scale transformations. The results
are shown in Fig. 3. Overall, our higher-order algorithms
show remarkable tolerance to the number of outliers and
significantly outperform previous higher-order as well as
second-order methods in terms of both matching score
and matching accuracy. The better performance of higherorder approaches shows their robustness against geometric
transformations compared to second-order ones. In the first
experiment without scaling as shown in Fig. 3 (a), most
higher-order methods and MPM perform very well without being deteriorated much by outliers. When the scale
change is present in the data as shown in Fig. 3(b), our
algorithms significantly outperform the current state-of-theart. In particular, both MPM and IPFP fail in this case since
they use pairwise distance similarities which are clearly not
invariant under scaling. Previous higher-order algorithms
like RRWHM or TM also show quite rapidly decreasing
performance as the number of outliers increases compared
to our methods. This experiment also features the difference
between two groups of our algorithms, i.e. the ones using
MPM as subroutine including BCAGM+MP and AdaptBCAGM3+MP perform the best in both settings. Thus, one
can state that our approach transfers the robustness of MPM
as a second-order method to the higher-order case and thus
is able to deal with scaling which is difficult for secondorder methods.

1

Accuracy

2

80

Baseline

(c) 10 pts vs. 30 pts

100

300
200
100
0

20

40

60

80

100

Baseline

(d) 10 pts vs. 30 pts

Fig. 5: CMU house dataset with different number of points in two images.
(Best viewed in color.)

For the deformation setting, σ is varied from 0 to challenging 0.4 while the number of inliers nin is fixed to 20.
Since we want to test the tolerance of all methods only with
respect to deformation, we do not add any outliers nor scale
the points in this setting. This type of experiment has been
used in previous works [12], [25], where σ is varied from 0
to 0.2. Fig. 4 shows the result where our algorithms are quite
competitive to RRWHM which has been shown to be robust
to deformation before [25]. It is interesting to note that while
MPM is very robust against outliers (without scaling), it
is outperformed by other methods when the deformation
is large. Notably, our algorithms BCAGM+MP and AdaptBCAGM3+MP using MPM as subroutine are not affected by
this slight weakness of MPM.
Considering all experiments in this section, our algorithms outperform or are on par with all previous approaches in terms of robustness to outliers, deformation
and scaling. Moreover, our proposed third-order algorithms
show on average quite similar performance compared to
their fourth-order counterparts [36].
5.3

CMU House Dataset

The CMU house dataset has been widely used in previous
work [10], [15], [25], [50] to evaluate matching algorithms.
In this dataset, 30 landmark points are manually tracked
over a sequence of 111 images, which are taken from the
same object under different view points. In this experiment,
“baseline” denotes the distance of the frames in the sequence and thus correlates well with the difficulty to match
the corresponding frames.
We match all possible image pairs with “baseline” of
10, 20, 30, . . . , 100 frames and compute the average matching accuracy for each algorithm. The algorithms are evaluated in three settings. In the first experiment, we match
30 points to 30 points. Then we make the problem significantly harder by randomly removing points from one

-1.5

-2

-0.5

TM
RRWHM
BCAGM
BCAGM+MP
BCAGM3
BCAGM3+MP
Adapt-BCAGM3
Adapt-BCAGM3+MP

Seconds (log10)

Seconds (log10)

11

-2.5

20

40

60

80

-1
-1.5
-2
-2.5

100

Baseline

(a) House

0

5

10

15

# Outliers

(b) Willow dataset

Fig. 6: Running time of higher-order methods (Best viewed in color.)

image, motivated by a scenario where one has background
clutter and thus not all points can be matched. This results
in two matching experiments, namely 10 points to 30 points,
and 20 points to 30 points. For the choice of σs in the affinity
tensor for second-order methods, we follow [10], [50] by
setting σs = 2500.
The experimental results are shown in Fig. 5. While most
algorithms perform rather well on the 30 to 30 task, our
methods perform significantly better than all other methods
in the more difficult tasks, thus showing as for the synthetic
datasets that our methods are quite robust to different kinds
of noise in the matching problem. Notably, in the 10 to 30
task, Adapt-BCAGM3+MP outperforms BCAGM+MP for
large baselines.
5.4

Willow Object Class Dataset

In this last experiment we use the Willow Object Class
dataset [1] previously introduced in Section 5.1. We use exactly the same setting as in Section 5.1. We test all algorithms
on 2000 randomly selected image pairs from each object
class, where the number of outliers in one image is varied
from 0 to 20. For a given number of outliers, we report in
Fig. 7 the average performance of all algorithms over all the
2000 image pairs. It can be seen that both our third-order
and fourth-order methods consistently outperform the stateof-the-art for all the object classes even so the differences
are a bit smaller. In particular, for the Duck class, AdaptBCAGM3+MP performs slightly better than BCAGM+MP
even though this is hard to see (see the Appendices for
detailed comparisons) because the result has been averaged
already over all pairs of images.
5.5 Comparison between third-order methods and
fourth-order methods
This section presents a detailed comparison of our thirdorder approaches and fourth-order approaches. In particular, we first compare the non-adaptive third-order algorithms (Algorithm 1 and Algorithm 2) to their non-adaptive
fourth-order counterparts [36]. Second, we extend the adaptive scheme derived for third-order tensors in Section
3.3 to fourth-order tensors, and denote the corresponding
adaptive fourth-order algorithms as Adapt-BCAGM, AdaptBCAGM+IPFP and Adapt-BCAGM+MP with the same notation as before. We then compare these algorithms with
the adaptive third-order methods from Algorithm 3 and Algorithm 4. The purpose is to test for advantages and disad-

vantages of third-order approaches versus fourth-order ones
under both settings – with and without adaptive updates.
For each pair of algorithms, we report in Table 1 the
number of test cases where one algorithm performs better than the other one in terms of matching score and
how much improvement it achieves on average over the
other method for the cases where the performance is not
equal. More precisely, let fi and gi be the matching score
of BCAGM3 and BCAGM respectively on the Synthetic
dataset. Then “BCAGM3 > BCAGM” denotes the cases
where BCAGM3 achieves a better objective value than
BCAGM, and we report this number in the first column
(No.), that is, |{i | fi > gi }| , and the average performance
gain in the second column (Avg(%)), that is,
X fi − gi
.
|{i | fi > gi }| i:f >g
gi

1

i

(14)

i

Table 1 shows that the standard fourth-order algorithms
[36] have a better overall performance than the standard
third-order ones in terms of matching score. However, the
differences are minor and obtaining better matching score
does not necessarily imply the same for matching accuracy,
especially when the ground truth is not a global optimum
of the hypergraph matching problem. Fig. 7 (a) shows an
example of this phenomenon where third-order methods
achieve better matching score than the fourth-order algorithms but their matching accuracy is very similar.
Table 2 shows the same kind of comparison but for
the adaptive schemes. One can see that the situation has
changed as the adaptive third-order schemes win more
frequently over their fourth-order counterparts. Again, the
differences are minor. While on the synthetic dataset the
adaptive fourth-order methods perform better, on the real
dataset, it is the opposite.
All in all, one can state that the standard fourth-order
algorithms are better than the standard third-order approaches in terms of achieving better matching score. This
indicates potential benefits of using fourth-order tensors
as well as the lifting step as done in [36]. However, the
differences between the adaptive third- and fourth-order
methods are negligible and are slightly on the side of the
adaptive third-order methods. As the third-order methods
are two times faster (see next section), our recommendation
is to use the adaptive third-order methods in practice.

5.6

Running time

Fig. 6 shows the running time of all hypergraph matching
algorithms. All our methods achieve competitive/smaller
running time than the state-of-the-art, in particular, they
are 4 to 16 times faster than TM while being competitive
to RRWHM in case of BCAGM3+MP or up to two times
faster in case of BCAGM3 and Adapt-BCAGM3. It can also
be seen that the proposed third-order approaches are two
times faster than their fourth-order counterparts on average.
Further experimental results on the running time of all
higher-order algorithms can be found in Tab. 6, Tab. 7 and
Tab. 8 in the appendix.

12

TM
RRWHM
BCAGM
BCAGM+MP
BCAGM3
BCAGM3+MP
Adapt-BCAGM3
Adapt-BCAGM3+MP

160

150

100

250

Matching Score

Matching Score

Matching Score

200

140
120
100

200

150

80
50
0

5

10

15

0

20

5

15

20

0

1

0.8

0.8

0.8

0.6
0.4
0.2

Accuracy

1

0.6

0

0.4

10

15

20

15

20

15

20

0.4
0.2

0
5

10

0.6

0.2

0

5

# Outliers

1

Accuracy

Accuracy

10

# Outliers

# Outliers

0
0

5

# Outliers

10

15

20

0

5

# Outliers

(a) Face

10

# Outliers

(b) Duck

(c) Winebottle

Fig. 7: Matching score and accuracy of higher-order algorithms on Willow Object Dataset [9]. (Best viewed in color.)
Synthetic

House

Face+Duck+Winebottle

Comparison

No.

Avg(%)

No.

Avg(%)

No.

Avg(%)

BCAGM3 > BCAGM

241

51.9

2

1.5

6277

30.0

BCAGM3 < BCAGM

292

347.0

7

18.5

7130

27.0

BCAGM3 = BCAGM

4567

0.0

1671

0.0

49593

0.0

BCAGM3+IPFP > BCAGM+IPFP

186

2.4

19

7.5

8960

24.0

BCAGM3+IPFP < BCAGM+IPFP

188

59.5

5

6.9

9834

41.0

BCAGM3+IPFP = BCAGM+IPFP

4726

0.0

1656

0.0

44206

0.0

BCAGM3+MP > BCAGM+MP

194

43.4

4

12.2

5695

32.0

BCAGM3+MP < BCAGM+MP

243

237.1

6

27.6

7461

56.0

BCAGM3+MP = BCAGM+MP

4663

0.0

1670

0.0

49844

0.0

TABLE 1: Comparison between the standard third-order methods (BCAGM3, BCAGM3+Ψ) and the standard fourth-order methods (BCAGM,
BCAGM+Ψ [36]) on different datasets. The first column (No.) shows the number of test cases where one algorithm achieves better matching
score than the other. The second column (Avg(%)) shows the average performance gain achieved by the better method in each case, which is
computed using Eq. (14).
Synthetic

House

Comparison

No.

Avg(%)

No.

Adapt-BCAGM3 > Adapt-BCAGM

269

47.2

Adapt-BCAGM3 < Adapt-BCAGM

239

51.6

Adapt-BCAGM3 = Adapt-BCAGM

4592

Adapt-BCAGM3+IPFP > Adapt-BCAGM+IPFP

Face+Duck+Winebottle

Avg(%)

No.

Avg(%)

3

1.6

6717

28.0

3

20.3

5621

23.0

0.0

1674

0.0

50662

0.0

207

2.3

20

6.9

10774

21.0

Adapt-BCAGM3+IPFP < Adapt-BCAGM+IPFP

109

55.6

3

1.9

6318

26.0

Adapt-BCAGM3+IPFP = Adapt-BCAGM+IPFP

4784

0.0

1657

0.0

45908

0.0

Adapt-BCAGM3+MP > Adapt-BCAGM+MP

217

39.0

4

12.2

6311

29.0

Adapt-BCAGM3+MP < Adapt-BCAGM+MP

177

59.0

4

10.4

5603

26.0

Adapt-BCAGM3+MP = Adapt-BCAGM+MP

4706

0.0

1672

0.0

51086

0.0

TABLE 2: Comparison between the adaptive third-order methods (Adapt-BCAGM3, Adapt-BCAGM3+Ψ) and the adaptive fourth-order methods
(Adapt-BCAGM, Adapt-BCAGM+Ψ) on different datasets. The first column (No.) shows the number of test cases where one algorithm achieves
better matching score than the other. The second column (Avg(%)) shows the average performance gain achieved by the better method in each
case, which is computed using Eq. (14).

13

6

D EMO OF MATCHING RESULTS

We show in this section a demo matching on realistic images. Different from the previous sections, we run a SIFT
detector on two realistic images, hence, there might be
outliers in both point sets. Since each point on one image
might not have a corresponding point on the other image,
the matching task becomes more difficult. Fig. 8 shows the
matching results. Our algorithms achieve more reasonable
matches than the previous higher-order approaches while
also achieving consistently better matching score than the
others. Further examples can be found in the appendix.

7

C ONCLUSION AND O UTLOOK

In this work, we present a new optimization framework
for hypergraph matching based on a multilinear reformulation. In particular, we extend our theory from fourthorder tensors [36] to third-order tensors, proving under
mild conditions the equivalence between the maximization
of a third-order multilinear form and the maximization
of its associated score function. Based on that, we propose Algorithm 1 and Algorithm 2 which both maximize
the third-order multilinear form directly over the set of
assignment constraints. These algorithms share the same
attributes with their fourth-order counterparts [36] in the
sense that they achieve monotonic ascent in the original
score function directly on the set of assignment matrices,
and both converge to a homogeneous solution in a finite
number of steps. Furthermore, we provide several sufficient
conditions, weaker than the convexity constraint used in
[36], to ensure the equivalence between the optimization
problems. We also propose two adaptive schemes which
yield slight improvements over the standard approaches.
Overall, our algorithms, including the fourth-order and the
third-order ones, not only achieve superior performance
in terms of matching score but also yield competitive or
significantly better matching accuracy than previous work.
This holds even when there are a large number of outliers,
deformation, scaling or other form of noise in the data.
Moreover, our methods are competitive or better in terms
of runtime compared to previous work. Our new adaptive
third-order methods are competitive to the corresponding
adaptive fourth-order ones but are on average two times
faster and thus we recommend them for usage in practice.
For future work it might be interesting to go beyond
injective mappings which require every point on one image
to be mapped to exactly one point in the second image.
Non-injective mappings are of interest if not all points need
to be matched as for instance when outliers are present on
both sides [30], [40]. Our Theorem 3.3 allows the use of new
constraint sets reflecting non-injective mappings, which is
an interesting topic of further research. Another possible
generalization is to go beyond third-order (this paper) or
fourth-order [36] hypergraphs which would require a generalization of Theorem 3.3. At the moment it is unclear to us
if this is possible.

(a) Image 1

(b) Image 2

(c) TM [15] (363)

(d) RRWHM [25] (441)

(e) BCAGM [36] (599)

(f) BCAGM+MP [36] (600)

(g) Adapt-BCAGM3 (600)

(h) Adapt-BCAGM3+MP (600)

Fig. 8: Demo of matching results. Matching score is reported for each
method. (Best viewed in color.)

and Cristian Caloian for helping us to correct grammatical
mistakes.

R EFERENCES
[1]
[2]
[3]
[4]

ACKNOWLEDGMENTS

[5]

This work has been supported by the ERC grant NOLEPRO.
We sincerely thank the reviewers for their helpful comments

[6]

http://www.di.ens.fr/willow/research/graphlearning/. 9, 11
C. Arora and A. Globerson. Higher order matching for consistent
multiple target tracking. In ICCV, 2013. 1
S. Arya, D. M. Mount, N. S. Netanyahu, R. Silverman, and A. Wu.
An optimal algorithm for approximate nearest neighbor searching.
Journal of the ACM, 45:891–923, 1998. 9
E. Boros and P. L. Hammer. Pseudo-boolean optimization. Discrete
Applied Mathematics, 123:155–225, 2002. 2
R. E. Burkhard, M. Dell’Amico, and S. Martello. Assignment
problems. SIAM, Philadelphia, 2012. 5
M. Carcassoni and E. R. Hancock. Spectral correspondence for
point pattern matching. Pattern Recognition, 36:193–204, 2003. 1

14

[7]

[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]

[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]

Z.-Q. Cheng, Y. Chen, M. R.R., Y.-K. Lai, and A. Wang. Supermatching: Feature matching using supersymmetric geometric constraints. IEEE Transactions on Visualization and Computer Graphics,
10:1885–1894, 2012. 1
M. Chertok and Y. Keller. Efficient high order matching. PAMI,
32:2205–2215, 2010. 1, 2
M. Cho, K. Alahari, and J. Ponce. Learning graphs to match. In
ICCV, 2013. 1, 9, 12
M. Cho, J. Lee, and K. M. Lee. Reweighted random walks for
graph matching. In ECCV, 2010. 2, 8, 9, 10, 11
M. Cho and K. M. Lee. Progressive graph matching: Making a
move of graphs via probabilistic voting. In CVPR, 2012. 1
M. Cho, J. Sun, O. Duchenne, and J. Ponce. Finding matches in
a haystack: A max-pooling strategy for graph matching in the
presence of outliers. In CVPR, 2014. 2, 6, 8, 9, 10, 20, 21
D. Conte, P. Foggia, C. Sansone, and M. Vento. Thirty years of
graph matching. PAMI, 18:265–298, 2004. 2
T. Cour, P. Srinivasan, and J. Shi. Balanced graph matching. In
NIPS, 2007. 1, 2, 8, 9
O. Duchenne, F. Bach, I. Kweon, and J. Ponce. A tensor-based
algorithm for high-order graph matching. PAMI, 33:2383–2395,
2011. 1, 2, 3, 8, 9, 10, 13, 20, 21, 22
O. Duchenne, A. Joulin, and J. Ponce. A graph-matching kernel
for object categorization. In ICCV, 2011. 1
A. Filatov, A. Gitis, and I. Kil. Graph-based handwritten digit
string recognition. In Proc. Third Int’l Conf. Document Analysis and
Recognition, 1995. 1
P. Foggia, G. Percannella, and M. Vento. Graph matching and
learning in pattern recognition in the last 10 years. Pattern
Recognition and Artificial Intelligence, 28, 2014. 2
R. Hartley and A. Zisserman. Tutorial i: The algebraic approach to
invariance. In ESPRIT Workshop, 1994. 1
J. Hays, M. Leordeanu, A. A. Efros, and Y. Liu. Discovering texture
regularity as a higher-order correspondence problem. In ECCV,
2006. 1
J. E. Hopcroft and J. K. Wong. Linear time algorithm for isomorphism of planar graphs (preliminary report). Proc. Sixth Ann. ACM
Symp. Theory of Computing, pages 172–184, 1974. 2
V. Kolmogorov and C. Rother. Minimizing nonsubmodular functions with graph cuts-a review. PAMI, 29:1274–1279, 2007. 2
H. W. Kuhn. The hungarian method for the assignment problem.
Naval Research Logistics Quarterly, 2:83–97, 1955. 5
J. Lee, M. Cho, and K. M. Lee. A graph matching algorithm using
data-driven markov chain monte carlo sampling. In ICPR, 2010.
1, 2
J. Lee, M. Cho, and K. M. Lee. Hyper-graph matching via
reweighted random walks. In CVPR, 2011. 1, 2, 8, 9, 10, 13, 20, 21,
22
R. Lee and J. Liu. An oscillatory elastic graph matching model
for recognition of offline handwritten chinese characters. In
Knowledge-Based Intelligent Information Engineering Systems, 1999.
1
M. Leordeanu and M. Hebert. A spectral technique for correspondence problems using pairwise constraints. In ICCV, 2005. 1, 2,
8
M. Leordeanu and M. Hebert. Unsupervised learning for graph
matching. In CVPR, 2009. 1, 9
M. Leordeanu, M. Hebert, and R. Sukthankar. An integer projected
fixed point method for graph matching and map inference. In
NIPS, 2009. 1, 2, 6, 8, 20, 21
W. Lian and L. Zhang. Point matching in the presence of outliers
in both point sets: A concave optimization approach. In CVPR,
2014. 13
Z.-Y. Liu, H. Qiao, and L. Xu. An extended path following
algorithm for graph-matching problem. PAMI, 34:1451–1456, 2012.
2
D. G. Lowe. Object recognition from local scale-invariant features.
In ICCV, 1999. 1
B. Luo and E. Hancock. Alignment and correspondence using
singular value decomposition. Advances in Pattern Recognition,
1876:226–235, 2000. 1
J. Maciel and J. P. Costeira. A global solution to sparse correspondence problems. PAMI, 25:187–199, 2003. 5
J. L. Mundy, A. Zisserman, and D. Forsyth. Applications of Invariance in Computer Vision. Springer, 1993. 1
Q. Nguyen, A. Gautier, and M. Hein. A flexible tensor block
coordinate ascent scheme for hypergraph matching. In CVPR,
2015. 1, 2, 3, 4, 5, 8, 9, 10, 11, 12, 13, 18, 19, 20, 21, 22
C. Schellewald and C. Schnörr. Probabilistic subgraph matching
based on convex relaxation. In EMMCVPR, 2005. 1
A. Sharma, R. Horaud, J. Cech, and E. Boyer. Topologically-

[39]

[40]
[41]
[42]
[43]
[44]

[45]
[46]
[47]
[48]
[49]
[50]

robust 3d shape matching based on diffusion geometry and seed
growing. In CVPR, 2011. 1
R. Singh, J. Xu, and B. Berger. Pairwise global alignment of protein
interaction networks by matching neighborhood topology. Proc.
11th Intl Conf. Research in Computational Molecular Biology, pages
16–31, 2007. 1
Y. Suh and K. M. Lee. Subgraph matching using compactness prior
for robust feature correspondence. In CVPR, 2015. 13
W. Taylor. Protein structure comparison using bipartite graph
matching and its application to protein structure classification.
Molecular and Cellular Proteomics, 1:334–339, 2002. 1
L. Torresani, V. Kolmogorov, and C. Rother. Feature correspondence via graph matching: Models and global optimization. In
ECCV, 2008. 1, 2
H. F. Wang and E. R. Hancock. Correspondence matching using
kernel principal components analysis and label consistency constraints. Pattern Recognition, 39:1012–1025, 2006. 1
Y. Wang, F. Makedon, and J. Ford. A bipartite graph matching
framework for finding correspondences between structural elements in two proteins. In Proc. IEEE 26th Ann. Intl Conf. Eng. in
Medicine and Biology Society, 2004. 1
J. Yan, C. Zhang, H. Zha, W. Liu, X. Yang, and S. M. Chu. Discrete
hyper-graph matching. In CVPR, 2015. 2
M. Zaslavskiy, F. Bach, and J. Vert. A path following algorithm for
the graph matching problem. PAMI, 31:2227–2242, 2009. 1, 2, 5
R. Zass and A. Shashua. Probabilistic graph and hypergraph
matching. In CVPR, 2008. 1, 2, 8
Y. Zeng, C. Wang, Y. Wang, X. Gu, D. Samaras, and N. Paragios. Dense non-rigid surface registration using high-order graph
matching. In CVPR, 2010. 1, 2
F. Zhou and F. De la Torre. Factorized graph matching. In CVPR,
2012. 1, 2, 5, 8
F. Zhou and F. De la Torre. Deformable graph matching. In CVPR,
2013. 1, 2, 5, 9, 10, 11
Quynh Nguyen is currently a PhD student in
the Machine Learning Group of Saarland University. He received his BSc degree in Information Technology from Vietnam National University in Hanoi in 2011. In 2014, he got his
Msc degree in Computer Science from Saarland
University through the support of the IMPRS for
Computer Science in Saarbruecken, Germany.
His research interests are in machine learning,
optimization and computer vision.
Francesco Tudisco is a Postdoctoral fellow
in Mathematics at Saarland University, within
the Machine Learning Group. He received his
diploma and his Ph.D. in Mathematics from the
University of Rome ”Tor Vergata”. His main scientific interests include Operator Theory, Numerical Linear Algebra, Network Analysis and Numerical Optimization.

Antoine Gautier is a PhD student in the machine learning group of Saarland university since
2014. He received a B.Sc. (2011) and a M.Sc.
(2013) degree in applied mathematics from
EPFL, Lausanne, Switzerland. In 2010, he spent
a year at the Humboldt university of Berlin for
an Erasmus exchange. His main research interests are higher-order structures, optimization,
nonlinear spectral theory and their applications
in computational science.
Matthias Hein is professor of mathematics and
computer science at Saarland University. He received his diploma in physics from the University of Tübingen in 2002 and his PhD in computer science from the University of Darmstadt
in 2005. His research interests are in machine
learning, statistics and optimization and the application of machine learning to problems in
bioinformatics, computer vision and other fields
in computer science and the natural sciences.

15

A PPENDIX A

The statement 2) also implies the following inequalities

Theorem 3.3. Let D ⊂ Rn be a compact set, F : Rn × Rn ×
Rn → R a symmetric third-order multilinear form and
S : Rn → R the associated third-order score function
such that S(x) = F (x, x, x), then we have the following
chain of implications: 1) ⇒ 2) ⇒ 3) ⇒ 4).

0 ≤ F(x, y − z, y − z)
0 ≤ F(y, z − x, z − x)
0 ≤ F(z, x − y, x − y)
which can be expanded as

2

1) ∇ S(x) is positive semidefinite for all x ∈ D.
2) F(x, y − z, y − z) ≥ 0 for all x, y, z ∈ D.
3) It holds for all x, y, z ∈ D
F(x, y, z) ≤

max
u∈{x,y,z}

2F(x, y, z) ≤ F(x, y, y) + F(x, z, z)
2F(y, z, x) ≤ F(y, z, z) + F(y, x, x)
2F(z, x, y) ≤ F(z, x, x) + F(z, y, y).

F(u, u, u).

4) The optimization of the multilinear form is equivalent
to the optimization of its associated score function

By summing up these inequalities and taking into account
the symmetry of F, one obtains

6 F(x, y, z) ≤

F(y, z, z) + F(z, x, x) + F(z, y, y). (16)

max F(x, y, z) = max F(x, x, y) = max F(x, x, x).

x,y,z∈D

x,y∈D

x∈D

F(x, y, y) + F(x, z, z) + F(y, x, x) +

Combining (15) and (16) to achieve

6 F(x, y, z) ≤ 2F(x, x, x) + 2F(y, y, y) + 2F(z, z, z),

Proof:

1) =⇒ 2) : We first recall the definition of the score
function as

•

S(x) = F(x, x, x) =

n
X

Fijk xi xj xk .

which implies
F(x, y, z) ≤ max{F(x, x, x), F(y, y, y), F(z, z, z)}.
•

3) =⇒ 4) : From the statement 3)

i,j,k=1

F(x, y, z) ≤

One can compute the Hessian of S at x as follows

∇2 S(x)


ij

=6

n
X

Fijk xk

max
u∈{x,y,z}

F(u, u, u),

one takes the maximum of both sides over D to achieve


max F(x, y, z) ≤
max
max F(u, u, u)

∀ i, j = 1, . . . , n.

k=1

x,y,z∈D

x,y,z∈D

=

Using our notation, the Hessian can be rewritten as

u∈{x,y,z}

max F(x, x, x).
x∈D

However, it holds for any compact set D ⊂ Rn that

2

∇ S(x) = 6 F(x, · , · ).
The Hessian ∇2 S(x) is is positive semidefinite for all x ∈
D if and only if t, ∇2 S(x) t ≥ 0 for all x ∈ D, t ∈ Rn ,
which is equivalent to F(x, t, t) ≥ 0 for all x ∈ D, t ∈ Rn
and thus,

max F(x, y, z) ≥ max F(x, x, y) ≥ max F(x, x, x)

x,y,z∈D

x,y∈D

x∈D

thus, the equality must hold for the whole chain, that is,

max F(x, y, z) = max F(x, x, y) = max F(x, x, x).

x,y,z∈D

x,y∈D

x∈D

F(x, y − z, y − z) ≥ 0 ∀ x, y, z ∈ D.



2) =⇒ 3) : Given x, y, z ∈ D, the statement 2) implies
that F(x, x − y, x − y) ≥ 0 and F(y, x − y, x − y) ≥ 0,
thus,

•

0

≤ F(x, x − y, x − y) + F(y, x − y, x − y)
= F(x, x, x) + F(y, y, y) − F(x, x, y) − F(x, y, y)

By switching the roles of x, y, z, we can derive the following inequalities

≤
≤
F(y, y, z) + F(y, z, z) ≤

F(x, x, y) + F(x, y, y)

F(x, x, x) + F(y, y, y)

F(x, x, z) + F(x, z, z)

F(x, x, x) + F(z, z, z)

Fα (x, y, z) = F(x, y, z) + α G(x, y, z).
qP
n
2
For all α ≥ 27
max
j,k=1 (Fijk ) it holds
4
i=1,2,...,n

1) Fα is a symmetric third-order multilinear form.
2) ∇2 Sα (x) is positive semidefinite for all x ∈ M.
3) The new problem is equivalent to the original one

arg max Sα (x) = arg max S(x).

F(y, y, y) + F(z, z, z)

x∈M

x∈M

Proof:
1) Since F and G are symmetric multilinear forms, Fα must
be symmetric as well.
2) In the following, we define a new tensor Gε as

Summing up these equations gives us
F(x, x, y) + F(x, y, y) + F(x, x, z) +
F(x, z, z) + F(y, y, z) + F(y, z, z)

≤ 2F(x, x, x) + 2F(y, y, y) + 2F(z, z, z).

Proposition 3.4. Let F : Rn × Rn × Rn → R be a symmetric
third-order multilinear form and G : Rn × Rn × Rn → R
as defined in Eq. (8). We consider the new multilinear
form Fα defined by

(15)

Gε =

n
X
i=1

eεi ⊗ eεi ⊗ eεi

16

where eεi = ε1 + (1 − ε)ei for i = 1, 2, . . . , n and 0 <
ε < 1. Let Gε be the associated multilinear form
Gε (x, y, z)

n
X

=

(Gε )pqr xp yq zr
p,q,r=1
n
X
heεi , xi heεi , yi heεi , zi ,
i=1

=

(17)

and the new score function is defined as

Combining this with (18) yields the following sufficient
condition
v
uX
u n
1
t
max
α ≥
(Fijk )2 .
ε(1 − ε)2 i=1,2,...,n j,k=1
In order to minimize this bound, one solves the following one-dimensional optimization problem

min

Sα,ε (x) = F(x, x, x) + α Gε (x, x, x).

0<ε<1

One first observes that the new score function is
related with the original one via Sα = Sα, 1 .
3
For any given 0 < εq< 1, we will prove that
Pn
1
2 is a suffimax
α ≥ ε(1−ε)
2
j,k=1 (Fijk )
i=1,2,...,n
for ∇2 Sα,ε (x)

cient condition
to be positive semidefinite at every x ∈ M. Indeed, ∇2 Sα,ε (x) =
6 (F(x, · , · ) + αGε (x, · , · ))  0 for all x ∈ M if and
only if F(x, y, y) + αGε (x, y, y) ≥ 0 for all x ∈ M and
y ∈ Rn . One has

|F(x, y, y)| =

n
X

which leads to ε = 1/3. Finally, applying the
above result
for ε = 1/3, one obtains that α ≥
qP
n
27
2
max
j,k=1 (Fijk ) is a sufficient condition for
4
i=1,2,...,n

∇2 Sα (x) to be positive semidefinite at every x in M.
3) From Eq. (17), one has the relation G(x, x, x) =
Gε=1/3 (x, x, x). It holds for all x ∈ M
Gε (x, x, x)

=



≤

=

≤

Fijk yj yk 

=

xi

Fijk yj yk

j,k=1

v
uX
u n
2
t
(Fijk )2 kxk1 kyk2

Sα (x)

=
=

i=1
n
X

≥
=

(19)

3

=
=

F(x, x, x) + α G(x, x, x)

S(x) + const,

arg max Sα (x) = arg max S(x).
x∈M

x∈M



Theorem 3.5. Let F be
a symmetric third-order tensor and
qP
n
2
α ≥ 27
max
j,k=1 (Fijk ) . Then the following
4

2

h(1 − ε)ei + ε1, xi heεi , yi

i=1,2,...,n

holds for Algorithm 1:


2
(1 − ε)xi + ε h1, xi heεi , yi

ε h1, xi heεi , yi

2

i=1
n
X

3

2

heεi , xi heεi , yi

i=1
n
X

3

which implies the modified problem is still equivalent
to the original problem, that is,

j,k=1

i=1
n
X

3

(1 − ε)3 kxk3 + nε3 kxk1 +

2

(18)
for all x ∈ M, y ∈ R . Moreover, it holds for all x ∈ M
and y ∈ Rn that

=

3

However, kxk1 = kxk2 = kxk3 = n1 = const for all
x ∈ M , thus Gε (x, x, x) is constant on M for any given
0 < ε < 1. This results in

n

Gε (x, y, y)

(1 − ε)xi + ε kxk1

3(1 − ε)ε2 kxk1 .

j,k=1

Therefore, a sufficient condition is to have
v
uX
u n
2
α Gε (x, y, y) ≥ max t
(Fijk )2 kxk1 kyk2
i=1,2,...,n

n
X

3

j,k=1

max

(1 − ε)xi + ε h1, xi

2

j,k=1

i=1,2,...,n

n
X

3(1 − ε)2 ε kxk2 kxk1 +

v
v
uX
uX
u n
u n
xi t
(Fijk )2 t
yj2 yk2

n
X

heεi , xi heεi , xi heεi , xi

i=1

=

n
X

i=1

≤

n
X

i=1



j,k=1

i=1
n
X

n
X

xi 

i=1
n
X

=

i=1

Fijk xi yj yk

i,j,k=1
n
X

1
ε(1 − ε)2

ε kxk1 (1 − ε)yi + ε h1, yi

2

1) The sequence Fα (xk , yk , zk ) for k = 1, 2, . . . is strictly
monotonically increasing or terminates.
2) The sequence of scores S(um ) for m = 1, 2, . . . is
strictly monotonically increasing or terminates. For every m, um ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations.

i=1

≥ ε(1 − ε

2

2

kxk1 kyk2 .

Proof:

17

1) It follows from the definition of steps 4) − 6) in Algorithm 1,
Fα,k

:=
≤
≤
≤

Fα (xk , yk , zk )
Fα (x̃k+1 , yk , zk )
Fα (x̃k+1 , ỹk+1 , zk )
Fα (x̃

k+1

, ỹ

k+1

, z̃

k+1

2) The sequence of scores S(um ) for m = 1, 2, . . . is
strictly monotonically increasing or terminates. For every m, um ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations.

(20)

) =: F̃α,k+1 .

Either F̃α,k+1 > Fα,k in which case

xk+1 = x̃k+1 , yk+1 = ỹk+1 , zk+1 = z̃k+1
and

Proof: The proof is similar to that of Theorem 3.5.

Proposition 3.8. Let F : Rn × Rn × Rn → R be a symmetric
third-order multilinear form and G : Rn × Rn × Rn → R
as defined in Eq. (8). For each α ∈ R, let Fα be defined as
Fα (x, y, z) = F(x, y, z) + α G(x, y, z).

Fα (xk+1 , yk+1 , zk+1 ) > Fα (xk , yk , zk ),
or F̃α,k+1 = Fα,k and the algorithm enters step 7). Since
the Hessian of Sα is positive semidefinite at every point
on M for the chosen value of α, we get by the inequality
(6) in Theorem 3.3
F̃α,k+1

=
=

If the strict inequality does not hold then the termination condition at step 12) of the algorithm is met.
Otherwise, we get strict ascent
F̃α,k+1

Then the following holds
1) Fα is a symmetric multilinear form.
2) For all non-homogeneous tuple (x, y, z) ∈ M (i.e. x 6=
y or y 6= z or z 6= x), the following inequality holds

= Fα (x̃k+1 , ỹk+1 , z̃k+1 )
≤
max
Fα (v, v, v)
v∈{x̃k+1 ,ỹk+1 ,z̃k+1 }
Fα (um+1 , um+1 , um+1 )
Sα (um+1 ).



Fα (x, y, z) ≤

max
u∈{x,y,z}

Fα (u, u, u)

if and only if
F(x, y, z) −

α ≥ Λ(x, y, z) :=

max
u∈{x,y,z}

F(u, u, u)

G(x, x, x) − G(x, y, z)

.

3) For every α ∈ R, it holds

< Fα (um+1 , um+1 , um+1 )
= Sα (um+1 ) = Fα (xk+1 , yk+1 , zk+1 ).

arg max Sα (x) = arg max S(x).
x∈M

x∈M

Proof:

2) From

Sα (um+1 )

=
=

Fα (um+1 , um+1 , um+1 )
Fα (xk+1 , yk+1 , zk+1 ),

it follows that Sα (um ) for m = 1, 2, . . . is a subsequence of Fα (xk , yk , zk ) for k = 1, 2, . . . and thus it
holds either Sα (um ) = Sα (um+1 ) in which case the
algorithm terminates or Sα (um ) < Sα (um+1 ). However, by Proposition 3.4, the additional term which has
been added to S is constant on M , that is Sα (x) =
S(x) + const for all x ∈ M. It follows that either
S(um ) = S(um+1 ) and the algorithm terminates or
S(um ) < S(um+1 ).
3) The algorithm yields a strictly monotonically increasing
sequence S(um ), m = 1, 2, . . . or it terminates. Since
there is only a finite number of possible assignment
matrices, the algorithm has to terminate after a finite
number of iterations.

Theorem 3.7. Let F q
be a symmetric third-order tensor and
Pn
27
2
α ≥ 4 max
j,k=1 (Fijk ) . Let Ψ be an algorithm
i=1,2,...,n

for the QAP which yields monotonic ascent w.r.t. the
current iterate. Then the following holds for Algorithm
2:
1) The sequence Fα (xk , yk , yk ) for k = 1, 2, . . . is strictly
monotonically increasing or terminates.

1) Since G is a symmetric multilinear form, Fα must be
symmetric as well.
2) It holds for all x, y, z ∈ M
Fα (x, y, z) ≤

⇔ Fα (x, y, z) ≤
⇔ Fα (x, y, z) ≤

max
u∈{x,y,z}

max
u∈{x,y,z}

max
u∈{x,y,z}

Fα (u, u, u)

{F(u, u, u) + α G(u, u, u)}
{F(u, u, u)} + α G(x, x, x).

Note that we have used in the last step the fact that
G(u, u, u) = G(x, x, x) = G(y, y, y) = G(z, z, z) for all
x, y, z ∈ M (see Eq. (19)). By substituting Fα (x, y, z) =
F(x, y, z) + α G(x, y, z), the above inequality chain is
equivalent to

α(G(x, x, x) − G(x, y, z)) ≥ F(x, y, z) − max F(u, u, u).
u∈{x,y,z}

Notice that G(x, x, x) > G(x, y, z) (see proof below) for
every non-homogeneous tuple {x, y, z} ∈ M, thus, one
can divide both sides by (G(x, x, x) − G(x, y, z)) to get
F(x, y, z) −

α≥

max
u∈{x,y,z}

F(u, u, u)

G(x, x, x) − G(x, y, z)

.

**) Assume {x, y, z} ∈ M is a non-homogeneous tuple,
we will prove for any 0 < ε < 1 that
Gε (x, x, x) > Gε (x, y, z).

(21)

18

be attained at z = u. Combining (22) and (23), one obtains
for all α ≥ α0

From Eq. (17), one has

:=
=

Gε (x, y, z)
n
X
heεi , xi heεi , yi heεi , zi
i=1
n
X

u

=

arg max {Fα0 (u, u, z) + (α − α0 ) G(u, u, z)}

=

arg max {F(u, u, z) + α G(u, u, z)}

=

arg max Fα (u, u, z).

z∈M
z∈M



(1 − ε)xi + ε h1, xi · (1 − ε)yi + ε h1, yi ·

z∈M

i=1





=

(1 − ε)zi + ε h1, zi
n
X


(1 − ε)xi + εn1 · (1 − ε)yi + εn1 ·
i=1

(1 − ε)zi + εn1
=

A PPENDIX B



3(1 − ε)ε2 n31 + (1 − ε)3

n
X

B.1

xi yi zi +

Table 3 compares non-adaptive methods from Algorithm 1
and Algorithm 2 with adaptive methods from Algorithm 3
and Algorithm 4. Table 4 shows the same comparison for
fourth-order algorithms.

i=1

(1 − ε)2 εn1

n
X

(xi yi + yi zi + xi zi ) + nε3 n31 .

i=1

P
Pn
xi yi zi ≤ P ni=1 xi xi xi ,
Since
Pn x, y, z ∈ M
Pn, one has Pi=1
n
xi , ni=1 yi ziP≤
i yi ≤
i=1 x
i=1 xiP
i=1 yi yi =
P
n
n
n
n
=
x
x
,
and
x
z
≤
z
zi = n 1 =
1
i
i
i
i
i
i=1
i=1
i=1
Pn
i=1 xi xi . Therefore, Gε (x, y, z) ≤ Gε (x, x, x). The
inequality becomes equality if and only if x = y = z.
Applying this result for ε = 1/3 finishes the proof.
3) It follows from Eq. (19) that G(x, x, x) is constant for all
x ∈ M. Thus, it holds for every α ∈ R and x ∈ M that

Sα (x)

arg max Sα (x) = arg max S(x).
x∈M


Lemma 3.6. Suppose Algorithm 1 runs with α = α0 for some
α0 ∈ R+ . If x̃k+1 = ỹk+1 = z̃k+1 =: u holds at some
iteration k then for all α ≥ α0 , it holds

u = arg max Fα (z, u, u).
z∈M

Proof: Recall from Algorithm
1 the multilinear form G
Pn
defined as G(x, y, z) = i=1 hēi , xi hēi , yi hēi , zi , and one
has the relation
Fα (x, y, z) = F(x, y, z) + α G(x, y, z).
By substituting x̃k+1 = ỹk+1 = z̃k+1 = u and α = α0 into
step (6) of Algorithm 1, we have

=

arg max Fα0 (u, u, z).

(22)

z∈M

Besides, one can show that

u = arg max G(u, u, z).

Table 5 shows the comparison between the proposed adaptive third-order methods with the state-of-the-art fourthorder approaches [36].
B.3

B.4

which implies

u

B.2 Adaptive third-order methods vs. Standard fourthorder methods

Running time

Tab. 6, Tab. 7 and Tab. 8 show the running time of all higherorder algorithms in seconds on the previous experiments.

:= Fα (x, x, x)
= F(x, x, x) + α G(x, x, x)
= S(x) + const,

x∈M

Adaptive methods vs Non-adaptive methods

(23)

z∈M

Indeed, it follows from Eq. (21) that G(u, u, u) > G(u, u, z)
for all z ∈ M and z 6= u, which implies the maximum must

Demo of matching results

Fig. 9 and Fig. 10 show the examples of matching results on
CMU House dataset. Fig. 11 and Fig. 12 shows the examples
of matching results on realistic images where our algorithms
return more reasonable matches and also achieve higher
matching scores than the state-of-the-art higher-order approaches.

19

Synthetic

House

Face+Duck+Winebottle

Comparison

No.

Avg(%)

No.

Avg(%)

No.

Avg(%)

Adapt-BCAGM3 > BCAGM3

185

321.9

7

10.4

4141

16.0

Adapt-BCAGM3 < BCAGM3

0

0.0

0

0.0

6

3.0

Adapt-BCAGM3 = BCAGM3

4915

0.0

1673

0.0

58853

0.0

Adapt-BCAGM3+IPFP > BCAGM3+IPFP

167

20.1

3

9.9

7566

32.0

Adapt-BCAGM3+IPFP < BCAGM3+IPFP

0

0.0

0

0.0

0

0.0

Adapt-BCAGM3+IPFP = BCAGM3+IPFP

4933

0.0

1677

0.0

55434

0.0

Adapt-BCAGM3+MP > BCAGM3+MP

168

146.7

2

61.9

4234

51.0

Adapt-BCAGM3+MP < BCAGM3+MP

0

0.0

0

0.0

0

0.0

Adapt-BCAGM3+MP = BCAGM3+MP

4932

0.0

1678

0.0

58766

0.0

TABLE 3: Adaptive third-order methods (Adapt-BCAGM3, Adapt-BCAGM3+Ψ) vs. non-adaptive third-order methods (BCAGM3, BCAGM3+Ψ). The
first column (No.) shows the number of test cases where one algorithm achieves better matching score than the other. The second column (Avg(%))
shows the average performance gain achieved by the better method in each case, which is computed using Eq. (14).

Synthetic

House

Face+Duck+Winebottle

Comparison

No.

Avg(%)

No.

Avg(%)

No.

Avg(%)

Adapt-BCAGM > BCAGM

52

4.7

2

0.8

1448

8.0

Adapt-BCAGM < BCAGM

0

0.0

0

0.0

4

13.0

Adapt-BCAGM = BCAGM

5048

0.0

1678

0.0

61548

0.0

Adapt-BCAGM+IPFP > BCAGM+IPFP

53

8.1

1

5.1

2130

14.0

Adapt-BCAGM+IPFP < BCAGM+IPFP

0

0.0

0

0.0

0

0.0

Adapt-BCAGM+IPFP = BCAGM+IPFP

5047

0.0

1679

0.0

60870

0.0

Adapt-BCAGM+MP > BCAGM+MP

50

4.4

0

0.0

876

8.0

Adapt-BCAGM+MP < BCAGM+MP

0

0.0

0

0.0

0

0.0

Adapt-BCAGM+MP = BCAGM+MP

5050

0.0

1680

0.0

62124

0.0

TABLE 4: Adaptive fourth-order methods (Adapt-BCAGM, Adapt-BCAGM+Ψ) vs. non-adaptive fourth-order methods (BCAGM, BCAGM+Ψ [36]).

Synthetic

House

Face+Duck+Winebottle

Comparison

No.

Avg(%)

No.

Avg(%)

No.

Avg(%)

Adapt-BCAGM3 > BCAGM

295

44.2

5

1.3

7467

26.0

Adapt-BCAGM3 < BCAGM

233

52.8

3

20.3

5394

24.0

Adapt-BCAGM3 = BCAGM

4572

0.0

1672

0.0

50139

0.0

Adapt-BCAGM3+IPFP > BCAGM+IPFP

239

3.7

20

7.1

11699

21.0

Adapt-BCAGM3+IPFP < BCAGM+IPFP

99

61.1

3

1.9

5854

28.0

Adapt-BCAGM3+IPFP = BCAGM+IPFP

4762

0.0

1657

0.0

45447

0.0

Adapt-BCAGM3+MP > BCAGM+MP

246

35.2

4

12.2

6722

28.0

Adapt-BCAGM3+MP < BCAGM+MP

169

61.4

4

10.4

5491

26.0

Adapt-BCAGM3+MP = BCAGM+MP

4685

0.0

1672

0.0

50787

0.0

TABLE 5: Adaptive third-order methods (Adapt-BCAGM3, Adapt-BCAGM3+Ψ) vs. non-adaptive fourth-order methods (BCAGM, BCAGM+Ψ [36]).

nout = 0

40

80

120

160

200

TM

0.0582

0.4969

0.9056

1.3159

1.7609

2.2812

RRWHM

0.0066

0.0537

0.1345

0.2255

0.3912

0.4648

BCAGM

0.0039

0.0552

0.1465

0.2674

0.4280

0.6756

BCAGM+MP

0.0057

0.0747

0.2009

0.4515

0.8428

1.5472

BCAGM3

0.0020

0.0247

0.0636

0.1178

0.2000

0.3327

BCAGM3+MP

0.0025

0.0344

0.0972

0.2238

0.4309

0.8144

Adapt-BCAGM3

0.0019

0.0247

0.0636

0.1182

0.2019

0.3421

Adapt-BCAGM3+MP

0.0025

0.0346

0.0981

0.2256

0.4412

0.8355

TABLE 6: Average running time (in seconds) of higher-order algorithms on the synthetic outlier settings. (see Fig. 3 for their matching score and
matching accuracy.)

20

baseline = 10

20

30

40

50

60

70

80

90

100

TM

0.0401

0.0404

0.0406

0.0406

0.0415

0.0419

0.0421

0.0412

0.0405

0.0412

RRWHM

0.0021

0.0022

0.0022

0.0023

0.0033

0.0054

0.0085

0.0091

0.0107

0.0083

BCAGM

0.0024

0.0024

0.0024

0.0025

0.0026

0.0027

0.0026

0.0025

0.0024

0.0024

BCAGM+MP

0.0038

0.0039

0.0038

0.0039

0.0043

0.0044

0.0044

0.0043

0.0042

0.0042

BCAGM3

0.0013

0.0013

0.0013

0.0013

0.0014

0.0013

0.0013

0.0013

0.0012

0.0013

BCAGM3+MP

0.0017

0.0017

0.0017

0.0017

0.0018

0.0019

0.0020

0.0019

0.0018

0.0019

Adapt-BCAGM3

0.0013

0.0013

0.0012

0.0013

0.0013

0.0013

0.0013

0.0013

0.0012

0.0012

Adapt-BCAGM3+MP

0.0016

0.0017

0.0017

0.0017

0.0019

0.0019

0.0020

0.0019

0.0018

0.0019

TABLE 7: Average running time (in seconds) of higher-order algorithms on the CMU House dataset. (see Fig. 5 for their matching score and
matching accuracy.)

nout = 0

5

10

15

TM

0.0914

0.1885

0.3219

0.3567

RRWHM

0.0220

0.0499

0.2326

0.1682

BCAGM

0.0062

0.0123

0.0212

0.0285

BCAGM+MP

0.0098

0.0247

0.0490

0.0654

BCAGM3

0.0032

0.0061

0.0103

0.0142

BCAGM3+MP

0.0047

0.0117

0.0232

0.0343

Adapt-BCAGM3

0.0031

0.0059

0.0101

0.0140

Adapt-BCAGM3+MP

0.0046

0.0119

0.0246

0.0370

TABLE 8: Average running time (in seconds) of higher-order algorithms on Willow dataset. (see Fig. 7 for their matching score and matching
accuracy.)

(a) Input: 10 pts vs 30 pts, baseline = 50

(b) IPFP [29] (44.8)

(c) MPM [12] (112.2)

(d) TM [15] (217.6)

(e) RRWHM [25] (285.9)

(f) BCAGM [36] (329.6)

(g) BCAGM+MP [36] (329.6)

(h) Adapt-BCAGM3 (329.6)

(i) Adapt-BCAGM3+MP (329.6)

Fig. 9: Demo of matching results on CMU House dataset with small baseline (baseline = 50). a) Input images. Yellow dots denote inlier points, blue
dots denote outlier points. b) c) Matching results of previous second-order methods. d) e) Matching results of previous higher-order methods. f) g)
h) i) Matching results of our higher-order approaches. The yellow/red lines indicate correct/incorrect matches. Matching score is reported for each
method. (Best viewed in color.)

21

(a) Input: 10 pts vs 30 pts, baseline = 80

(b) IPFP [29] (86.6)

(c) MPM [12] (34)

(d) TM [15] (80.1)

(e) RRWHM [25] (108.8)

(f) BCAGM [36] (194.4)

(g) BCAGM+MP [36] (194.4)

(h) Adapt-BCAGM3 (194.4)

(i) Adapt-BCAGM3+MP (194.4)

Fig. 10: Demo of matching results on CMU House dataset with large baseline (baseline = 80). a) Input images. Yellow dots denote inlier points,
blue dots denote outlier points. b) c) Matching results of previous second-order methods. d) e) Matching results of previous higher-order methods.
f) g) h) i) Matching results of our higher-order approaches. The yellow/red lines indicate correct/incorrect matches. Matching score is reported for
each method. (Best viewed in color.)

(a) Image 1

(b) Image 2

(c) TM [15] (135)

(d) RRWHM [25] (174)

(e) BCAGM [36] (195)

(f) BCAGM+MP [36] (226)

(g) Adapt-BCAGM3 (218)

(h) Adapt-BCAGM3+MP (226)

Fig. 11: Demo of matching results. Matching score is reported for each method. (Best viewed in color.)

22

(a) Image 1

(b) Image 2

(c) TM [15] (1071)

(d) RRWHM [25] (1639)

(e) BCAGM [36] (1751)

(f) BCAGM+MP [36] (1751)

(g) Adapt-BCAGM3 (1751)

(h) Adapt-BCAGM3+MP (1751)

Fig. 12: Demo of matching results. Matching score is reported for each method. (Best viewed in color.)

