A Quadratic Penalty Method for Hypergraph Matching

TL;DR

This paper introduces a quadratic penalty method for hypergraph matching, reformulating it as a sparse tensor optimization problem and demonstrating effective support set recovery with theoretical guarantees and efficient algorithms.

Contribution

It presents a novel quadratic penalty approach for hypergraph matching, including theoretical analysis and practical algorithms for support set identification.

Findings

Supports exact support set recovery under certain conditions

The proposed method is efficient in accuracy and speed

Numerical results validate the effectiveness of the approach

Abstract

Hypergraph matching is a fundamental problem in computer vision. Mathematically speaking, it maximizes a polynomial objective function, subject to assignment constraints. In this paper, we reformulate the hypergraph matching problem as a sparse constrained tensor optimization problem. The optimality conditions are characterized based on the sparse constrained optimization theory. By dropping the sparsity constraint, we show that the resulting relaxation problem can recover the global minimizer of the original problem. The critical step in solving the original problem is to identify the location of nonzero entries (referred as support set) in a global minimizer. Inspired by such observations, we penalize the equality constraints and apply the quadratic penalty method to solve the relaxation problem. Under reasonable assumptions, we show that the support set of the global minimizer in hypergraph matching problem can be correctly identified when the number of iterations is sufficiently large. A projected gradient method is applied as a subsolver to solve the quadratic penalty subproblem. Numerical results demonstrate that the exact recovery of support set indeed happens, and the proposed algorithms are efficient in terms of both accuracy and speed.

Equations

$$\begin{array}[]{ll}\Pi^{1}=\{X\in\{0,1\}^{n_{1}\times n_{2}}:\sum\limits_{l_{2}=1}^{n_{2}}X_{l_{1}l_{2}}=1,\ l_{1}=1,\ldots,n_{1}\},\end{array}$$

$$\begin{array}[]{ll}\Pi^{2}=&\{X\in\{0,1\}^{n_{1}\times n_{1}}\,:\ \sum\limits_{l_{2}=1}^{n_{1}}X_{l_{1}l_{2}}=1,\ l_{1}=1,\ldots,n_{1};\\ &\sum\limits_{l_{1}=1}^{n_{1}}X_{l_{1}l_{2}}=1,\ l_{2}=1,\ldots,n_{1}\}.\end{array}$$

$$X_{l_{1}l_{2}}=\left\{\begin{array}[]{ll}1,&\hbox{ if }l_{1}\in V_{1}\hbox{ is assigned to }l_{2}\in V_{2};\\ 0,&\hbox{ otherwise.}\end{array}\right.$$

$$\max_{X\in\Pi^{1}}\ \ \sum_{\tiny\begin{array}[]{c}(l_{1},j_{1},k_{1})\in E_{1}\\ (l_{2},j_{2},k_{2})\in E_{2}\end{array}}\mathcal{B}_{l_{1}l_{2}j_{1}j_{2}k_{1}k_{2}}X_{l_{1}l_{2}}X_{j_{1}j_{2}}X_{k_{1}k_{2}}.$$

$$x:=(\bar{x}_{1}^{T},\ldots,\bar{x}_{n_{1}}^{T})^{T},\ \hbox{with }X=\left[\begin{array}[]{ccc}X_{11}&\cdots&X_{1n_{2}}\\ \vdots&\ddots&\vdots\\ X_{n_{1}1}&\cdots&X_{n_{1}n_{2}}\end{array}\right]:=\left[\begin{array}[]{c }\bar{x}_{1}^{T}\\ \vdots\\ \bar{x}_{n_{1}}^{T}\end{array}\right].$$

$${\mathcal{A}}_{ljk}={\mathcal{B}}_{l_1{l}_2{j}_1{j}_2{k}_1{k}_2},$$

$$l=(l_1-1)n_2+l_2,j=(j_1-1)n_2+j_2,k=(k_1-1)n_2+k_2.$$

$$\begin{array}[]{rl}\min\limits_{x\in\mathbb{R}^{n}}&\quad f(x):=-\frac{1}{6}\mathcal{A}x^{3}\\ \hbox{s.t.}&\quad e^{T}\bar{x}_{i}=1,\ i=1,\ldots,n_{1},\\ &\quad x\in\{0,1\},\end{array}$$

$${\mathcal{B}}_{l_1{l}_2{j}_1{j}_2{k}_1{k}_2}={\mathcal{B}}_{l_1^{\prime }{l}_2^{\prime }{j}_1^{\prime }{j}_2^{\prime }{k}_1^{\prime }{k}_2^{\prime }}.$$

$${\mathcal{A}}_{ljk}={\mathcal{A}}_{l^{\prime }{j}^{\prime }{k}^{\prime }}.$$

$$x_{-i}=({\bar{x}}_1^T,\dots ,{\bar{x}}_{i-1}^T,{\bar{x}}_{i+1}^T,\dots ,{\bar{x}}_{n_1}^T{)}^T,I(i,n_2)=\{(i-1)n_2+1,\dots ,i{n}_2\}.$$

$$f(x)$$

$$f^i({\bar{x}}_i,x_{-i})={\bar{x}}_i^T{\nabla }_{{\bar{x}}_i}f(x).$$

$$(\nabla f(x){)}_l=-\sum _{l<j<k}{\mathcal{A}}_{ljk}{x}_j{x}_k-\sum _{j<l<k}{\mathcal{A}}_{ljk}{x}_j{x}_k-\sum _{j<k<l}{\mathcal{A}}_{ljk}{x}_j{x}_k.$$

$$f^i({\bar{x}}_i,x_{-i})=-\sum _{l\in I(i,n_2)}\left(\sum _{l<j<k}{\mathcal{A}}_{ljk}{x}_l{x}_j{x}_k+\sum _{j<l<k}{\mathcal{A}}_{ljk}{x}_l{x}_j{x}_k+\sum _{j<k<l}{\mathcal{A}}_{ljk}{x}_l{x}_j{x}_k\right).$$

$$\begin{array}[]{rl}\min\limits_{x\in\mathbb{R}^{n}}&f(x)\\ \mathrm{s.t.}&e^{T}\bar{x}_{i}=1,\ i=1,\ldots,n_{1},\\ &x\geq 0,\quad\|x\|_{0}\leq n_{1}.\end{array}$$

$$\begin{array}[]{rl}\min\limits_{x\in\mathbb{R}^{n}}&f(x)\\ \mathrm{s.t.}&e^{T}\bar{x}_{i}=1,\ i=1,\ldots,n_{1},\\ &\hat{e}_{i}^{T}x=1,\ i=1,\ldots,n_{1},\\ &x\geq 0,\quad\|x\|_{0}\leq n_{1},\end{array}$$

$$\begin{array}[]{rl}\min\limits_{x\in\mathbb{R}^{n}}&f(x)\\ \mathrm{s.t.}&e^{T}\bar{x}_{i}=1,\ i=1,\ldots,n_{1},\ \ x\geq 0.\end{array}$$

$$\left\{\begin{array}[]{ll}\nabla_{\bar{x}_{i}}f(x)-\lambda_{i}e-\bar{\mu}_{i}=0,\\ \bar{x}_{i}\geq 0,\ \bar{\mu}_{i}\geq 0,\ \bar{x}_{i}^{T}\bar{\mu}_{i}=0,\\ e^{T}\bar{x}_{i}-1=0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\nabla_{\bar{x}_{i}}f(x)-\lambda_{i}e\geq 0,\ \bar{x}_{i}\geq 0,\ (\nabla_{\bar{x}_{i}}f(x)-\lambda_{i}e)^{T}\bar{x}_{i}=0,\\ e^{T}\bar{x}_{i}-1=0,\end{array}\right.$$

$${\mathcal{I}}_i(x)=\{p:({\bar{x}}_i{)}_p=0\},{\Gamma }_i(x)=\{p:({\bar{x}}_i{)}_p>0\}.$$

$$\left\{\begin{array}[]{lll}(\nabla_{\bar{x}_{i}}f(x))_{p}-\lambda_{i}=0,&\ (\bar{x}_{i})_{p}>0,&\ p\in\Gamma_{i}(x),\\ (\nabla_{\bar{x}_{i}}f(x))_{p}-\lambda_{i}\geq 0,&\ (\bar{x}_{i})_{p}=0,&\ p\in\mathcal{I}_{i}(x),\\ e^{T}\bar{x}_{i}-1=0,&&\end{array}\right.$$

$$f^i({\bar{x}}_i,x_{-i})$$

$$(\bar{y}^{1}_{i})_{p}=\left\{\begin{array}[]{ll}1,&\hbox{ if }p=p_{0};\\ 0,&\hbox{otherwise},\end{array}\right.\ \text{and}\quad\bar{y}^{1}_{i^{\prime}}=\left\{\begin{array}[]{ll}\bar{y}^{1}_{i},&\hbox{ if }i^{\prime}=i;\\ \bar{y}^{0}_{i^{\prime}},&\hbox{otherwise}.\end{array}\right.$$

$${\nabla }_{{\bar{x}}_i}f(y^0)={\nabla }_{{\bar{x}}_i}f(y^1).$$

$$f(y^1)-f(y^0)$$

$$\Vert y^r{\Vert }_0<\dots <\Vert y^1{\Vert }_0<\Vert y^0{\Vert }_0.$$

$$(\bar{x}^{*}_{i})_{p}=\left\{\begin{array}[]{ll}1,&\quad\hbox{ if }p=p_{i},\\ 0,&\quad\hbox{otherwise,}\end{array}\right.$$

$$\begin{array}[]{rl}\min\limits_{x\in\mathbb{R}^{n}}&\hat{f}(x)\\ \mathrm{s.t.}&e^{T}\bar{x}_{i}=\alpha_{i},\ i=1,\ldots,n_{1},\ \ x\geq 0,\end{array}$$

$$\begin{array}[]{rl}\min\limits_{x\in\mathbb{R}^{n}}&\hat{f}(x)\\ \mathrm{s.t.}&e^{T}\bar{x}_{i}=\alpha_{i},\ i=1,\ldots,n_{1},\ \ x\geq 0,\\ &\|x\|_{0}\leq n_{1}.\end{array}$$

Full text

1 Introduction

Recently, hypergraph matching has become a popular tool in establishing correspondence between two sets of points. It is a central problem in computer vision, and has been used to solve several applications, including object detection berg2005shape , image retrieval yan2005efficient , image stitching zaragoza13 ; zaragoza14 , and bioinformatics wu2012prl .

From the point of view of graph theory, hypergraph matching belongs to bipartite matching. Traditional graph matching models only use point-to-point features or pair-to-pair features, which can be solved by linear assignment algorithms jiang2007 ; maciel2003 or quadratic assignment algorithms egozi2013 ; Jiang2016 ; lee2011 ; litman2014 ; yan2015 , respectively. To use more geometric information such as angles, lines, and areas, triple-to-triple graph matching was proposed in 2008 zass2008probabilistic , and was further studied in duchenne2011tensor ; lee2011hyper ; nguyen2016efficient . Since three vertices are associated with one edge, it is also termed as hypergraph matching. Numerical experiments in literature duchenne2011tensor ; lee2011hyper ; nguyen2016efficient ; zass2008probabilistic show that hypergraph matching is more efficient than traditional graph matching. The aim of this paper is to study the hypergraph matching problem in both theory and algorithm.

The mathematical model of hypergraph matching is to maximize a multi-linear objective function subject to the row permutation constraints for $n_{1}\leq n_{2}$

[TABLE]

or permutation constraints for $n_{1}=n_{2}$

[TABLE]

We call a matrix satisfying (1.1) or (1.2) a binary assignment matrix. Optimization problems over binary assignment matrices are known to be NP-hard due to the combinatorial property.

Most existing algorithms for hypergraph matching relax the binary constraints into bound constraints and solve a continuous optimization problem. For instance, the probabilistic Hypergraph Matching method (HGM) zass2008probabilistic reformulated the constraints as the intersection of three convex sets, and successively projected the variables onto the sets until convergence. The Tensor Matching method (TM) duchenne2011tensor solved the optimization problem using the power iteration algorithm. The Hypergraph Matching method via Reweighted Random Walks (RRWHM) lee2011hyper dealt with the problem by walking among two feasible vectors randomly. Different from the above algorithms, Block Coordinate Ascent Graph Matching (BCAGM) nguyen2016efficient applied a block coordinate ascent framework, where they kept the binary constraints, and proposed to reformulate the multi-linear objective function into a linear one and solve it using linear assignment algorithms. All the existing algorithms require the equality constraints in (1.1) or (1.2) to be satisfied strictly at each iteration. In fact, we only expect that one of the elements is significantly larger than the others in each row or column of $X$. That is, the equality constraints are only soft constraints, which allow violations to some extent. Therefore, we penalize the equality constraint violations as part of the objective function in our algorithm.

The hypergraph matching problem can also be reformulated equivalently as a nonlinear optimization problem with sparse constraint. During the last few years, in the optimization community, there has been significant progress on solving sparse constrained nonlinear problems, particularly on dealing with optimality conditions and numerical algorithms in different situations. Recent development in optimality conditions can be found in Pan2017 , where based on decomposition properties of the normal cones, the authors characterized different kinds of stationary points and performed detailed investigations on relations of local minimizers, global minimizers and several types of stationary points. Other related work includes Bauschke2014Restricted ; Burdakov2015MATHEMATICAL ; C2016Constraint ; Li2015 ; Pan2015On . The related algorithms can be summarized into two approaches. One is the direct approach, aiming at dealing with the sparse constraint directly, such as the hard-thresholding type based algorithms Beck2012Sparsity ; Pan2016A and the $\ell_{0}$ penalty based algorithms Lu2012Sparse . The other one is the relaxation approach such as the $\ell_{p}$ regularization based algorithms ChenLeiLuYe ; Jiang2016 . In particular, an efficient $\ell_{p}$ regularization algorithm was proposed in Jiang2016 , which deals with problems over the permutation matrix constraints (1.2). It can be applied to solve the hypergraph matching problem subject to (1.2).

Motivation. Noting that hypergraph matching is essentially a mixed integer programming, most existing methods relax the integer constraints as box constraints, and solve the relaxed continuous optimization problem. A natural question is: what is the relation between hypergraph matching and the relaxation problem? Furthermore, the key step in solving this problem is actually to identify the support set of the global minimizer. None of the existing algorithms has taken this fact into account. This leads to the second question: can we make use of this insight to design our algorithm?

Our Contributions. In this paper, by reformulating hypergraph matching equivalently as a sparse constrained optimization problem, we study it from the following aspects.

• •

Relaxation problem. By dropping the sparse constraint, we show that the relaxation problem can recover the solution of the original problem in the sense that the former problem shares at least one global minimizer with the latter one (Theorem 3.1). This result highly depends on the special structures of hypergraph matching. Furthermore, we show that Theorem 3.1 can be extended to more general problems (Corollary 2). For any global minimizer of the relaxation problem, we propose a procedure to reduce its sparsity until a global minimizer of the original problem is reached.

• •

Quadratic penalty method. Our aim is to identify the support set of a global minimizer of the original problem, thus the equality constraints are not necessary to be satisfied strictly. This motivates us to penalize the equality constraint violations, and solve the relaxation problem by a quadratic penalty method. We show that under reasonable assumptions, the support set of a global minimizer of the original problem can be recovered exactly, when the number of iteration is sufficiently large (Theorems 4.2 and 4.3).

• •

Projected gradient method. For the quadratic penalty subproblem, which is a nonlinear problem with simple box constraints, we choose one of the active set based methods called the projected gradient method as a subsolver. The advantage of the active set based method is that it well fits our motivation, which is to identify the support set of the solution rather than to look for the magnitude. Numerical results demonstrate that the exact recovery of the support set indeed happens, and the proposed algorithm is particularly suitable for large-scale problems.

Organization. The rest of the paper is organized as follows. In Section 2, we introduce the reformulation of the hypergraph matching problem, and discuss several preliminary properties. In Section 3, we study the properties of the relaxation problem by dropping the sparse constraint. In Section 4, we study the quadratic penalty method by penalizing the equality constraint violations and establish the convergence results in terms of support set under different situations. An existing projected gradient method is also discussed to solve the quadratic penalty subproblem. Numerical experiments are reported in Section 5. Final conclusions are drawn in Section 6.

Notations. For $x\in\mathbb{R}^{n}$, define the active set as $\mathcal{I}(x)=\{l:\,x_{l}=0\}$ and the support set as $\Gamma(x)=\{l:\,x_{l}>0\}$. We also use $\mathcal{I}^{k}$ and $\Gamma^{k}$, and $\mathcal{I}^{*}$ and $\Gamma^{*}$ to denote the corresponding sets at $x^{k}$ and $x^{*}$, respectively. Let $|\mathcal{I}|$ be the number of elements in the set $\mathcal{I}$. $\|x\|$ denotes the $\ell_{2}$ norm of $x$, $\|x\|_{0}$ the number of nonzero entries in $x$, and $\|x\|_{\infty}$ the infinity norm of $x$.

2 Problem Reformulation

In this section, we will reformulate hypergraph matching as a sparse constrained optimization problem, and discuss several preliminary properties.

2.1 Hypergraph matching problem

In this part, we will give the mathematical formulation for hypergraph matching, including its objective function and constraints.

Consider two hypergraphs $G_{1}=\{V_{1},E_{1}\}$, and $G_{2}=\{V_{2},E_{2}\}$, where $V_{1}$ and $V_{2}$ are sets of points with $|V_{1}|=n_{1}$, $|V_{2}|=n_{2}$, and $E_{1}$, $E_{2}$ are sets of hyperedges. In this paper, we always suppose that $n_{1}\leq n_{2}$, and each point in $V_{1}$ is matched to exactly one point in $V_{2}$, while each point in $V_{2}$ can be matched to arbitrary number of points in $V_{1}$. That is, we focus on (1.1). For each hypergraph, we consider three-uniform hyperedges. Namely, the three points involved in each hyperedge are different, for example, $(l_{1},j_{1},k_{1})\in E_{1}$. Our aim is to find the best correspondence (also referred to as ‘matching’) between $V_{1}$ and $V_{2}$ with the maximum matching score.

Let $X\in\mathbb{R}^{n_{1}\times n_{2}}$ be the assignment matrix between $V_{1}$ and $V_{2}$, i.e.,

[TABLE]

Two hyperedges $(l_{1},j_{1},k_{1})\in E_{1}$ and $(l_{2},j_{2},k_{2})\in E_{2}$ are said to be matched if $l_{1},j_{1},k_{1}\in V_{1}$ are assigned to $l_{2},j_{2},k_{2}\in V_{2}$, respectively. It can be represented equivalently by $X_{l_{1}l_{2}}X_{j_{1}j_{2}}X_{k_{1}k_{2}}=1$. Let $\mathcal{B}_{l_{1}l_{2}j_{1}j_{2}k_{1}k_{2}}$ be the matching score between $(l_{1},j_{1},k_{1})$ and $(l_{2},j_{2},k_{2})$. Then $\mathcal{B}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{1}\times n_{2}\times n_{1}\times n_{2}}$ is a sixth order tensor. Assume $\mathcal{B}$ is given, satisfying $\mathcal{B}_{l_{1}l_{2}j_{1}j_{2}k_{1}k_{2}}\geq 0$ if $(l_{1},j_{1},k_{1})\in E_{1}$ and $(l_{2},j_{2},k_{2})\in E_{2}$, and $\mathcal{B}_{l_{1}l_{2}j_{1}j_{2}k_{1}k_{2}}=0$, otherwise.

Given hypergraphs $G_{1}=\{V_{1},E_{1}\}$, $G_{2}=\{V_{2},E_{2}\}$, and the matching score $\mathcal{B}$, the hypergraph matching problem takes the following form

[TABLE]

Note that (2.1) is a matrix optimization problem, which can be reformulated as a vector optimization problem as follows.

Let $n=n_{1}n_{2}$, $x\in\mathbb{R}^{n}$ be the vectorization of $X$, that is

[TABLE]

Here, $\bar{x}_{i}\in\mathbb{R}^{n_{2}}$ is the $i$-th block of $x$. In the following, for any vector $z\in\mathbb{R}^{n}$, we always assume it has the same partition as $x$. Define $\mathcal{A}\in\mathbb{R}^{n\times n\times n}$ as

[TABLE]

where

[TABLE]

Consequently, (2.1) can be reformulated as

[TABLE]

where $e\in\mathbb{R}^{n_{2}}$ is a vector with all entries equal to one, and $\mathcal{A}x^{3}:=\sum\limits_{l,j,k}\mathcal{A}_{ljk}x_{l}x_{j}x_{k}$.

2.2 Preliminary properties

In this subsection, we will discuss several properties of $\mathcal{A}$, $\mathcal{B}$, and $f(x)$. We begin with properties of $\mathcal{B}$.

Proposition 1

• (i)

$\mathcal{B}_{l_{1}l_{2}j_{1}j_{2}k_{1}k_{2}}\geq 0$* for all $l_{1},j_{1},k_{1}=1,\ldots,n_{1}$ and $l_{2},j_{2},k_{2}=1,\ldots,n_{2}$;*

• (ii)

If $(l_{1},j_{1},k_{1})\in E_{1}$, then $l_{1}$, $j_{1}$, and $k_{1}$ are distinct. If $(l_{2},j_{2},k_{2})\in E_{2}$, then $l_{2}$, $j_{2}$, and $k_{2}$ are also distinct;

• (iii)

For any permutation operator $\pi$, suppose $\pi(l_{1},j_{1},k_{1})=(l_{1}^{\prime},j_{1}^{\prime},k_{1}^{\prime})$ and $\pi(l_{2},j_{2},k_{2})=(l_{2}^{\prime},j_{2}^{\prime},k_{2}^{\prime})$. There is

[TABLE]

The above properties of $\mathcal{B}$ result in the following properties of $\mathcal{A}$ directly.

Proposition 2

• (i)

$\mathcal{A}_{ljk}\geq 0$, for all $l,j,k=1,\ldots,n$;

• (ii)

For nonzero entries of $\mathcal{A}$, say $\mathcal{A}_{ljk}$, $x_{l}$, $x_{j}$ and $x_{k}$ come from different blocks of $x$;

• (iii)

Suppose $(l^{\prime},j^{\prime},k^{\prime})$ is any permutation of $(l,j,k)$. Then

[TABLE]

In other words, $\mathcal{A}$ is nonnegative and symmetric.

Proof. (i) follows directly from the nonnegativity of $\mathcal{B}$. In terms of (ii), by the definition of $\mathcal{A}$, there exist $(l_{1},j_{1},k_{1})$ and $(l_{2},j_{2},k_{2})$ such that (2.2) and (2.3) hold. Further, we know that $x_{l}$ is the $l_{2}$-th entry in the $l_{1}$-th block of $x$, i.e., $x_{l}=(\bar{x}_{l_{1}})_{l_{2}}$. Similarly, $x_{j}=(\bar{x}_{j_{1}})_{j_{2}}$ and $x_{k}=(\bar{x}_{k_{1}})_{k_{2}}$. By (ii) in Proposition 1, $l_{1},j_{1},k_{1}$ are distinct, which implies that $x_{l}$, $x_{j}$, $x_{k}$ come from different blocks of $x$. In terms of (iii), since $\pi(l_{1},j_{1},k_{1})=(l_{1}^{\prime},j_{1}^{\prime},k_{1}^{\prime})$, $\pi(l_{2},j_{2},k_{2})=(l_{2}^{\prime},j_{2}^{\prime},k_{2}^{\prime})$, again by the definition of $\mathcal{A}$ and (2.3), there is $B_{l_{1}^{\prime}l_{2}^{\prime}j_{1}^{\prime}j_{2}^{\prime}k_{1}^{\prime}k_{2}^{\prime}}=A_{l^{\prime}j^{\prime}k^{\prime}}$. Together with (2.5) and (2.2), there is (2.6). ∎

Different from other nonlinear problems, the homogenous polynomial $f(x)$ enjoys special structures. To see this, for the $i$-th block $\bar{x}_{i}$, denote

[TABLE]

Rewrite $f(x)$ as follows:

[TABLE]

Proposition 3

• (i)

For each block $\bar{x}_{i}$, $i\in\{1,\ldots,n_{1}\}$, $f(x)$ is a linear function of $\bar{x}_{i}$, i.e., $\nabla_{\bar{x}_{i}}f(x)$ is independent of $\bar{x}_{i}$;

• (ii)

[TABLE]

Proof. In terms of (i), by the definition of $\mathcal{A}$, we only need to consider the term $\mathcal{A}_{ljk}x_{l}x_{j}x_{k}$, where $\mathcal{A}_{ljk}$ is nonzero. Due to (ii) in Proposition 2, $\mathcal{A}_{ljk}x_{l}x_{j}x_{k}$ is linear in each related block $\bar{x}_{l_{1}}$, $\bar{x}_{j_{1}}$, and $\bar{x}_{k_{1}}$. Therefore, $f(x)$ is a linear function of $\bar{x}_{i}$, $i=1,\ldots,n_{1}$.

In terms of (ii), the elements of gradient $\nabla f(x)$ take the following form

[TABLE]

Rewrite $f^{i}(\bar{x}_{i},x_{-i})$ in (2.7) as

[TABLE]

Hence, $f^{i}(\bar{x}_{i},x_{-i})=\sum\limits_{l\in I(i,n_{2})}x_{l}(\nabla f(x))_{l}$, which gives (2.8). ∎

Equation (2.8) will be useful in Section 3.

2.3 Sparse constrained optimization problem

Problem (2.4) is a 0-1 mixed integer programming, which is one of Karp’s 21 NP-complete problems Karp1972 . In this subsection, we will reformulate (2.4) into a sparse constrained optimization problem.

By direct computations, (2.4) can be reformulated as the following sparse constrained minimization problem

[TABLE]

To see this, for each $x$ satisfying the equality constraints, we have $\|x\|_{0}\geq n_{1}$. Together with $\|x\|_{0}\leq n_{1}$, we actually have $\|x\|_{0}=n_{1}$.

In particular, if $n_{1}=n_{2}$, by the permutation constraints (1.2), problem (2.9) reduces to the following hypergraph matching problem

[TABLE]

where $\hat{e}_{i}=((e^{n_{1}}_{i})^{T},\ldots,(e^{n_{1}}_{i})^{T})^{T}\in\mathbb{R}^{n},$ and $e^{n_{1}}_{i}$ is the $i$-th column of the $n_{1}$-by-$n_{1}$ identity matrix.

Remark 1

Note that the dimension of $x$ is $n=n_{1}n_{2}$, which can be large even for moderate $n_{1}$ and $n_{2}$. For instance, if $n_{1}=100$ and $n_{2}=100$, then $n=10^{4}$, and the number of elements in $\mathcal{A}$ will be around $10^{12}$. Hence, algorithms capable of dealing with large-scale problems are highly in demand.

Remark 2

Problem (2.9) is essentially a [math]-$1$ mixed integer programming. Each feasible point is actually an isolated feasible point, which means that it is a strict local minimizer and of course is a stationary point of (2.9). For a theoretical verification from the optimality point of view, please see Theorems 1 and 3 in an earlier version of our paper CuiLiQiYan2017 .

3 Relaxation Problem of (2.9)

In this section, we will study the relaxation problem (3.1) and its connections with the original problem (2.9).

By dropping the sparse constraint in (2.9), we obtain the following problem (referred to as the relaxation problem)

[TABLE]

As we will show later in Theorem 3.1, although we drop the sparse constraint, the relaxation problem (3.1) still admits a global minimizer with sparsity $n_{1}$ due to the special structures of (2.9). That is, the relaxation problem (3.1) recovers a global minimizer of (2.9).

Let $\lambda\in\mathbb{R}^{n_{1}},\mu\in\mathbb{R}^{n}$ be the Lagrange multipliers corresponding to the constraints in (3.1). The KKT conditions of (3.1) are

[TABLE]

which are equivalent to

[TABLE]

for all $i=1,\ldots,n_{1}$. Define the active set and the support set for the $i$-th block $\bar{x}_{i}$ as

[TABLE]

The KKT conditions can be reformulated as

[TABLE]

for all $i=1,\ldots,n_{1}$. The above analysis gives the following lemma.

Lemma 1

Let $x\in\mathbb{R}^{n}$ be a stationary point of (3.1), and $\lambda\in\mathbb{R}^{n_{1}}$ be the Lagrange multiplier corresponding to the equality constraints. For all $i=1,\ldots,n_{1}$, we have

• (i)

$\lambda_{i}=\min\limits_{p\in\{1,,\ldots,n_{2}\}}(\nabla_{\bar{x}_{i}}f(x))_{p},\hbox{ and }(\nabla_{\bar{x}_{i}}f(x))_{p}=\lambda_{i}$* for $p\in\Gamma_{i}(x)$;*

• (ii)

$f^{i}(\bar{x}_{i},x_{-i})=\lambda_{i}$.

Proof. (i) can be obtained directly from the KKT conditions (3.3).

In terms of (ii), by (2.8), there is

[TABLE]

where the last equality is due to $e^{T}\bar{x}_{i}=1$. This completes the proof.∎

Theorem 3.1

There exists a global minimizer $x^{*}$ of (3.1) such that $\|x^{*}\|_{0}=n_{1}$. Furthermore, $x^{*}$ is a global minimizer of (2.9).

Proof. Without loss of generality, let $y^{0}$ be a global minimizer of (3.1) with $\|y^{0}\|_{0}>n_{1}$. Find the first block of $y^{0}$, denoted as $\bar{y}^{0}_{i}$, such that $\|\bar{y}^{0}_{i}\|_{0}>1$. Now we choose one index $p_{0}$ from $\Gamma_{i}(y^{0}):=\{p:(\bar{y}^{0}_{i})_{p}>0\}$, and define a new point $y^{1}=((\bar{y}^{1}_{1})^{T},\ldots,(\bar{y}^{1}_{n_{1}})^{T})^{T}$ as follows:

[TABLE]

Then $y^{1}$ is a feasible point for (3.1), and satisfies $y^{1}_{-i}=y^{0}_{-i}$. Furthermore, by Proposition 3, $\nabla_{\bar{x}_{i}}f(x)$ is a function of $x_{-i}$, there is

[TABLE]

Next, we will show that $f(y^{1})=f(y^{0})$. Indeed,

[TABLE]

This gives that $y^{1}$ is a feasible point with $f(y^{1})=f(y^{0})$. In other words, $y^{1}$ is another global minimizer of (3.1) with $\|y^{1}\|_{0}<\|y^{0}\|_{0}$. If $\|y^{1}\|_{0}=n_{1}$, let $x^{*}:=y^{1}$. Otherwise, by repeating the above process, we can obtain a finite sequence $y^{0},y^{1},\ldots,y^{r}$, which are all feasible points for (3.1) satisfying

[TABLE]

Note that there are $n_{1}$ blocks in $y^{0}\in\mathbb{R}^{n}$. After at most $n_{1}$ steps, the process will stop. In other words, $1\leq r\leq n_{1}$. The final point $y^{r}$ will satisfy $\|y^{r}\|_{0}=n_{1}$. One can obtain a global minimizer $x^{*}:=y^{r}$ of (3.1) with $n_{1}$ nonzero elements.

Next, we will show that $x^{*}$ is also a global minimizer of (2.9). Note that the feasible region of (2.9) is a subset of the feasible region of (3.1). $\|x^{*}\|_{0}=n_{1}$ implies that $x^{*}$ is also a feasible point for (2.9). Together with the fact that $f(x^{*})$ attains the global minimum of (3.1), we conclude that $x^{*}$ is a global minimizer of (2.9). ∎

Theorem 3.1 shows that $\|x^{*}\|_{0}=n_{1}$ is a necessary and sufficient condition for a global minimizer $x^{*}$ of (3.1) to be a global minimizer of (2.9). We highlight this relation in the following corollary.

Corollary 1

A global minimizer $x^{*}$ of (3.1) is a global minimizer of (2.9) if and only if $\|x^{*}\|_{0}=n_{1}$.

A special case of Theorem 3.1 is $|\Gamma_{i}(x^{*})|=1$, for each $i\in\{1,\ldots,n_{1}\}$. Then the global minimizer $x^{*}$ of (3.1) is a global minimizer of (2.9).

Remark 3

From the proof of Theorem 3.1, one can start from any global minimizer $y^{0}$ of (3.1) to reach a point $x^{*}$, which is a global minimizer of both (2.9) and (3.1). We only need to choose one index as the location of nonzero entry in each block $\bar{y}_{i}^{0}$. Assume $p_{i}$ is chosen from $\Gamma_{i}(y^{0})$. Let $\Gamma^{*}_{i}=p_{i}$. This will give the support set in the $i$-th block, which in turn determines the global minimizer $x^{*}$ of (2.9) by

[TABLE]

for each $p\in\{1,\ldots,n_{2}\}$ and $i\in\{1,\ldots,n_{1}\}$. One particular method to choose $p_{i}$ is to choose the index with the largest value within the block. This is actually the projection of $y^{0}$ onto the feasible set of (2.9). Here, we summarize the process in Algorithm 1.

Note that HGM zass2008probabilistic also solves the relaxation problem (3.1), whereas TM duchenne2011tensor and RRWHM lee2011hyper solve the relaxation problem with the permutation constraints (1.2). However, none of them analyzes the connections between the original problem and the relaxation problem in terms of global minimizers. On contrast, the result in Theorem 3.1 reveals for the first time the connections between the original problem (2.9) and the relaxation problem (3.1), which is one of the main differences of our work from existing algorithms for hypergraph matching.

Theorem 3.1 reveals an interesting connection between the original problem (2.9) and the relaxation problem (3.1) in terms of global minimizers. The result heavily relies on the property of $f(x)$ in Proposition 3, as well as the equality constraints in (2.9). It can be extended to the following general case.

Corollary 2

Consider

[TABLE]

where $\alpha_{i}>0$, and $\bar{x}_{i}\in\mathbb{R}^{m_{i}}$ with $m_{i}$ being positive integers satisfying $\sum_{i=1}^{n_{1}}m_{i}=n$. Suppose that $\hat{f}(x)$ satisfies Proposition 3. Then there exists a global minimizer $x^{*}$ of (3.5) such that $\|x^{*}\|_{0}=n_{1}$. Furthermore, $x^{*}$ is a global minimizer of the following problem

[TABLE]

4 The Quadratic Penalty Method

In this section, we will consider the quadratic penalty method for the relaxation problem (3.1). It contains three parts. The first part is devoted to motivating the quadratic penalty problem and its preliminary properties. The second part mainly focuses on the quadratic penalty method and the convergence in terms of the support set. In the last part, we apply an existing projected gradient method for the quadratic penalty subproblem.

4.1 The quadratic penalty problem

Note that (3.1) is a nonlinear problem with separated simplex constraints, which can be solved by many traditional nonlinear optimization solvers such as fmincon in MATLAB. As mentioned in Section 1, existing algorithms for hypergraph matching require the equality constraints in (3.1) to be satisfied strictly. On contrast, our aim here is actually to identify the support set of a global minimizer of (3.1) rather than the magnitude. Once the support set is found, we can follow the method in Remark 3 to obtain a global minimizer of (2.9). Inspired by such observations, we penalize the equality constraint violations as part of the objective function. This is another main difference of our method from existing algorithms. It leads us to the following quadratic penalty problem

[TABLE]

where $\sigma>0$ is a penalty parameter. However, this problem is not well defined in general, since for a fixed $\sigma$ the global minimizer will approach infinity. We can add an upper bound to make the feasible set bounded. This gives the following problem

[TABLE]

where $M\geq 1$ is a given number. (4.1) is actually the quadratic penalty problem of the following problem

[TABLE]

which is equivalent to (3.1).

Having introduced the quadratic penalty problem (4.1), next we will analyze the properties of (4.1) and its connection with the relaxation problem (3.1).

The Lagrangian function of (4.1) is

[TABLE]

where $w$ and $\nu$ are the Lagrange multipliers corresponding to the inequality constraints in (4.1). The KKT conditions are

[TABLE]

for each $i\in\{1,\ldots,n_{1}\}$. In particular, for a stationary point $x$ of (4.1), let $\mathcal{I}_{i}(x)$ and $\Gamma_{i}(x)$ be defined by (3.2). Define

[TABLE]

The KKT conditions are equivalent to the following, for each $i\in\{1,\ldots,n_{1}\}$,

[TABLE]

Define the violations of the equality constraints $h\in\mathbb{R}^{n_{1}}$ as

[TABLE]

There is

[TABLE]

The above analysis can be stated in the following lemma.

Lemma 2

Let $x\in\mathbb{R}^{n}$ be a stationary point of (4.1). We have $h_{i}\geq 0$ for all $i=1,\ldots,n_{1}$.

Proof. For each $i$, consider two cases. If $\mathcal{I}_{i}(x)\cup{\widehat{\Gamma}}_{i}(x)\neq\emptyset$, by (4.2), there exists $p\in\mathcal{I}_{i}(x)\cup{\widehat{\Gamma}}_{i}(x)$ such that $\sigma h_{i}\geq-(\nabla_{\bar{x}_{i}}f(x))_{p}$. By the nonnegativity of the entries in $\mathcal{A}$ and $x$, there is $-\nabla_{\bar{x}_{i}}f(x)\geq 0$ and $h_{i}\geq 0$. If $\mathcal{I}_{i}(x)\cup{\widehat{\Gamma}}_{i}(x)=\emptyset$, then $|\overline{\Gamma}_{i}(x)|=n_{2}$. In other words, $(\bar{x}_{i})_{p}=M$ for all $p=1,\ldots,n_{2}$. Then $h_{i}=e^{T}\bar{x}_{i}-1=n_{2}M-1\geq 0$. ∎

Let $u\in\mathbb{R}^{n}$ and $c\in\mathbb{R}^{n_{1}}$ be defined by

[TABLE]

and

[TABLE]

where $\bar{u}_{i}$ is the $i$-th block of $u$. It follows from the nonnegativity of $\mathcal{A}$ that $c\geq 0$. The following lemma describes the relation between the penalty parameter $\sigma$ and the violations of the equality constraints.

Lemma 3

For each stationary point $x$ of (4.1), there is

[TABLE]

where $h_{i}$ is defined by (4.3), and $c_{i}$ is defined by (4.5).

Proof. Note that $x\in[0,M]$. By the definition $f(x)=-\frac{1}{6}\mathcal{A}x^{3}$, we have $-(\nabla f(x))_{l}\leq M^{2}u_{l}.$ Together with (4.4), there is \sigma h_{i}$$\leq M^{2}\max_{p\in\{1,\ldots,n_{2}\}}(\bar{u}_{i})_{p} $=c_{i}$. The proof is complete. ∎

Lemma 4

For each feasible point $x\in\mathbb{R}^{n}$ of (3.1), it is a stationary point of (4.1) if and only if for all $i=1,\ldots,n_{1}$, there is

[TABLE]

Proof. Let $x$ be a feasible point for (3.1). There is $e^{T}\bar{x}_{i}-1=0$, $i=1,\ldots,n_{1}$. If $x$ is a stationary point of (3.1), by the KKT conditions (4.2), we have

[TABLE]

Consequently, $(\nabla_{\bar{x}_{i}}f(x))_{p}\geq-\sigma(e^{T}\bar{x}_{i}-1)=0$. On the other hand, $(\nabla_{\bar{x}_{i}}f(x))_{p}\leq 0$ due to the nonnegativity of entries in $\mathcal{A}$ and $x$. Therefore, $(\nabla_{\bar{x}_{i}}f(x))_{p}=0$ for all $p\in\mathcal{I}_{i}(x)$. For $p\in\Gamma_{i}(x)$, there is $(\nabla_{\bar{x}_{i}}f(x))_{p}=0.$ This gives (4.6).

Conversely, for a feasible point $x$ for (2.9), if (4.6) holds, the first two conditions in (4.2) hold by $\bar{x}_{i}^{T}e-1=0$, $i=1,\ldots,n_{1}$. For the third condition in (4.2), consider two cases. If $\overline{\Gamma}(x)=\emptyset,$ the result is trivial. Otherwise, there is $(\nabla_{\bar{x}_{i}}f(x))_{p}\leq 0$ due to the nonnegativity of entries in $\mathcal{A}$ and $x$. The third condition holds automatically. In both two cases, $x$ satisfies (4.2). That is, $x$ is a stationary point of (4.1). ∎

4.2 A quadratic penalty method for (3.1)

Having investigated the properties of the quadratic penalty problem, we then solve (3.1) by the traditional quadratic penalty method, i.e., by solving (4.1) sequentially. At each iteration, $x^{k}$ is a global minimizer of the following problem

[TABLE]

The quadratic penalty method is given in Algorithm 2.

The following theorem addresses the convergence of the quadratic penalty method, which can be found in classic optimization books such as (NumericalOpt, , Theorem 17.1) and (SunYuan, , Corollary 10.2.6). Therefore, the proof is omitted.

Theorem 4.1

Let $\{x^{k}\}$ be generated by Algorithm 2, and $\lim_{k\to\infty}\sigma_{k}=+\infty$. Then any accumulation point of the generated sequence $\{x^{k}\}$ is a global minimizer of (3.1).

Due to Theorem 4.1, in following analysis, we always assume the following holds.

Assumption 1

Let $\{x^{k}\}$ be generated by Algorithm 2 and $\lim_{k\rightarrow\infty}\sigma_{k}=+\infty$. Denote $K$ as a subset of $\{1,2,\ldots\}$. Assume that $\lim_{k\rightarrow\infty,k\in K}x^{k}=z$, and $z$ is a global minimizer of (3.1).

The next theorem mainly addresses the relation between the support set of $x^{k}$ and that of the global minimizer of (2.9). Recall that for $x^{k}$, there is

[TABLE]

Theorem 4.2

Suppose that Assumption 1 holds. If there exists a positive integer $k_{0}$, such that $\|x^{k}\|_{0}=n_{1}$ for all $k\geq k_{0},\ k\in K$, then there is a positive integer $k_{1}\geq k_{0}$ such that the support set of $z$ can be identified correctly. That is,

[TABLE]

Furthermore, $z$ is a global minimizer of (2.9).

Proof. First, we show $|\Gamma(z)|=n_{1}$. Noting that $z$ is a global minimizer of (4.1), we have

[TABLE]

Since $\lim_{k\to+\infty,k\in K}x^{k}=z$, there exists a positive integer $k^{\prime}$ such that for $k\geq k^{\prime},\ k\in K$, there is

[TABLE]

This implies that $\Gamma(z)\subseteq\Gamma^{k}$. It follows from the assumption that $|\Gamma^{k}|=n_{1}$ for all $k\geq k_{0},\ k\in K$. Consequently, we have $|\Gamma(z)|=n_{1}$. Therefore, $\Gamma(z)=\Gamma^{k}$ holds for $k\geq k_{1}:=\max\{k_{0},k^{\prime}\}$. The second part holds following the second part of Theorem 3.1. The proof is finished.∎

Theorem 4.2 indicates that we do not need to drive $\sigma_{k}$ to infinity since only the support set of $z$ is needed. If the conditions in Theorem 4.2 hold, then we can stop the algorithm when the number of elements in $\mathcal{I}^{k}$ keeps unchanged for several iterations. However, if there is $\|x^{k}\|_{0}>n_{1}$, we need more notations to analyze the connections.

Let $\mathcal{J}^{k}_{i}$ be the set of indices corresponding to the largest values in the $i$-th block $\bar{x}^{k}_{i}$, $p^{k}_{i}$ be the smallest index in $\mathcal{J}^{k}_{i}$, and $\mathcal{J}^{k}$ be the set of indices containing the largest values in each block of $x^{k}$, i.e.,

[TABLE]

Similarly, we define

[TABLE]

Theorem 4.3

Suppose that Assumption 1 holds.

• (i)

If $\|z\|_{0}=n_{1}$, then there exists an integer $k_{0}>0$, such that $\Gamma(z)=\mathcal{J}^{k}$ for all $k\geq k_{0},\ k\in K$;

• (ii)

If $\|z\|_{0}>n_{1}$ and $|\mathcal{J}_{i}(z)|=1$ for all $i=1,\ldots,n_{1}$, then there exists a global minimizer $x^{*}$ of (2.9) and a positive integer $k_{0}$, such that for all $k\geq k_{0},k\in K$, there is $\Gamma^{*}=\mathcal{J}^{k}$;

• (iii)

If $\|z\|_{0}>n_{1}$ and $|\mathcal{J}_{i}(z)|>1$ for one $i=1,\ldots,n_{1}$, then there exists a global minimizer $x^{*}$ of (2.9), a subsequence $\{x^{k}\}_{k\in K^{\prime}}$ and a positive integer $k_{0}$, such that for all $k\geq k_{0},k\in K^{\prime}$, there is $\Gamma^{*}=\mathcal{J}^{k}$.

Proof. With Theorem 3.1 and $\|z\|_{0}=n_{1}$, $z$ must be a global minimizer of (2.9). By the definition of $\Gamma(z)$ and $\mathcal{I}(z)$, there exists an integer $k_{0}>0$, such that for all $k\geq k_{0},k\in K$, there is $z_{l}>z_{l^{\prime}}$ for ${l^{\prime}}\in\mathcal{I}(z)$ and $l\in\Gamma(z)$. This gives $\mathcal{J}^{k}_{i}=\mathcal{J}_{i}(z)$ and (i).

In terms of (ii), $|\mathcal{J}_{i}(z)|=1$ implies that for $k\in K$ sufficiently large, there is

[TABLE]

Consequently, there is $\mathcal{J}^{k}_{i}=\mathcal{J}_{i}(z)$. Now let $y^{0}:=z$. Similar to the arguments in the proof of Theorem 3.1, we construct $y^{1}$ by choosing $p_{0}=\mathcal{J}_{i}(z)$. Then we can obtain a finite sequence $y^{0},y^{1},\ldots,y^{r}$ with

[TABLE]

After at most $n_{1}$ steps, the process will stop. In other words, $1\leq r\leq n_{1}$. At the final point $y^{r}$ will satisfy $\|y^{r}\|_{0}=n_{1}$. One can find a global minimizer $x^{*}:=y^{r}$ of problem (3.1) with sparsity $n_{1}$. Further, $x^{*}$ is also a global minimizer of (2.9) and satisfies

[TABLE]

Consequently, (ii) holds.

For (iii), suppose there exists an index $q_{1}$ such that $|\mathcal{J}_{q_{1}}(z)|>1$. Consequently, there exists $p_{1}\in\mathcal{J}^{k}_{q_{1}}$, such that for $k\in K$ sufficiently large, there are infinite number of $k$ satisfying $\mathcal{J}_{q_{1}}^{k}=p_{1}$. Denote the corresponding subsequence as $\{x^{k}\}_{k\in K_{1}}$, where $K_{1}\subset K$. Similarly, for $|\mathcal{J}_{q_{2}}(z)|>1$, we can find an infinite number of $k\in K_{2}\subseteq K_{1}$ such that $\mathcal{J}_{q_{2}}^{k}=p_{2}$. Repeating the process until for all blocks, there exists an integer $k_{0}>0$, such that $|\mathcal{J}_{i}^{k}|=1$, $i=1,\ldots,n_{1}$, for all $k\in K_{t}\subseteq K_{t-1}\ldots\subseteq K_{1}$, $k\geq k_{0}$. Let $K^{\prime}:=K_{t}$. Now similar to Remark 3, for all $i=1,\ldots,n_{1}$, we define $x^{*}$ as follows:

[TABLE]

Then we find a global minimizer of (3.1) such that $\|x^{*}\|_{0}=n_{1}$. For $k\geq k_{0}$, $k\in K^{\prime}$, there is $\mathcal{J}_{i}^{k}=\mathcal{J}_{i}^{*}$, $i=1,\ldots,n_{1}$. Consequently, $x^{*}$ is also a global minimizer of (2.9). Hence, (iii) holds. This completes the proof. ∎

Theorems 4.2 and 4.3 state that there is always a subsequence of $x^{k}$ whose support set will coincide with the support set of one global minimizer of (2.9). Consequently, it provides a method to design the termination rule for Algorithm 2.

4.3 A projected gradient method for the subproblem (4.1)

In this subsection, we will use a projected gradient method to solve the subproblem.

Note that the subproblem (4.1) is a nonlinear problem with simple box constraints. Various methods can be chosen to solve (4.1), one of which is the active set based method. We prefer such type of method because it quite fits our motivation to identify the support set of the global minimizer of (2.9) rather than the magnitude. The strategy of identifying the active set is therefore crucial in solving (4.1). We choose a popular approach proposed in Bertsekas1982Bertsekas , and modify it into the resulting projected gradient method, as shown in Algorithm 3. Other typical projected gradient methods in calamai1987 ; dai2006 can also be used.

Remark 4

Note that the projected gradient method is only guaranteed to converge to a stationary point. Based on Lemma 2, the sum of each block in the stationary point is larger than or equal to one. In other words, at least one entry in each block is larger than zero. This will partly explain the numerical observation that the magnitudes of the returned solution by our algorithm clearly fall into two parts: the estimated active part, which is close to zero, and the estimated nonzero part. The latter part is actually the estimated support set where the true support set of global minimizers of (2.9) lies in. Moreover, based on Remark 3, one could identify the support set of a global minimizer of (2.9) easily. On the other hand, noting that the quadratic penalty problem (4.1) is in general nonconvex, it is usually not easy to find a global minimizer. Fortunately, our numerical results demonstrate that in many cases, the projected gradient method can return a solution with accurate support set.

Note that the relaxation problem (3.1) does not take any sparsity into account. However, as shown in Theorem 3.1 at least one of the global minimizers of the relaxation problem (3.1) is a global minimizer of the original problem (2.9). By the quadratic penalty method, we can indeed identify the support set of one global minimizer of (2.9) under reasonable assumptions.

Remark 5

We focused on the problem (2.9) so far. One may wonder whether the theoretical results can be extended to (2.10). It turns out that the extension is not trivial and the analysis becomes more challenging and complicated due to the equality constraints $e^{T}\bar{x}_{i}-1=0$ and $\hat{e}_{i}^{T}x-1=0$. We leave it as a topic to study in future. However, as we will demonstrate in the numerical part, the algorithm designed here can also be applied to solving the relaxation problem of (2.10).

5 Numerical Results

In this section, we will evaluate the performance of our algorithm and compare it with several state-of-the-art approaches for hypergraph matching.

5.1 Implementation issues

Our algorithm is termed as QPPG, which is the abbreviation of Quadratic Penalty Projected Gradient method. Basically, we run Algorithm 2 (referred to as outer iterations) and solve the subproblem (4.7) by calling Algorithm 3 (referred to as inner iterations). In practice, we only execute an inexact version of Algorithm 2 by one step. QPPG2 means that Algorithm 2 is applied to permutation constraints (1.2). For TM duchenne2011tensor , RRWHM lee2011hyper , HGM zass2008probabilistic , and BCAGM nguyen2016efficient , we use the authors’ MATLAB codes and C++ mex files. Our algorithm is implemented in MATLAB (R2015a), while tensor vector multiplications are computed with C++ mex files. All the experiments are preformed on a Dell desktop with Intel dual core i7-4770 CPU at 3.40 GHz and 8GB of memory running Windows 7.

In Algorithm 2, set $\sigma_{0}=10$ and the initial point $x^{0}$ as the vector with all entries equal to one. Update $\sigma_{k}$ as

[TABLE]

where $h_{i}^{k}=e^{T}\bar{x}_{i}^{k}-1$ and $\underline{h}^{k}$ is the maximal value of $\sum_{i}|h_{i}^{k}|$ for five consecutive steps. We stop Algorithm 2 if one of the following conditions is satisfied: (a) $|\Gamma^{k}|$ is less than $1.2n_{1}$; (b) $|\Gamma^{k}|$ stays unchanged for ten consecutive steps. As for the output, each $x^{k}$ returned by different algorithms is projected to its nearest binary assignment matrix by Algorithm 1 except HGM and BCAGM, which output a binary assignment matrix directly). The parameters in Algorithm 3 are $\hbox{Tol}=10^{-5}$, $\rho=10^{-6}$, $\epsilon=10^{-2}$, and $\beta=0.5$. $\eta^{j}$ is chosen as $\eta^{j}=\frac{n_{1}}{\|g_{\bar{I}_{j}}(x^{j})\|_{\infty}}$. The positive definite diagonal matrices $N$ and $Z^{j}$ are set to be the identity matrix.

Generating Tensor $\mathcal{A}$. Note that $\mathcal{A}\in\mathbb{R}^{n\times n\times n}$ contains $n^{3}$ elements. Fortunately, in hypergraph matching, as analyzed in Proposition 2, $\mathcal{A}$ has special structures. Further, $\mathcal{A}$ is also sparse. There are three steps to generate $\mathcal{A}$. The first step is to construct hyperedges $E_{1}$ and $E_{2}$, where each hyperedge connects three different points. The hyperedges in $E_{1}$ are generated by randomly selecting three points in $V_{1}$. We fix $|E_{1}|$ as $n$. $E_{2}$ contains the nearest triples to elements in $E_{1}$, and is generated following the nearest neighbour query approach in duchenne2011tensor ; nguyen2016efficient . The second step is to generate $\mathcal{B}$. Note that the number of nonzero entries in $\mathcal{B}$ are at most $|E_{1}||E_{2}|$, which will be large even for moderate $|E_{1}|$ or $|E_{2}|$. In fact, for each hyperedge in $E_{1}$, we only use $s$ nearest hyperedges in $E_{2}$ to construct $\mathcal{B}$. In other words, $\mathcal{B}$ is calculated by

[TABLE]

where $f_{l_{1}j_{1}k_{1}}$ and $f_{l_{2}j_{2}k_{2}}$ are feature vectors determined by hyperedges $(l_{1},j_{1},k_{1})$ and $(l_{2},j_{2},k_{2})$, and $\gamma=\frac{1}{\text{mean}(\|f_{l_{1}j_{1}k_{1}}-f_{l_{2}j_{2}k_{2}}\|)}$ 111$\hbox{mean}(\|f_{l_{1}j_{1}k_{1}}-f_{l_{2}j_{2}k_{2}}\|)=\frac{\sum_{\mathcal{B}_{l_{1}l_{2}j_{1}j_{2}k_{1}k_{2}}>0}\|f_{l_{1}j_{1}k_{1}}-f_{l_{2}j_{2}k_{2}}\|}{\hbox{number of }B_{l_{1}l_{2}j_{1}j_{2}k_{1}k_{2}}>0}$. is a normalization parameter. Here, for each $(l_{1},j_{1},k_{1})\in E_{1}$, the $s$ nearest neighbours are the $s$ smallest solutions for $\min_{(l_{2},j_{2},k_{2})\in E_{2}}\|f_{l_{1}j_{1}k_{1}}-f_{l_{2}j_{2}k_{2}}\|$. Then $\mathcal{A}$ can be obtained according to (2.2). The number of nonzero elements is $O(sn)$, which is linear in $n$. Therefore, $\mathcal{A}$ is a sparse tensor.

We evaluate the numerical performance mainly from the following three aspects: (1) ‘Accuracy’: denoting the ratio of successful matching, calculated by

[TABLE]

(2) ‘Matching Score’: calculated by $\frac{1}{6}\mathcal{A}(x_{B}^{k})^{3}$, where $x_{B}^{k}$ is the nearest binary assignment vector of $x^{k}$ generated by Algorithm 1; (3) ‘Running Time’: the total CPU time in seconds. For each algorithm (except BCAGM), we only count the computing time for solving the optimization problem. However, BCAGM has to compute all elements in $\mathcal{A}$ to obtain results with high accuracy. Therefore, the running time for BCAGM contains two parts: generating $\mathcal{A}$ with all elements and solving the optimization problem.

Role of Sparsity of $\mathcal{A}$. To see this, we test different values of $s$ on the examples from the CMU house dataset222Downloaded from http://vasc.ri.cmu.edu/idb/html/motion/house/, which has been widely used in literature duchenne2011tensor ; lee2011hyper ; nguyen2016efficient ; zhou2015local . For all examples, there is $n_{1}=30$ and $n_{2}=30$. We take all 111 pictures with labels from 0 to 110, which are the same house taken from slightly different viewpoints. That is, two houses with close labels are similar. For each picture with label $v_{1}$, we match it with $v_{1}+60$. In other words, matching picture $v_{1}$ with $v_{1}+60$ is a test problem. Then we change $v_{1}$ from [math] to $50$ to produce $51$ test examples. To save time of generating input data $\mathcal{A}$, only elements with $l\leq j\leq k$ are computed in $\mathcal{A}$, and the time consumed is denoted by ‘GTensor’. The average results for the test examples are reported in Figure 5.1. One can see that CPU time for generating tensor is not neglectble comparing with CPU time for solving the problem. On the other hand, the accuracy stays almost unchanged for $s\geq 100$. Note that the matching score will be larger when $s$ increases. It is reasonable as a denser $\mathcal{A}$ will result in a larger objective function. Therefore, we set $s=100$ in all the following tests.

Role of Upper Bound $M$. To see the role of $M$, numerical tests are performed on the synthetic data following the approach in duchenne2011tensor ; nguyen2016efficient . Firstly, $n_{1}$ points in $V_{1}$ are sampled following the standard normal distribution $\mathcal{N}(0,1)$. Secondly, points in $V_{2}$ are computed by $V_{2}=TV_{1}+\epsilon$, where $T\in\mathbb{R}^{n_{1}\times n_{1}}$ is a transformation matrix, and $\epsilon\in\mathbb{R}^{n_{1}}$ is the Gaussian noise. We choose $n_{1}=n_{2}$ ranging from 20 to 100, and $M$ from 1 to 10000. All experiments are executed for 100 times, and the average results are reported in Figure 5.2.

We can see that $M=1000$ or $M=10000$ produces competitive results, while $M\leq 100$ is not good for large problems in terms of both accuracy and CPU time. A possible reason is that small $M$ might lead to less flexibility for the entries in $x$. Hence, in the following results, we choose $M=10000$.

5.2 Performance of QPPG and QPPG2

In this subsection, we will illustrate the performance of our algorithm with synthetic data discussed above. We set $n_{1}=n_{2}=30$. Figure 5.3 shows the information while running Algorithm 2, including the accuracy, matching score and size of support set at $x^{k}$.

From Figure 5.3, one can find that $|\Gamma^{k}|$ keeps unchanged in the first few steps, and then drops rapidly from $n_{1}^{2}$ to $n_{1}$, while both accuracy and matching score reach their maximum value within five steps. It shows the potential of our algorithm for identifying the exact support set quickly, even during the process of iteration. This motivates us to stop our algorithm when $|\Gamma^{k}|$ is small enough, or stay unchanged for several iterations.

We also report the magnitude of entries in $x^{k}$ at several selected steps of QPPG in Figure 5.4. The algorithm stops at $k=39$. One can see that $|\Gamma^{k}|$ is decreasing. At the final step, the solution is sparse. This coincides with Remark 4, i.e., the magnitudes of the returned solution by our algorithm clearly fall into two parts: the estimated active part, which is usually close to zero, and the estimated nonzero part, which is the support set we are looking for.

5.3 CMU house dataset

In this subsection, we will test our algorithms on the CMU house dataset. Similar to Section 5.1, we try to match picture $v_{1}$ with $v_{2}$. As $v_{1}$ and $v_{2}$ change, we deal with different hypergraph matching test problems. For a fixed value $v=|v_{1}-v_{2}|$, we set $v_{1}=0,\ldots,110-v$ and $v_{2}=v,\ldots,110$. The total number of test examples is $111-v$. We test these examples, and plot the average results for each $v$ in Figure 5.5. One can see that most algorithms (except HGM) achieve good performance in terms of both accuracy and matching score. In terms of CPU time, QPPG and QPPG2 are competitive with HGM and TM, and faster than other methods.

We also compare QPPG with other algorithms on CMU house dataset with $n_{1}=20$ and $n_{2}=30$. The results are obtained in a similar way as that for Figure 5.5, and are shown in Figure 5.6. One can see that QPPG performs well in both accuracy and matching score. As for CPU time, all the algorithms are competitive since the maximum time is about 0.06s. Figure 5.7 shows the matching results for two houses with $v_{1}=0$ and $v_{2}=60$.

5.4 Large dimensional synthetic data

In this section, large dimensional problems in the fish dataset333Downloaded from http://www.umiacs.umd.edu/$\sim$zhengyf/PointMatching.htm are used to test our algorithms. We use all 100 examples in the subfolder res_fish_def_1. For each example, $V_{1}$ is the set of target fish, and $V_{2}$ is the set of deformation fish. The number of points in each set is around 100. Our task is to match the two sets. We select $n_{1}=n_{2}=10,20,\ldots,100$ points randomly from each fish (for fish with less than 100 points, we use all the points). The average results are shown in Figure 5.8. It can be seen that our algorithm is competitive with other methods in terms of accuracy, matching score and CPU time. One of the matching results is shown in Figure 5.9.

Furthermore, synthetic data explained in Section 5.1 is also used to test these algorithms. All the algorithms are tested except BCAGM, as their codes run into memory troubles for large-scale problems. We choose $n_{1}=n_{2}$ from 50 to 300, and repeat the tests for 100 times. The average results are reported in Figure 5.10. One can see that QPPG and QPPG2 perform comparably well with RRWHM in terms of both accuracy and matching score for $n_{1}$ less than or equal to 200. For $n_{1}$ greater than or equal to $250$, the running time for QPPG and QPPG2 increases slowly as $n_{1}$ increases, which implies that the proposed algorithm can deal with large-scale problems while returning good matching results.

6 Conclusions

In this paper, we reformulated hypergraph matching as a sparse constrained optimization problem. By dropping the sparse constraint, we showed that the relaxation problem has at least one global minimizer, which is also the global minimizer of the original problem. Aiming at seeking for the support set of the global minimizer of the original problem, we allowed violations of the equality constraints by penalizing them in a quadratic form. Then a quadratic penalty method was applied to solve the relaxation problem. Under reasonable assumptions, we showed that the support set of the global minimizer in hypergraph matching can be identified correctly without driving the penalty parameter to infinity. Numerical results demonstrated the high accuracy of the support set returned by our method.

Acknowledgements

The authors would like to thank Dr. Yafeng Liu from Academy of Mathematics and Systems Science, Dr. Bo Jiang from Nanjing Normal University, and Dr. Lili Pan from Beijing Jiaotong University for discussions and insightful comments on this paper. We are also grateful to two anonymous reviewers for their valuable comments, which further improved the quality of this paper.