Continuous optimization methods for the graph isomorphism problem Open Access

Abstract

The graph isomorphism problem looks deceptively simple, but although polynomial-time algorithms exist for certain types of graphs such as planar graphs and graphs with bounded degree or eigenvalue multiplicity, its complexity class is still unknown. Information about potential isomorphisms between two graphs is contained in the eigenvalues and eigenvectors of their adjacency matrices. However, symmetries of graphs often lead to repeated eigenvalues, so that associated eigenvectors are determined only up to basis rotations, which complicates graph isomorphism testing. We consider orthogonal and doubly stochastic relaxations of the graph isomorphism problem, analyze the geometric properties of the resulting solution spaces, and show that their complexity increases significantly if repeated eigenvalues exist. By restricting the search space to suitable subspaces, we derive an efficient Frank–Wolfe based continuous optimization approach for detecting isomorphisms. We illustrate the efficacy of the algorithm with the aid of various highly symmetric graphs.

1. Introduction

Graphs and networks play a central role in many applications ranging from computer science, engineering, chemistry, and biology to the social sciences. Determining whether or not two graphs are isomorphic, i.e., have exactly the same structure, is called the graph isomorphism problem [1–3]. It has been shown that the problem can be solved in polynomial time if the graphs are planar [4] or have bounded degree [5] or eigenvalue multiplicity [6]. Furthermore, Babai recently proved that the graph isomorphism problem is solvable in quasi-polynomial time [7] and summarizes its current status in [3] as ‘not expected to be NP-complete, yet not known to be solvable in polynomial time.’ A historical perspective and a concise overview of different techniques for solving the graph isomorphism problem can be found in [8].

It is well-known that spectral properties of graphs contain crucial information about potential isomorphisms. If the graphs have distinct eigenvalues, then checking whether a permutation exists that transforms the eigenvectors of one graph into the eigenvectors of the other graph is already sufficient. Repeated eigenvalues, however, complicate graph isomorphism testing since the eigenvectors are then only determined up to basis rotations and it is not possible anymore to directly compare their entries. The graph isomorphism problem can be regarded as a combinatorial optimization problem, which requires minimizing a cost function over the set of all permutation matrices. Different relaxations of the graph isomorphism problem have been proposed in the past, where the set of permutation matrices is either replaced by the set of orthogonal [9–11] or doubly stochastic [12–15] matrices. The applicability of such relaxations, however, is not well understood [13]. It has been shown in [13] that the convex relaxation is guaranteed to find the exact isomorphism, provided that the graphs are asymmetric and friendly. This has been generalized to asymmetric graphs with unfriendly eigenvectors that have certain sparsity patterns in [14]. The convex relaxation is called convex exact if it does not admit any additional solutions. For which types of symmetries the relaxation is (and in general is not) exact was analyzed in [15]. We extend these results and consider in particular potentially highly symmetric graphs. We first show how repeated eigenvalues increase the set of feasible solutions of the orthogonal and convex relaxations of the graph isomorphism problem, analyze the geometric properties of these spaces, and then propose an efficient algorithm for detecting graph isomorphisms using a concave reformulation of the relaxed optimization problem that penalizes non-binary matrices. Numerical results for a set of guiding examples and benchmark problems illustrate the efficacy of the proposed algorithm. A method that is similar in spirit is described in [12], where a convex–concave relaxation along with a path following approach is used to solve graph matching and quadratic assignment problems. It was recently shown in [16] that if the signal-to-noise ratio is sufficiently large, then the semidefinite relaxations of these problems are exact and recover the true solutions. We focus on the graph isomorphism problem, which allows us to exploit additional information encoded in the null spaces of associated matrices, and construct a direct approach to solve the resulting relaxed optimization problems.

2. The graph isomorphism problem

Let |$ A $| and |$ B $| be the adjacency matrices of two weighted undirected graphs and ⁠, respectively, where and are the sets of vertices and and the sets of edges. We assume that ⁠, i.e., the graphs share the same vertices. Our goal is to check if the two graphs are isomorphic.

 

Given the permutation |$ \pi $|⁠, we can construct the permutation matrix |$ P $| and vice versa using

If and are isomorphic, we write and define

to be the set of all isomorphisms. An isomorphism from a graph to itself is called an automorphism and the set of automorphisms forms a group under composition (or matrix multiplication), which we will denote by ⁠. Graphs with nontrivial automorphism groups are called symmetric. Examples of highly symmetric graphs are shown in Fig. 1, see also Example 4.18.

We can formulate the graph isomorphism problem as an optimization problem of the form

where |$ \left \lVert \cdot \right \rVert _{F} $| denotes the Frobenius norm, so that if and only if |$ c_{\mathcal{P}}^{} = 0 $|⁠. The solution of the optimization problem, however, is in general not unique if the graphs have symmetries. Since

minimizing the Frobenius norm in (1) is equivalent to maximizing the trace of |$ A^\top \, P^\top \, B \, P $|⁠, which illustrates that the graph isomorphism problem or, more generally, the graph matching problem is closely related to the quadratic assignment problem [17, 18]. In what follows, we will analyze different relaxations of the graph isomorphism problem.

3. Orthogonal relaxation

Let |$ I_{n} $| be the |$ n \times n $| identity matrix and let

denote the set of all orthogonal matrices, which forms a nonlinear disconnected manifold (since |$ \det (Q) = \pm 1 $|⁠). Instead of solving the combinatorial optimization problem (1), we now consider the relaxed problem

which is also known as the two-sided orthogonal Procrustes problem [19, 20].

 

If |$ \lambda _{A}^{(k)} $| is a distinct eigenvalue, then |$ U_{A}^{(k)} $| is the associated eigenvector. If, on the other hand, |$ \lambda _{A}^{(k)} $| is a repeated eigenvalue, then the columns of |$ U_{A}^{(k)} $| form an orthogonal basis of the corresponding eigenspace. Repeated eigenvalues can be caused by symmetries of the graph. Assume that |$ \big (\lambda _{A}^{(k)}, u_{A}^{(k)}\big ) $| is an eigenvalue–eigenvector pair and let ⁠, i.e., |$ P \, A = A \, P $|⁠, then

and |$ P \, u_{A}^{(k)} $| is also an eigenvector associated with |$ \lambda _{A}^{(k)} $|⁠. If these two eigenvectors are linearly independent, this implies that the eigenspace is at least two-dimensional. The larger the automorphism group, the higher the likelihood that the associated permuted eigenvectors are linearly independent. It was shown in [21, 22] that if all eigenvalues of a graph are distinct, then every non-trivial automorphism has order |$ 2 $|⁠. Nevertheless, repeated eigenvalues are not necessarily related to symmetries. The graph in Fig. 2(a), for example, is asymmetric but has repeated eigenvalues. Since it is in practice difficult to determine the size of the automorphism group, a related property called friendliness, which can be verified easily, was proposed in [13].

 

Every friendly graph is automatically asymmetric as shown in [13]. The other direction, however, is not true. The Frucht graph [23] illustrated in Fig. 2(b), for instance, is asymmetric but not friendly, see also [14]. In fact, any asymmetric regular graph will be highly unfriendly since is an eigenvector and the remaining |$ n-1 $| eigenvectors are orthogonal to it. In [11], we defined a stronger property, which can sometimes be used to disambiguate the signs of unfriendly eigenvectors.

 

This definition allows us to assign a canonical sign to unambiguous eigenvectors. An ambiguous eigenvector is automatically unfriendly since ⁠, which implies ⁠. The vector |$ u_{A}^{(k)} = [1, -2, 0, 1]^\top $|⁠, for instance, is unfriendly but not ambiguous.

 

That is, isospectral graphs have the same eigenvalues and also their multiplicities are identical. Isomorphic graphs are clearly isospectral.

 

 

Note that we did not explicitly use the orthogonality of the matrices |$ R^{(k)} $|⁠, which is only required so that |$ Q^{*}(R) $| remains orthogonal. A similar proof can be found in [24]. For a more detailed discussion of the equivariance properties of self-adjoint matrices, see also [25].

 

A simple way to check whether two graphs are potentially isomorphic is to first solve the orthogonal Procrustes problem (2). If |$ c_{\mathcal{O}}^{}> 0 $|⁠, the graphs cannot be isomorphic. However, if |$ c_{\mathcal{O}}^{} = 0 $|⁠, this only implies that the graphs are isospectral (since |$ A $| and |$ B $| are in this case similar) but not necessarily isomorphic (i.e., |$ c_{\mathcal{P}}^{} $| might still be greater than zero). The graph isomorphism problem can thus be reformulated as: does the set of feasible solutions |$ Q^{*}(R) $| contain at least one permutation matrix? Finding such a matrix, assuming the graphs are isomorphic, is a difficult nonlinear optimization problem. One potential way to obtain a valid solution of the graph isomorphism problem is to penalize orthogonal matrices that are not permutation matrices. This approach was used in [10] to solve the graph matching problem.

4. Doubly stochastic relaxation

We will now consider a different relaxation of the graph isomorphism problem. Let

be the set of doubly stochastic matrices, which is the convex hull of the set of all permutation matrices |$ \mathcal{P}(n) $|⁠. This is also known as Birkhoff’s theorem [26]. The doubly stochastic relaxation of the graph isomorphism problem is given by

We again have the property that the graphs cannot be isomorphic if |$ c_{\mathcal{D}}^{}> 0 $|⁠.

4.1 Solutions of the relaxed problem

The convex relaxation of the problem admits additional solutions as the following lemma shows.

 

 

It is important to note here that all convex combinations of isomorphisms are solutions of the relaxed problem, but not all doubly stochastic matrices that minimize (3) are automatically convex combinations of isomorphisms as the following example shows.

 

This is related to the aforementioned friendliness property, but most of the eigenvectors are in this case not ambiguous, so that we can eliminate spurious solutions as we will show below.

4.2 Vectorization of the problem

We now rewrite the optimization problem (3) as

using the cyclic property of the trace and the fact that |$ A $| and |$ B $| are symmetric. Note that here the first term and last term are not independent of |$ X $| as they are for the orthogonal relaxation since in general |$ X^\top X \ne I_{n} \ne X \, X^\top $|⁠.

 

 

 

If |$ A $| is symmetric, we have |$ \operatorname{tr}(A \, X^\top B \, X) = x^\top (A \otimes B) \, x $|⁠. The row and column sum constraints can be reformulated as |$ C \, x = d $|⁠, with

and ⁠, where |$ e_{i} \in{\mathbb{R}}^{n} $| denotes the |$ i $|th unit vector. This defines a |$ (2 \, n - 1) $|-dimensional affine subspace.

 

 

This is a convex but in general not strictly convex optimization problem. In particular, all matrices |$ X $| defined in Lemma 4.1 are optimal solutions.

 

 

4.3 Isospectral graphs

If |$ A $| and |$ B $| are isospectral, Theorem 4.6 immediately implies that the rank of |$ H $| is at most |$ n (n-1) $| since at least |$ n $| of the eigenvalues of |$ H $| are zero. If the matrices are isospectral and also have repeated eigenvalues, the rank will be accordingly lower.

 

 

 

 

 

We can now interpret the relaxed graph isomorphism problem (3) as follows: Find a stochastic matrix |$ X, $| so that its vectorization |$ x $| is contained in the null space of |$ H $|⁠. This also shows that if the graphs have repeated eigenvalues, the search space is much larger since the dimension of the null space increases as shown in Corollary 4.7.

 

 

The structure of the solutions is, as expected, similar to the one obtained for the orthogonal relaxation considered in Section 3. So far, we have not taken the constraint that |$ X $| must be doubly stochastic into account. Let

where |$ w_{A}^{(k)}, w_{B}^{(k)} \in{\mathbb{R}}^{\mu _{A}^{(k)}} $|⁠, then this requirement translates to

Condition (5) is closely related to the definition of friendliness, which requires that all entries of |$ w_{A} $| and |$ w_{B} $| are nonzero and that |$ S $| is a diagonal matrix. Additionally, solutions have to be non-negative, i.e.,

Matrices that satisfy (5) but not (6) are sometimes called pseudo-stochastic. Since |$ S $| is a block-diagonal matrix, we can decompose (5) into |$ m $| separate equations. This results in

If |$ S^{(k)} $| is an orthogonal matrix, then these equations are trivially satisfied (see also Section 3), but there are in general additional solutions. This is in particular the case if |$ w_{A}^{(k)} = w_{B}^{(k)} = 0 $|⁠.

 

 

If the unique optimal solution is a permutation matrix, then ⁠. A similar result was also shown in [13]. Geometrically, this means that the affine subspace given by the constraints |$ C \, x = d $| and the null space of |$ H $| intersect in exactly one point |$ x \ge 0 $|⁠.

If |$ A $| and |$ B $| are isospectral but not friendly, then the intersection of the affine subspace defined by |$ C \, x = d $| and the null space of |$ H $| might again form an affine subspace. Let |$ x_{s} $| be one solution of the augmented inhomogeneous system of linear equations

then all vectors of the form |$ x = x_{s} + x_{h} $|⁠, where |$ x_{h} $| is contained in the null space of the augmented matrix, are pseudo-stochastic solutions. We have

since the dimension of the intersection of the two (affine) subspaces must be lower than the sum of their dimensions.

 

This illustrates again how symmetries and repeated eigenvalues complicate graph isomorphism testing since the solution space increases.

 

By exploiting additional properties of the eigenvectors, we can in this case eliminate spurious solutions and obtain again a unique solution.

 

4.4 Penalizing non-binary matrices

The question now is how we can penalize doubly stochastic matrices that are not permutation matrices. We use Lemma 4.9 to turn the cost function into a constraint and optimize a different function instead.

 

 

A similar penalty term was used in the convex-to-concave projection approach in [15]. Since is an eigenvector of any regular graph ⁠, the matrix is of the form |$ X^{*}(S) $|⁠, where only the entry |$ S_{11} $| is nonzero. This solution does not contain any useful information about potential isomorphisms.

 

4.5 Frank–Wolfe algorithm

In order to solve the optimization problem defined in Lemma 4.15, we use the Frank–Wolfe algorithm [27, 28], which for a constrained convex optimization problem of the form

defined on a convex set |$ \mathcal{C} $| consists of the following two steps:

1. Given a feasible initial condition |$ x^{(k)} \in \mathcal{C} $|⁠, compute an optimal solution |$ y^{*} $| of the linear problem

   $$ \begin{align*}& \min_{y \in \mathcal{C}} \big\langle y, \nabla f(x^{(k)}) \big\rangle. \end{align*} $$
2. Define |$ x^{(k+1)} $| to be a convex combination of |$ x^{(k)} $| and |$ y^{*} $|⁠, i.e.,

   $$ \begin{align*}& x^{(k+1)} = (1 - \gamma) \, x^{(k)} + \gamma \, y^{*}, \end{align*} $$

   where |$ 0 \le \gamma \le 1 $| is chosen so that |$ f\big (x^{(k+1)}\big ) $| is minimized.

The cost function |$ f(x) = -x^\top x $| we seek to minimize is not convex but strictly concave and the global minima are permutation matrices. Since we are just interested in showing that at least one permutation matrix exists that solves the graph isomorphism problem, the non-uniqueness of the optimal global solution is not a problem. The gradient of the cost function is simply |$ \nabla f(x) = -2 \, x $| and, ignoring the constant factor 2, in each iteration of the Frank–Wolfe algorithm, we solve

to obtain a new solution |$ x^{(k+1)} $| as described above, where |$ x^{(0)} $| is the solution of the relaxed problem defined in Lemma 4.5 (if no solutions exist, then ⁠). Since |$ x^{(k)} $| and |$ y^{*} $| are feasible solutions, the convex combination |$ x^{(k+1)} $| also satisfies the constraints. That is, the Frank–Wolfe algorithm is guaranteed to compute a doubly stochastic matrix whose vectorization is contained in the null space of |$ H $|⁠.

 

The results illustrate that the new cost function penalizes doubly stochastic matrices that are not permutation matrices and forces the Frank–Wolfe solver into the vertices of the convex set.

 

Even for highly symmetric graphs the Frank–Wolfe algorithm is able to find one of the global minima of the cost function. In general, the solver will not converge after just one iteration. For the Frucht graph, for instance, see Example 4.2, the algorithm converges after two iterations, although the solution is in this case uniquely defined.

5. Numerical results

In this section, we will present numerical results for two challenging highly symmetric graphs.

5.1 Paley graph

We first consider the Paley graph [29] of size |$ n = 29 $| shown in Fig. 5(a), whose automorphism group contains |$ 406 $| permutations, see [30] for a detailed description of the automorphism groups of Paley graphs. The eigenvalues of the graph are |$ \lambda _{A} = \left [14, \frac{1}{2}\sqrt{29} - \frac{1}{2}, -\frac{1}{2}\sqrt{29} - \frac{1}{2}\right ] $| with multiplicities |$ \mu _{A} = [1, 14, 14] $|⁠. That is, |$ \operatorname{rank}(H) = 448 $|⁠, whereas |$ n^{2} = 841 $|⁠. Despite the many symmetries, the Frank–Wolfe solver quickly converges to a permutation matrix as shown in Fig. 5(b–e). The number of iterations depends on the random permutation we choose to generate ⁠.

5.2 Biggs–Smith graph

The automorphism group of the Biggs–Smith graph [32] of size |$ n = 102 $|⁠, see Fig. 6(a), contains 2448 elements. The adjacency matrix |$ A $| is shown in Fig. 6(b) and its spectrum in Fig. 6(c). The multiplicities of the eigenvalues are |$ \mu _{A} = [1, 9, 18, 16, 17, 16, 9, 16] $|⁠, which implies that |$ \operatorname{rank}(H) = 8860 $|⁠. The Frank–Wolfe solver converges after just one iteration to a permutation matrix as shown in Fig. 6(d–e).

6. Conclusion

The main advantage of the proposed approach is that the graph isomorphism problem is turned into a continuous optimization problem, which in general tend to be easier to solve than combinatorial optimization problems. Although orthogonal and doubly stochastic relaxations have been considered before, we extended these techniques to strongly regular graphs and showed how symmetries and repeated eigenvalues increase the dimensions of the spaces of feasible solutions. We believe that continuous relaxations will not only open new avenues for deriving efficient numerical algorithms for the graph isomorphism problem, but also for gaining insights into the geometric properties of the sets of feasible solutions and the inherent symmetries of graphs.

Although the Frank–Wolfe solver typically converges within a few iterations, the algorithm sometimes gets stuck in local minima, which are often convex combinations of just a few optimal permutation matrices. We also implemented SQP-based solvers, but they suffer from the same problem. Small random perturbations of the initial condition in a feasible direction seem to alleviate this problem. Also using different (i.e., more aggressive) penalty functions, adding barrier functions or continuation methods might help avoid local minima. A detailed analysis of the convergence properties and improved search-direction strategies are beyond the scope of this paper and will be considered in future work.

Further open problems include: How sensitive is the algorithm with respect to numerical errors and is it possible to develop sparse solvers? What is the worst-time complexity of the algorithm and for which types of graphs is the problem guaranteed to be solvable in polynomial time? Can we exploit properties of the orthogonal and convex relaxations at the same time? Another interesting question is whether this approach can be extended to the graph matching problem by including not only the vectors contained in the null space of |$ H $|⁠. The derived optimization problem formulation might also be amenable to quantum annealers, provided that the constrained problem can be rewritten as a quadratic unconstrained binary optimization problem.

Acknowledgements

We would like to thank the reviewers for the helpful comments and suggestions.

Funding

QuantERA II programme funded by the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 101017733 to P.G.).

Data availability

The Matlab code and examples that support the findings presented in this paper are available at https://github.com/sklus/FrankWolfeGI.

Footnotes

References