{'fraction_non_alphanumeric': 0.04675239577360963, 'fraction_numerical': 0.015906298632156604, 'mean_word_length': 4.129915966386554, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Efficient Multilinear Optimization Framework for Hypergraph Matching

Hypergraph matching has recently become a popular approach for solving correspondence problems in computer vision as it allows to integrate higher-order geometric information. Hypergraph matching can be formulated as a third-order optimization problem subject to the 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.

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-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 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.

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 x∈M x T Ax while the hypergraph matching problem is formulated as a higher-order assignment problem (HAP) max x∈M n 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.

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.

HYPERGRAPH MATCHING FORMULATION
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 be two point sets with n 1 = |V | ≤ n 2 = |V |. The matching problem consists of finding a binary assignment matrix X ∈ {0, 1} n1×n2 such that X ij = 1 if v i ∈ V matches v j ∈ V , and X ij = 0 otherwise. The set of all possible assignment matrices that assign each vertex of V to exactly one vertex of V is given by Fig. 1: Illustration of constructing the affinity tensor. Each entry is computed by comparing the angles of two triangles formed by three candidate correspondences (i 1 , i 2 ), (j 1 , j 2 ) and (k 1 , k 2 ).
Suppose that we have a function F : j3) . For example, F is calculated as the similarity of the corresponding angles of the two triangles i 1 j 1 k 1 and i 2 j 2 k 2 as shown in Fig. 1. Then, the score S(X) of a matching X ∈ M is defined as follows [15] For ease of notation, we introduce a linear ordering in V ×V so that each X ∈ M can be rewritten as a vector x ∈ {0, 1} n with n = n 1 n 2 and F as a tensor in R n×n×n . With this convention, the score function S : R n → R can be computed in the following way (1) We recall that a third-order tensor T ∈ R n×n×n is called symmetric if its entries T ijk are invariant under any permutation of {i, j, k}. In particular, the tensor T ijk := x i x j x k 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 symmetric 1 . In this paper, we consider only terms of order 3, i.e. F ijk = 0 if i = j, i = k or j = k. However, note that it is possible to integrate first-order terms on the main diagonal F iii , and pairwise potentials F ijj for i = 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.

MATHEMATICAL FOUNDATIONS FOR MULTILIN-EAR OPTIMIZATION FRAMEWORK
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 : R n × R n × R n → R which is defined as 1. Note that F can always be symmetrized (without changing the score function) as followsF ijk = 1 3!
. . , n and S 3 is the permutation group of three elements.
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(x 1 , x 2 , x 3 ) = F(x σ(1) , x σ(2) , x σ(3) ) for all x 1 , x 2 , x 3 ∈ R n and σ a permutation of {1, 2, 3}. Note that if the tensor F is symmetric then F is also symmetric. In the following, we write F( · , y, z) to denote a vector in R n such that F( · , y, z) i = n j,k=1 F ijk y j z k for all 1 ≤ i ≤ n, and write F( · , · , z) to denote a matrix in R n×n such that F( · , · , z) ij = n k=1 F ijk z k 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, In [36] we propose to solve instead the equivalent problem defined as where the new symmetric fourth-order tensor F 4 is obtained by lifting the third order tensor F via 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 R n . Theorem 3.1 (Nguyen et al. [36]). Let F 4 be a symmetric fourth-order tensor and suppose the corresponding score function S 4 : R n → R is convex. Then it holds for any compact constraint set D ⊂ R n , = max x,y,z,t∈D 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 F 4 (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 in [36] as non-trivial third-order score functions are not convex on R n as the following lemma shows. (1). If S is not constant zero, then S : R n → 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 R n 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 ⊂ R n 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 ⊂ R n be a compact set, F : R n × R n × R n → R a symmetric third-order multilinear form and S : R n → R the associated third-order score function such that S(x) = F (x, x, x), then we have the following chain of implications: 3) It holds for all x, y, z ∈ D F(x, y, z) ≤ max u∈{x,y,z} F(u, u, u).
4) The optimization of the multilinear form is equivalent to the optimization of its associated score function max x,y,z∈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 R n . 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 new score function defined as S α (x) = S(x) + α x 2 2 , 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 A homogeneous modification would be S α (x) = S(x) + α x   3 . In this paper, we consider a family of modifications which includes the one above by perturbing the set of basis vectors asē i = ε1+(1−ε)e i with 1 being the all-ones vector. Note that x 3 3 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ē i = 1 3 1 + 2 3 e i we define the symmetric third-order tensor G as and the associated multilinear form as Proposition 3.4. Let F : R n × R n × R n → R be a symmetric third-order multilinear form and G : R n × R n × R n → 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).
3) The new problem is equivalent to the original one arg max x∈M S α (x) = arg max x∈M S(x).

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 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 x,y,z∈M which leads to Algorithm 1 and the other for solving which leads to Algorithm 2. Both variants come along with theoretical guarantees including strict monotonic ascent and finite convergence.

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 α ≥ 27 Then the following holds for Algorithm 1: 1) The sequence F α (x k , y k , z k ) for k = 1, 2, . . . is strictly monotonically increasing or terminates.
2) The sequence of scores S(u m ) for m = 1, 2, . . . is strictly monotonically increasing or terminates. For every m, u m ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations. 10: 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 get no further ascent and only then we set α = 27 4 max i=1,2,...,n n j,k=1 (F ijk ) 2 . It turns 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 + . Ifx k+1 =ỹ k+1 =z k+1 =: u holds at some iteration k then for all α ≥ α 0 , it holds 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).

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 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, where y k is the current iterate and the nonnegative symmetric A ∈ R n×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 be a symmetric third-order tensor and α ≥ 27 Let Ψ be an algorithm for the QAP which yields monotonic ascent w.r.t. the current iterate. The following holds for Algorithm 2: 1) The sequence F α (x k , y k , y k ) for k = 1, 2, . . . is strictly monotonically increasing or terminates.
2) The sequence of scores S(u m ) for m = 1, 2, . . . is strictly monotonically increasing or terminates. For every m, u m ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations.
is a sub-routine method that provides monotonic ascent for the QAP, i.e. z, Az ≥ y k , Ay k 2: Output: else 15: x k+1 =x k+1 , y k+1 =ỹ k+1 16: end if 17: 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.

A Homotopy Tensor Block Coordinate Ascent Scheme
Both algorithms consist of two phases -the first phase uses α = 0 as initialization and the second phase uses α = 27 However, if the value of α 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. x k = y k = z k , 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 : R n × R n × R n → R be a symmetric third-order multilinear form and G : R n × R n × R n → R as defined in Eq. (8). For each α ∈ R, let F α be defined as Then the following holds 1) F α is a symmetric multilinear form.
2) For all non-homogeneous tuple (x, y, z) ∈ M (i.e. x = y or y = z or z = x), the following inequality holds if and only if 3) For every α ∈ R, it holds .
then the inequality (6) is satisfied for all x, y, z ∈ M , for which, it follows from Theorem 3.3 max x,y,z∈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. 10: x k+1 = y k+1 = z k+1 = u m+1 11: m = m + 1 12: else if α does not satisfy Eq. (11) then 13: α = Λ(x k+1 ,ỹ k+1 ,z k+1 ) + ξ 14: x k+1 =x k+1 , y k+1 =ỹ k+1 , z k+1 =z k+1 15: else 16: return u m+1 17: end if 18: else 19: x k+1 =x k+1 , y k+1 =ỹ k+1 , z k+1 =z k+1 20: end if 21: 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 ,ỹ k+1 ,z k+1 ), then inequality (10) holds for the current tuple {x k+1 ,ỹ 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.

COMPUTATIONAL ANALYSIS
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 iteratẽ requires to compute the vector F rst y k s z k t , ∀ r = 1, 2, . . . , n.
Without any knowledge of y k and z k , 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 entries indexed by (r, s, t) for which y k s = 0 or y k t = 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 y k and z k are binary assignment matrices, they have exactly n 1 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 y k s > 0 and z k t > 0, we can access all the related entries F rst using a linked list or hashing method. With that, the above computation can be done in O min{T, n 2 1 K} , where K n = n 1 n 2 is the average number of tensor entries that are associated to each pair (s, t).
For solving the LAP, we use the Hungarian method 2 . The Hungarian method for the LAP has cubic complexity O max{n 1 , n 2 } 3 . Therefore, the total complexity of each iteration is O max{n 1 , n 2 } 3 + min{T, n 2 1 K} . For K = max{n 1 , n 2 } and n 2 1 K < T , the total complexity becomes O max{n 1 , n 2 } 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{n 1 , n 2 } 3 + min{T, n 2 1 K} . For the QAP, the iteratẽ , ∀ s, t = 1, 2, . . . , n.
Adopting the same technique as before, we first pass through all the nonzero components ofx k+1 , and at each index r wherex k+1 r > 0, we access all the associated tensor entries F rst via the use of linked list. This costs O (min{T, n 1 L}), where L (n 1 n 2 ) 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 matrix 3 . Then the computation ofỹ k+1 takes O (t min{T, n 1 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{n 1 , n 2 } 3 + min{T, n 2 1 K} + t min{T, n 1 L}). For L = max{n 1 , n 2 } 2 and n 1 L < T , the whole complexity is simplified to O(t max{n 1 , n 2 } 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.

EXPERIMENTS
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.
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 f i1,i2,i3 and f j1,j2,j3 denote the feature vectors of two triples (i 1 , i 2 , i 3 ) and (j 1 , j 2 , j 3 ), we compute the third-order affinity tensor as: 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 n 1 and n 2 points respectively, then the number of candidate correspondences is n = n 1 n 2 . Thus a third-order affinity tensor consists of n 3 3 entries, which is not only infeasible to compute in reasonable time but also requires huge memory resources even for moderate n 1 and n 2 . Therefore we consider only sparse tensors as in [15], [25]. In particular, we randomly sample n 2 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 where d P i1i2 and d Q j1j2 are the Euclidean distances between two points i 1 , i 2 ∈ P and j 1 , j 2 ∈ Q respectively, σ s is a normalization term that is specified later in each experiment.

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  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].

Synthetic Dataset
In this section we are trying to match point sets in R 2 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 n in 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 n out 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 n in , 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 Adapt-BCAGM3+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. For the deformation setting, σ is varied from 0 to challenging 0.4 while the number of inliers n in 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 Adapt-BCAGM3+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].

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 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.

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, Adapt-BCAGM3+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.

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, Adapt-BCAGM+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 f i and g i 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 | f i > g i }| , and the average performance gain in the second column (Avg(%)), that is, 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. 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.

DEMO 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.

CONCLUSION AND OUTLOOK
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.

ACKNOWLEDGMENTS
This work has been supported by the ERC grant NOLEPRO. We sincerely thank the reviewers for their helpful comments Let D ⊂ R n be a compact set, F : R n × R n × R n → R a symmetric third-order multilinear form and S : R n → 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.

4) The optimization of the multilinear form is equivalent
to the optimization of its associated score function max x,y,z∈D Proof: • 1) =⇒ 2) : We first recall the definition of the score function as One can compute the Hessian of S at x as follows F ijk x k ∀ i, j = 1, . . . , n.
By summing up these inequalities and taking into account the symmetry of F, one obtains 6 F(x, y, z) ≤ F(x, y, y) + F(x, z, z) + F(y, x, x) + F(y, z, z) + F(z, x, x) + F(z, y, y). (16) Combining (15) and (16)  However, it holds for any compact set D ⊂ R n that max x,y,z∈D thus, the equality must hold for the whole chain, that is, Proposition 3.4. Let F : R n × R n × R n → R be a symmetric third-order multilinear form and G : R n × R n × R n → R as defined in Eq. (8). We consider the new multilinear form F α defined by

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 where e ε i = ε1 + (1 − ε)e i for i = 1, 2, . . . , n and 0 < ε < 1. Let G ε be the associated multilinear form G ε (x, y, z) = n p,q,r=1 (17) and the new score function is defined as One first observes that the new score function is related with the original one via S α = S α, 1 3 . For any given 0 < ε < 1, we will prove that α ≥ 1 ε(1−ε) 2 max i=1,2,...,n n j,k=1 (F ijk ) 2 is a sufficient condition for ∇ 2 S α,ε (x) 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 ∈ R n . One has Therefore, a sufficient condition is to have for all x ∈ M, y ∈ R n . Moreover, it holds for all x ∈ M and y ∈ R n that Combining this with (18) yields the following sufficient condition In order to minimize this bound, one solves the following one-dimensional optimization problem min 0<ε<1 1 ε(1 − ε) 2 which leads to ε = 1/3. Finally, applying the above result for ε = 1/3, one obtains that α ≥ 27 4 max i=1,2,...,n n j,k=1 (F ijk ) 2 is a sufficient condition for ∇ 2 S α (x) to be positive semidefinite at every x in M.
3) From Eq. (17), one has the relation G(x, x, However, x 1 = x 2 2 = x 3 3 = n 1 = const for all x ∈ M , thus G ε (x, x, x) is constant on M for any given 0 < ε < 1. This results in which implies the modified problem is still equivalent to the original problem, that is, 1) The sequence F α (x k , y k , z k ) for k = 1, 2, . . . is strictly monotonically increasing or terminates.
2) The sequence of scores S(u m ) for m = 1, 2, . . . is strictly monotonically increasing or terminates. For every m, u m ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations.

Proof:
1) It follows from the definition of steps 4) − 6) in Algorithm 1, EitherF α,k+1 > F α,k in which case orF α,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 If the strict inequality does not hold then the termination condition at step 12) of the algorithm is met. Otherwise, we get strict ascent

2) From
it follows that S α (u m ) for m = 1, 2, . . . is a subsequence of F α (x k , y k , z k ) for k = 1, 2, . . . and thus it holds either S α (u m ) = S α (u m+1 ) in which case the algorithm terminates or S α (u m ) < S α (u m+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(u m ) = S(u m+1 ) and the algorithm terminates or S(u m ) < S(u m+1 ).
3) The algorithm yields a strictly monotonically increasing sequence S(u m ), 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 be a symmetric third-order tensor and α ≥ 27 4 max i=1,2,...,n n j,k=1 (F ijk ) 2 . Let Ψ be an algorithm for the QAP which yields monotonic ascent w.r.t. the current iterate. Then the following holds for Algorithm 2: 1) The sequence F α (x k , y k , y k ) for k = 1, 2, . . . is strictly monotonically increasing or terminates.
2) The sequence of scores S(u m ) for m = 1, 2, . . . is strictly monotonically increasing or terminates. For every m, u m ∈ M is a valid assignment matrix.
3) The algorithm terminates after a finite number of iterations.

Proof:
The proof is similar to that of Theorem 3.5.
Proposition 3.8. Let F : R n × R n × R n → R be a symmetric third-order multilinear form and G : R n × R n × R n → 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 = y or y = z or z = x), the following inequality holds 3) For every α ∈ R, it holds arg max x∈M S α (x) = arg max x∈M S(x).

Proof:
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) ≤ max u∈{x,y,z} 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 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) y, z) .
**) Assume {x, y, z} ∈ M is a non-homogeneous tuple, we will prove for any 0 < ε < 1 that From Eq. (17), one has G ε (x, y, z) . The inequality becomes equality if and only if x = y = z. Applying this result for ε = 1/3 finishes the proof.

Proof:
Recall from Algorithm 1 the multilinear form G defined as G(x, y, z) = n i=1 ē i , x ē i , y ē i , z , and one has the relation F α (x, y, z) = F(x, y, z) + α G(x, y, z).
Indeed, it follows from Eq. (21) that G(u, u, u) > G(u, u, z) for all z ∈ M and z = u, which implies the maximum must be attained at z = u. Combining (22) and (23) 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. Table 5 shows the comparison between the proposed adaptive third-order methods with the state-of-the-art fourthorder approaches [36].

B.3 Running time
Tab. 6, Tab. 7 and Tab. 8 show the running time of all higherorder algorithms in seconds on the previous experiments. 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.