# Permutation-Invariant Subgraph Discovery

Raghvendra Mall, Shameem A. Puthiya Parambath, Ting Yu, Sanjay Chawla
Qatar Computing Research Institute, Hamad Bin Khalifa University, Doha, Qatar
and
Han Yufei
INRIA

###### Abstract

We introduce Permutation and Structured Perturbation Inference (PSPI), a new
problem formulation which abstracts many graph matching tasks that arise in
systems biology. PSPI can be viewed as a robust formulation of the permutation
inference or graph matching, where the objective is to find a permutation
between two graphs under the assumption that a set of edges may have undergone
a perturbation due to an underlying cause. For example, suppose there are two
gene regulatory networks $\mathbf{X}$ and $\mathbf{Y}$ from a diseased and
normal tissue respectively. Then, the PSPI problem can be used to detect if
there has been a structural change between the two networks which can serve as
a signature of the disease. Besides the new problem formulation, we propose an
ADMM algorithm (STEPD) to solve a relaxed version of the PSPI problem. An
extensive case study on comparative gene regulatory networks (GRNs) is used to
demonstrate that STEPD is able to accurately infer structured perturbations
and thus provides a tool for computational biologists to identify novel
prognostic signatures. A spectral analysis confirms that STEPD can recover
small clique-like perturbations making it a useful tool for detecting
permutation-invariant changes in graphs.

## Introduction

Given two graphs $\mathbf{X}$ and $\mathbf{Y}$ defined on a common vertex set
$\mathcal{V}$, can we _infer_ a permutation matrix $\mathbf{P}$ and a
perturbation matrix $\mathbf{Z}$ s.t.
$\mathbf{Y}\sim\mathbf{P}\mathbf{X}\mathbf{P}^{\prime}+\mathbf{Z}$.

We refer to the above as the PSPI (Permutation and Structured Perturbation
Inference) problem, which can be considered analogous to the famous robust
principal component analysis problem introduced in (Candès et al. 2011), where
the objective is to decompose a rectangular matrix $\mathbf{L}$ into a sum of
a low rank matrix $\mathbf{M}$ and a structured sparse matrix $\mathbf{S}$.
However the PSPI problem is substantially harder due to the combinatorial
nature of the search over permutation matrices. There are $n!$ permutations on
a graph of size $n$.

PSPI is a new problem and abstracts many practical tasks. While our primary
motivation is from cancer research, PSPI can be applied in many network
analysis scenarios. To understand the nature of cancer, scientists often
compare the gene regulatory networks of healthy ($\mathbf{X}$) and diseased
samples ($\mathbf{Y}$) (Mall et al. 2018). Studies have shown that regulatory
networks undergo some amount of localized re-wirings as cancer progresses.
Given two such regulatory networks, it is important to detect not only the
correspondence between the nodes of the networks, but also structured
localized perturbations (Zhang, Horvath et al. 2005; Wallace, Martin, and Ambs
2011; Mall et al. 2017). Even though gene networks are vertex annotated, genes
are known to take over the functionality of other genes making a direct
difference comparison between two networks often misleading. Thus, detecting
perturbations in a graph without taking permutation into account is
fundamentally a mis-identification problem.

Another example can be found in neurology research. $\mathbf{X}$ and
$\mathbf{Y}$ represent two networks over the same set of neurons from healthy
and diseased brain tissue respectively. It was shown in (Demarin et al. 2016)
that neurons in the brain could “reorganize” after a traumatic event and take
over its role allowing normal functioning of an individual. Usually, the
combinations of genes and neurons that provide specific biological functions
have relatively fixed combinatorial patterns, compared to random mutation.
These combination patterns form structured cliques in the graphs. Matching the
graphs of genes and neurons helps identify these structured and localized
changes for downstream research.

Example: Consider the toy example given in Fig 1.

Figure 1: The left plot shows the original random geometric graph $\mathbf{X}$
and middle plot depicts the isomorphic adjacency matrix with random noise,
permutation (first 50 nodes) and structured perturbation $\mathbf{Y}$. The
perturbation corresponds to structured rewiring in $\mathbf{X}$. Our STEPD
approach can correctly infer the perturbation matrix corresponding to
$\hat{\mathbf{Z}}$.

There is an adjacency matrix ($\mathbf{X}$) of a random geometric graph of
$500$ nodes (leftmost) and a permuted and perturbed version of it
($\mathbf{Y}$). The first $50$ nodes in the isomorphic $\mathbf{Y}$ are
permuted w.r.t. $\mathbf{X}$ and it additionally contains some random noise
(in practice this happens due to noise in the acquisition process). The goal
of PSPI is: Given the paired networks, we infer the permutation matrix
$\mathbf{P}$ and the structured perturbation matrix $\mathbf{Z}$. The
rightmost image in Fig 1 shows the inferred $\mathbf{Z}$ which is referred as
$\hat{\mathbf{Z}}$.

The problem of PSPI in paired networks is closely related to graph matching.
In fact, the permutation inference corresponds exactly to the graph matching
problem (Conte et al. 2004). In graph matching, one is interested in finding
the correspondence (isomorphism) between the nodes of two graphs such that the
graphs are ‘structurally the same’. In practice, the observed graphs are
subject to random noise due to many factors including the noise in the
acquisition process and one is interested in approximate matching. In general,
graph matching is a difficult combinatorial problem, and the complexity class
of the corresponding decision problem is not yet known. It is worth noting
that computer vision is another domain where graph matching is a popular topic
(Cho, Lee, and Lee 2010; Dym, Maron, and Lipman 2017; Lê-Huu and Paragios
2017; Wang et al. 2018). They focus on recovering the node correspondences
between object components of different scene scenarios. These works do not
consider the impact of noise corruption (Cho, Lee, and Lee 2010; Dym, Maron,
and Lipman 2017; Lê-Huu and Paragios 2017; Wang et al. 2018), or only assume
the existence of random and sparse edge addition/deletion (Yan et al. 2015;
Jia et al. 2016). In contrast, our study focuses on identifying and
reconstructing _structured perturbations implanted in graphs_. The PSPI
problem which we investigate here is significantly different and to the best
our knowledge this is the first attempt to solve this problem.

Our major contributions are the following:

  1. 1.

We introduce a new problem Permutation and Structure Perturbation Inference
(PSPI), for inferring both permutations and structured perturbations in paired
networks.

  2. 2.

We formulate PSPI as a non-linear integer program whose relaxation turns out
to be a bi-convex program and propose an algorithm termed as STructurEd
Perturbation and permutation Detection (STEPD) based on ADMM principles to
solve the bi-convex program.

  3. 3.

We show STEPD performs better than several state-of-the-art graph matching
solutions on several simulated and benchmark Scale-Free and
Erdos-R$\grave{e}$nyi networks.

  4. 4.

We demonstrate that STEPD method can identify biologically relevant structural
differences between real-world networks of healthy and tumor patients.

## Related Work

Most of the prior work related to the PSPI problem is confined to the
permutation inference problem in graphs i.e. graph matching, which has been
extensively studied within the computer vision and bioinformatics. Due to the
combinatorial nature of the problem, exact methods based on full or partial
enumeration do not scale well, and approximate algorithms are used in
practice. We briefly discuss some of the approximation algorithms closely
related to our work. A good review of this topic can be found in (Conte et al.
2004).

In (Umeyama 1988), one of the early work on graph matching, the authors
proposed an algorithm based on spectral decomposition of the adjacency
matrices. The core idea lies in representing the nodes as the orthogonal
eigenvectors of the adjacency matrices in the eigenspace, the vector space
spanned by the eigenvectors. Then the problem of finding the optimal matching
reduces to the problem of finding the permutation matrix. In (Umeyama 1988),
authors proposed to find the permutation matrix by optimizing the euclidean
distance between the rows of the absolute value orthogonal eigenvectors using
techniques like the Hungarian method. In (Singh, Xu, and Berger 2007), the
authors proposed the IsoRank algorithm, similar to PageRank, to find the
maximum common subgraph between two protein networks by associating a
topological similarity score to the nodes, and solving the eigenvalue problem
using the power method.

One of the most commonly used techniques to solve the graph matching problem
approximately, including this work, is based on the relaxation of a problem
specific discrete objective function. The original combinatorial objective
function is relaxed to a tractable optimization problem which can be
efficiently solved. Often the element-wise $\ell_{1}$ or $\ell_{2}$
(Frobenius) norm of the matching error, defined as the number of the adjacency
disagreements between the two graphs i.e.
$\|\mathbf{P}\mathbf{X}-\mathbf{Y}\mathbf{P}\|_{p},p>0$, is used as the
objective function. In (Zaslavskiy, Bach, and Vert 2009), the authors proposed
two algorithms based on convex-concave relaxations. The convex relaxation
based algorithm called QCV is obtained by relaxing the Frobenius norm of the
matching error to the convex set of double stochastic matrix. In the second
PATH algorithm, a non-linear objective function is formed by linearly
combining the convex relaxation term and a concave relaxation term. The final
objective function is solved using the conditional gradient method. Recently,
in (Fiori et al. 2013) an algorithm was proposed for multi-modal graph
matching (MGM) based on convex relaxation of the Frobenius norm of the
matching error. A relaxed version of the problem was solved in (Vogelstein et
al. 2015) using projected gradient descent and conditional gradient in the
context of comparing brain images. Other graph matching methods include
relaxation labeling, replicator equations and tree search etc. Interested
readers are referred to the survey paper by (Conte et al. 2004) and the
references therein. Our proposed algorithm differs from the above discussed
approaches, as we are interested in inferring both structured perturbation and
permutations in paired networks.

Another relevant topic is graph classification. Previous works (Przulj 2006;
Shervashidze and M.Borgwardt 2006; Kashima, Tsuda, and Inokuchi 2003;
Borgwardt and Kriegel 2005; Perozzi, Al-Rfou, and Skiena 2014; Yanardag and
Vishwanathan. 2015; Grover and Leskovec 2016; Hamilton, Ying, and Leskovec
2017; Lee, Rossi, and Kong 2018; Kipf and Welling 2016) focused on deciding
whether two graphs contain components, such as sub-graphs, of similar
structures with random-walk based graph kernels or various neural network
architectures. The output is a binary classification decision. By contrast,
the PSPI problem focuses on more fine-grained information: it aims at
recovering the clique based perturbation and permutation, rather than
producing a binary classification output.

## Problem Description

Given two graphs of size $n$, $\mathcal{X}$ and $\mathcal{Y}$, we denote their
adjacency matrices as $\mathbf{X}$ and $\mathbf{Y}$ respectively. We assume
that the graphs are undirected with no self-loops and parallel edges.
Mathematically, the graph matching problem can be stated as finding a
permutation matrix $\mathbf{P}$ (bijection between nodes), such that the two
graphs are ‘edge preserving’ isomorphic i.e.
$\mathbf{Y}=\mathbf{P}\mathbf{X}\mathbf{P}^{\prime}$. In practice, the
observed graphs are subject to random noise and might not be perfectly
isomorphic. Formally, we can write the problem as
$\mathbf{Y}=\mathbf{P}\mathbf{X}\mathbf{P}^{\prime}+\mathbf{E}$ where
$\mathbf{E}$ represents the random noise. In a more realistic settings, in
graph matching one is interested in finding a permutation matrix which is
‘edge preserving’ in an ‘optimal sense’. The optimality is usually defined in
terms of the matching error which is defined as the number of adjacency
disagreements between $\mathbf{X}$ and $\mathbf{Y}$ quantified using the
Frobenius matrix norm. We use $\|\cdot\|$ to indicate matrix Frobenius norm.
Thus the objective function becomes
${\displaystyle\operatorname*{argmin}_{\mathbf{P}\in\mathcal{P}}\|\mathbf{Y}-\mathbf{P}\mathbf{X}\mathbf{P}^{\prime}\|^{2}}$
where $\mathcal{P}$ is the set of all permutation matrices. In general, the
combinatorial nature of the permutation matrix search makes the problem NP. A
common strategy is to relax the $\mathcal{P}$ to a compact set. The convex
hull of the permutation matrix corresponds to the set of double stochastic
matrix, also called Birkhoff polytope, and the original problem reduces to the
relaxed problem
${\displaystyle\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D}}\|\mathbf{Y}-\mathbf{P}\mathbf{X}\mathbf{P}^{\prime}\|^{2}}$
where $\mathcal{D}$ is the Birkhoff polytope i.e.
$\mathcal{D}=\\{\mathbf{A}\in\Re^{n\times n},\mathbf{A}\geq
0,\mathbf{A}\bm{1}=\bm{1},\mathbf{A}^{\prime}\bm{1}=\bm{1}\\}$, $\bm{1}$ being
the $n$-dimensional vector of 1s.

As stated in the introduction, in many real world networks, in addition to the
random noise represented by spurious edges, we have localized structured
perturbations which indicates ‘re-wirings‘ or new activity in specific regions
of the network. Given two graphs in the form of adjacency matrices
$\mathbf{X}$ and $\mathbf{Y}$, we aim to find the optimal matching between the
graphs and infer only the localized structured perturbations. We propose a
simple extension to the original graph matching problem, as
${\displaystyle\mathbf{Y}=\mathbf{P}\mathbf{X}\mathbf{P}^{\prime}+\mathbf{Z}+\mathbf{E}}$.
Here, $\mathbf{Z}$ captures the localized structured perturbations and can be
conceived as the adjacency matrix of a difference graph and hence symmetric.
We limit ourselves to the case where the structured perturbation consists of
only addition of new clusters of edges but not the removal. Hence,
$\mathbf{Z}$ is a symmetric binary matrix. The binary restriction on
$\mathbf{Z}$ makes the problem combinatorial in $\mathbf{Z}$, hence we relax
$\mathbf{Z}$ to the set of matrices such that entries are between $[0,1]$.

### Optimization Framework

Following (Zaslavskiy, Bach, and Vert 2009; Fiori et al. 2013), we propose an
optimization framework as problem solution. In case of $\mathbf{P}$, similar
to (Fiori et al. 2013), we use an objective function which disregards the
spurious mismatching edges and looks for group structure between
$\mathbf{P}\mathbf{X}$ and $\mathbf{Y}\mathbf{P}$ i.e. the objective function
encourages that $\mathbf{P}\mathbf{X}$ and $\mathbf{Y}\mathbf{P}$ to share the
same active set except for spurious edges. We enforce a group regularization
on $\mathbf{Z}$ to capture any localized group structure present in
$\mathbf{Z}$. We also assume that $\mathbf{Z}$ is a very sparse matrix, as the
perturbation occur at very specific regions of the network. Thus, our
objective function for $\mathbf{Z}$ takes into account the two aspects _(i)_
perturbations are localized but structured and _(ii)_ perturbation appears
only at a small number of regions of the network.

Combining the different objectives, our final objective function becomes,

$\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D},\mathbf{Z}}\sum_{i,j}\big{\|}({\bm{\alpha}}_{ij},{\bm{\beta}}_{ij})\big{\|}_{2}+\|\mathbf{P}\mathbf{X}+\mathbf{Z}\mathbf{P}-\mathbf{Y}\mathbf{P}\|^{2}+\nu\|\mathbf{Z}\|_{1}\\\
~{}~{}~{}~{}~{}~{}~{}~{}+\mu\sum_{i}\|\mathbf{Z}_{i}\|_{2}\text{~{}where~{}}{\bm{\alpha}}=\mathbf{P}\mathbf{X},~{}~{}{\bm{\beta}}=\mathbf{Y}\mathbf{P}$
(1)

where $\nu$ and $\mu$ are regularization co-efficients.

The group lasso on the pair $({\bm{\alpha}},{\bm{\beta}})$ encourages non-
random group of edges to be active, in particular the minimizer of the
objective function (1) with no perturbation is exactly
$\mathbf{P}\mathbf{X}=\mathbf{Y}\mathbf{P}$. The group lasso on $\mathbf{Z}$
promotes only localized structured perturbation to be active and the
$\ell_{1}$ regularization on $\mathbf{Z}$ encourages sparsity.

The optimization problem in Equation (1) is a bi-convex problem. The problem
is convex in $\mathbf{P}$, keeping $\mathbf{Z}$ constant and convex in
$\mathbf{Z}$ keeping $\mathbf{P}$ constant. Since, in general bi-convex
problems do not have closed-form solutions, iterative methods like alternating
direction methods are used in practice. In the next section, we propose an
algorithm to solve the optimization problem given in (1).

### ADMM formulation

Alternating Direction Method of Multipliers is an iterative method to solve
non-linear optimization problems in a Gauss-Seidel fashion (Boyd et al. 2011),
and naturally fits in our optimization framework. Moreover, ADMM based
algorithm is more appealing as the intermediate updates in ADMM scheme
naturally lend to a parallel implementation. Here, we propose an ADMM based
algorithm to solve (1).

Introducing auxiliary variables for $\mathbf{C}=\mathbf{Z}$ the augmented
Lagrangian for (1) becomes:

$\displaystyle\begin{split}\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D},\mathbf{Z}}&\sum_{i,j}\big{\|}({\bm{\alpha}}_{ij},{\bm{\beta}}_{ij})\big{\|}_{2}+\frac{1}{2}\|\mathbf{P}\mathbf{X}-\mathbf{Y}\mathbf{P}+\mathbf{Z}\mathbf{P}\|_{2}^{2}+\nu\|\mathbf{C}\|_{1}\\\
&~{}~{}~{}~{}+\mu\sum_{i}\|\mathbf{C}_{i}\|_{2}+\frac{\rho}{2}\|\mathbf{C}-\mathbf{Z}+\mathbf{D}\|_{2}^{2}+\\\
&~{}~{}~{}~{}~{}~{}~{}\frac{\rho}{2}\|{\bm{\alpha}}-\mathbf{P}\mathbf{X}+\mathbf{U}\|^{2}+\frac{\rho}{2}\|{\bm{\beta}}-\mathbf{Y}\mathbf{P}+\mathbf{V}\|^{2}\end{split}$

here $\mathbf{U},\mathbf{V}$ and $\mathbf{D}$ are related to the Lagrange
multiplier, and $\rho$ is the penalization parameter.

In ADMM, we iteratively update each variable, keeping others constant,
starting with primal variables
${\bm{\alpha}},{\bm{\beta}},\mathbf{C},\mathbf{P},\mathbf{Z}$ and followed by
the dual variables $\mathbf{U},\mathbf{V}$ and $\mathbf{D}$. The complete ADMM
based algorithm is given in Algorithm 1.

Input : Adjacency matrix $\mathbf{X}$ for $\mathbf{G}_{1}$ and $\mathbf{Y}$
for $\mathbf{G}_{2}$, penalty parameter $\rho$

1 Initialize
${\bm{\alpha}}=0,{\bm{\beta}}=0,{\bm{\gamma}}=0,\mathbf{Z}=0,\mathbf{P}=\frac{1}{n}\bm{1}^{\prime}\bm{1},\mathbf{U}=0,\mathbf{V}=0,\mathbf{W}=0$
;

2 repeat

3
$({\bm{\alpha}},{\bm{\beta}})=\operatorname*{argmin}_{{\bm{\alpha}},{\bm{\beta}}}~{}\sum_{i,j}\big{\|}({\bm{\alpha}}_{ij},{\bm{\beta}}_{ij}\big{\|}_{2}+$
$\frac{\rho}{2}\big{\|}{\bm{\alpha}}-\mathbf{P}\mathbf{X}+\mathbf{U}\big{\|}_{2}^{2}+\frac{\rho}{2}\big{\|}{\bm{\beta}}-\mathbf{Y}\mathbf{P}+\mathbf{V}\big{\|}_{2}^{2}$;

4
$\mathbf{C}=\operatorname*{argmin}_{\mathbf{C}}~{}\nu\big{\|}\mathbf{C}\big{\|}_{1}+\mu\sum_{i}\big{\|}\mathbf{C}_{i}\big{\|}_{2}+$
$\frac{\rho}{2}\big{\|}\mathbf{C}-\mathbf{Z}+\mathbf{D}\big{\|}_{2}^{2}$;

5
$\mathbf{P}=\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D}}\frac{1}{2}\big{\|}\mathbf{P}\mathbf{X}+\mathbf{Z}\mathbf{P}-\mathbf{Y}\mathbf{P}\big{\|}_{2}^{2}+$
$\frac{\rho}{2}\big{\|}{\bm{\alpha}}-\mathbf{P}\mathbf{X}+\mathbf{U}\big{\|}_{2}^{2}+\frac{\rho}{2}\big{\|}{\bm{\beta}}-\mathbf{Y}\mathbf{P}+\mathbf{V}\big{\|}_{2}^{2}$;

6
$\mathbf{Z}=\operatorname*{argmin}_{\mathbf{Z}}\frac{1}{2}\big{\|}\mathbf{P}\mathbf{X}+\mathbf{Z}\mathbf{P}-\mathbf{Y}\mathbf{P}\big{\|}_{2}^{2}+$
$\frac{\rho}{2}\big{\|}\mathbf{C}-\mathbf{Z}+\mathbf{D}\big{\|}_{2}^{2}$;

7 $\mathbf{U}=\mathbf{U}+{\bm{\alpha}}-\mathbf{P}\mathbf{X}$;

8 $\mathbf{V}=\mathbf{V}+{\bm{\beta}}-\mathbf{Y}\mathbf{P}$;

9 $\mathbf{D}=\mathbf{D}+\mathbf{C}-\mathbf{Z}$;

10

11until _stopping criterion is not satisfied_ ;

Output : $\mathbf{P},\mathbf{Z}$

Algorithm 1 ADMM Based Algorithm

We now take a closer look at each of these update steps, and propose efficient
methods to solve each of the sub-optimization problems, where the closed-form
solution is not available. We refer to the numbering in Algorithm 1 when the
sub-optimization problems are referred using numbers. The subproblem at line 1
is a well studied group lasso problem, and in (Yuan and Lin 2006), authors
propose a closed-form solution based on the soft-thresholding operator as
given by ${\displaystyle
S(a,\rho)=\Big{[}1-\frac{\rho}{\|a\|_{2}}\Big{]}_{+}a}$. In our settings, $a$
corresponds to the 2-dimensional vector of individual elements of
$\mathbf{P}\mathbf{X}$ and $\mathbf{Y}\mathbf{P}$. Similarly, the solution for
the subproblem at line 1 can be expressed in a closed-form solution. Similar
to the problem at line 1, this problem is a group lasso with an additional
lasso term on the variable matrix. As shown in (Ming et al. 2014), the problem
can be independently solved for each column (or for each row due to
symmetrization) of $\mathbf{C}$, and the closed-form solution can be written
in terms of the soft-thresholding operator as: ${\displaystyle
S(a,\mu,\nu)=\Big{[}\frac{\|a\|_{2}-\mu}{\nu\|a\|_{2}}\Big{]}_{+}a}$. Here $a$
represents the column vector of $\mathbf{C}$.

#### Solving for $\mathbf{P}$

The subproblem at line 1 is a constrained convex optimization problem, where
the domain is constrained to be in the set of double stochastic matrix
(Birkhoff polytope). Unfortunately, no closed-form solution exists for this
problem. Due to the high computation cost of interior-point methods, in
practice first-order methods are preferred to solve such problems. The
projected gradient descent algorithm (Luenberger and Ye 2015) is a popular
technique due to its guaranteed linear convergence rate. But the projection of
$\mathbf{P}$ to the set of double stochastic matrix does not have a closed-
form solution either. Hence employing gradient based methods will not be
computationally efficient. To solve for $\mathbf{P}$, we use linearized
version of the alternating direction method (ADM), recently proposed in (Lin,
Liu, and Su 2011). Linearized ADM approach is extremely useful for solving
ADMM subproblems where a closed-form solution does not exist. Interestingly,
linearization makes the auxiliary variables unnecessary, and hence there is no
need to update them. Moreover, Linearized ADM converges faster than the
traditional ADMM procedure (Lin, Liu, and Su 2011).

Linear ADM is a first order method, where we augment the first order Taylor
approximation of the objective function with the proximal operator. Its
primary advantage is that if the objective function is quadratic, the
augmented objective takes a simple projection on the domain constraint
($\mathcal{D}$ in our case), thus avoiding the iterative gradient updates. We
now give the update formula for $\mathbf{P}$, by analyzing the three terms
appearing in the optimization problem at line 1 independently. By introducing
auxiliary variables ${\bm{\sigma}}$ and ${\bm{\omega}}$, for
$\mathbf{P}\mathbf{X}$ and $\mathbf{Y}\mathbf{P}$, the proximal augmented
first order Taylor approximation at iteration $j+1$ of ADMM procedure becomes,

$\mathbf{P}_{j+1}=\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D}}\frac{1}{2}\|{\bm{\sigma}}-{\bm{\omega}}+\mathbf{Z}\mathbf{P}_{j}\|^{2}+\langle({\bm{\sigma}}-{\bm{\omega}}+\mathbf{Z}\mathbf{P}_{j}),\\\
\mathbf{P}-\mathbf{P}_{j}\rangle+\frac{\eta}{2}\|\mathbf{P}-\mathbf{P}_{j}\|^{2}$

where $\langle\,,\rangle$ is the Frobenius inner product defined as
$\langle\mathbf{A},\mathbf{B}\rangle=\mathbf{\operatorname{Tr}}(\mathbf{B}^{\prime}\mathbf{A})$,
and $\eta$ is the adaptive penalty parameter. Now, the above equation can be
equivalently re-written as (disregarding the constant terms)

$\displaystyle\mathbf{P}_{j+1}=\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D}}\frac{1}{2}\|\mathbf{P}-\big{(}\mathbf{P}_{j}+\tau\mathbf{Z}^{\prime}(\mathbf{M}-\mathbf{Z}\mathbf{P})\big{)}\|^{2}$
(2)

where $\mathbf{M}={\bm{\omega}}-{\bm{\sigma}}$. Note that we do not introduce
any Lagrangian multipliers corresponding to the auxiliary variables
${\bm{\omega}}$ and ${\bm{\sigma}}$. By initializing ${\bm{\omega}}$ to
${\bm{\alpha}}$ and ${\bm{\sigma}}$ to ${\bm{\beta}}$, the Lagrangian
variables $\mathbf{U}$ and $\mathbf{V}$ stand as the Lagrangian variables for
${\bm{\omega}}$ and ${\bm{\sigma}}$ as well. In a similar fashion, the second
and third terms in the line 1 can be equivalently re-written as:

$\displaystyle\mathbf{P}_{j+1}=\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D}}\frac{1}{2}\|\mathbf{P}-\big{(}\mathbf{P}_{j}+\tau(\mathbf{N}-\mathbf{P}\mathbf{X})\mathbf{X}^{\prime}\big{)}\|^{2}$
(3)

where $\mathbf{N}={\bm{\alpha}}+\mathbf{U}$, and

$\displaystyle\mathbf{P}_{j+1}=\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D}}\frac{1}{2}\|\mathbf{P}-\big{(}\mathbf{P}_{j}+\tau\mathbf{Y}^{\prime}(\mathbf{O}-\mathbf{Y}\mathbf{P})\big{)}\|^{2}$
(4)

where $\mathbf{O}={\bm{\beta}}+\mathbf{V}$

Now, combining Equations (2),(3),(4), our final optimization problem for
$\mathbf{P}$ becomes,

$\displaystyle\mathbf{P}_{j+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{P}\in\mathcal{D}}\frac{1}{2}\Big{\|}\mathbf{P}-\frac{\big{(}\mathbf{A}+\mathbf{B}+\mathbf{C})}{3}\Big{\|}^{2}$
(5)

where
$\mathbf{A}=\mathbf{P}_{j}+\tau\mathbf{Z}^{\prime}(\mathbf{M}-\mathbf{Z}\mathbf{P}),\mathbf{B}=\mathbf{P}_{j}+\tau(\mathbf{N}-\mathbf{P}\mathbf{X})\mathbf{X}^{\prime}$
and
$\mathbf{C}=\mathbf{P}_{j}+\tau\mathbf{Y}^{\prime}(\mathbf{O}-\mathbf{Y}\mathbf{P})$.

Equation (5) corresponds exactly to the euclidean projection of the matrix
$\nicefrac{{(\mathbf{A}+\mathbf{B}+\mathbf{C})}}{{3}}$ onto the set of double
stochastic matrix. The problem of projecting a matrix onto the Birkhoff
polytope is a well studied problem, and there exists very simple and efficient
algorithm as shown in (Sinkhorn and Knopp 1967). The algorithm proceeds by
alternately normalizing rows and columns of the given non-negative matrix. In
fact, it is established that Sinkhorn projection algorithm returns a double
stochastic matrix which is optimal according to KL-divergence i.e. the double
stochastic matrix has the lowest KL-divergence distance rather than euclidean
distance. A recent paper (Wang, Li, and Konig 2010) proposed a projection
algorithm which returns a projection which is optimal with respect to the
euclidean projection. We use the algorithm proposed in (Wang, Li, and Konig
2010) to project $\mathbf{P}$ onto the Birkhoff polytope.

#### Solving for $\mathbf{Z}$

The subproblem at line 1 corresponding to $\mathbf{Z}$ is also a constrained
optimization problem. Unlike for $\mathbf{P}$, here we can find a simple
closed-form solution for the projection. The problem is to find the closest
projection of a given matrix $\mathbf{Z}$ to the set of the matrices
$\mathcal{F}$. The set of matrices $\mathcal{F}$ is defined as the set of all
matrices whose entries lie in the interval [0,1]. The problem can be written
as
${\displaystyle\mathbf{Z}^{*}=\arg\min_{\mathbf{A}\in\mathcal{F}}||\mathbf{A}-\mathbf{Z}||^{2}}$.
This is a convex problem, where the objective becomes
${\displaystyle\|\mathbf{A}-\mathbf{Z}\|^{2}=\sum_{i,j}(\mathbf{A}_{ij}-\mathbf{Z}_{ij})^{2}}$
and thus, it is enough to find solutions for $\mathbf{A}_{ij}$ individually.

$\mathbf{A}_{ij}=\begin{cases}\mathbf{Z}_{ij},&\mbox{if}~{}0<\mathbf{Z}_{ij}<1\\\
0,&\mbox{if}~{}-\infty<\mathbf{Z}_{ij}<0\\\
1,&\mbox{if}~{}~{}1<\mathbf{Z}_{ij}<\infty\end{cases}$ (6)

We solve for $\mathbf{Z}$ using the projected gradient descent method. The
gradient of $\mathbf{Z}$ at the point $\mathbf{Z}_{j}$ can be estimated as,

${\bm{\nabla}}\mathbf{Z}\big{|}_{\mathbf{Z}_{j}}=0.5*(\mathbf{Z}_{j}*\mathbf{P}-({\bm{\beta}}-{\bm{\alpha}}))\mathbf{P}^{\prime}+\nu(\mathbf{Z}_{j}-\mathbf{D}-\mathbf{C})$

At the end of each gradient descent step, we project the result to the set of
all matrices whose entries lie in the interval [0,1], using the projection
rule defined in Equation 6. We project the resulting ${\mathbf{Z}}^{*}$ to the
set of symmetric matrix having minimum euclidean distance by using
$\Pi(\mathbf{Z}^{*})=\frac{{\mathbf{Z}}^{*}+{\mathbf{Z}}^{*}{{}^{\prime}}}{2}$.

## Experiments & Analysis

We evaluate the performance of our algorithm on several synthetic and real
world datasets. The objective of our experiments is three-fold: 1) including
the structured perturbation helps reducing the adjacency mismatches thus
obtaining lower matching error; 2) one can infer the structured perturbations
with high precision and recall in simulated networks where the true
$\mathbf{Z}$ is known; and 3) evaluate the structured perturbations as
significant clusters i.e. block-diagonals in a node-cluster membership based
re-ordered $\hat{\mathbf{Z}}$ using the unsupervised F-measure defined as the
harmonic mean of entropy and balance (Mall, Langone, and Suykens 2013).

The probability for cluster $i$ is defined as:
$p_{i}=\nicefrac{{|C_{i}|}}{{n}}$, where $C_{i}$ represents the nodes in
$\hat{\mathbf{Z}}$ belonging to structured perturbation cluster $i$. Entropy
is defined as: $E=\sum_{i=1}^{k}p_{i}\log(p_{i})$, which ranges from $0$, when
all the $n$ nodes are part of an inferred structured perturbation, to a
maximum of $\log(n)$, when each individual node is a structured perturbation
cluster. Normalized entropy is defined as: $NE=\nicefrac{{E}}{{\log(n)}}$ to
make its range $[0,1]$. Similarly, balance is defined as:
$B=\nicefrac{{(\sum_{i}^{k}\frac{|C_{i}|}{\max(|C_{1}|,|C_{2}|,\ldots,|C_{k}|)})}}{{k}}$.
Balance values are high when a few large clusters are identified as structured
perturbation in the inferred $\hat{\mathbf{Z}}$. Thus, there exist a trade-off
between $E$ and $B$. Our goal is to identify the $k$ that yields the optimal
unsupervised F-measure. This gives a good way to measure the performance on
real paired networks where true $\mathbf{Z}$ is unknown. Further details about
using F-measure can be found in (Mall, Langone, and Suykens 2013).

In the first set of experiments, following (Fiori et al. 2013), we compare the
matching error against state-of-the-art graph matching algorithms on different
synthetic graphs. We used the algorithm UNMY (Umeyama 1988), IsoRank (RANK)
(Singh, Xu, and Berger 2007), PATH (Zaslavskiy, Bach, and Vert 2009) and
multi-modal graph matching (MMG) (Fiori et al. 2013) as our baselines. Our
proposed algorithm is denoted as STEPD. In STEPD, the regularization co-
efficients ($\nu$ and $\mu$) are selected from the range
$\\{2^{-5},\cdots,1\\}$ as multiple of 2 that resulted in the best matching
error. It has been shown that ADMM works well without any tuning of the
penalty parameter $\rho$, and a value of 1 works well in most settings (Boyd
et al. 2011). The number of ADMM iteration is set to 150. The matching error
for STEPD is defined as
$\|\mathbf{P}\mathbf{X}+\mathbf{Z}\mathbf{P}-\mathbf{Y}\mathbf{P}\|$, as
$\mathbf{Z}$ is not noise whereas in baselines it is defined as
$\|\mathbf{P}\mathbf{X}-\mathbf{Y}\mathbf{P}\|$.

###  Data

We experimented with $5$ randomized Scale-Free (SF) graphs (with exponent
$\alpha$=$1.5$) and Erdos-R$\grave{e}$nyi graphs (ER) (with parameter
$r$=$0.15$) respectively and the reported results are the average over the $5$
runs. For each experiment, the number of vertices ($n$) was set to $500$, two
structured perturbations (by adding edges to form cliques) were introduced in
$\mathbf{Y}$, one of $50$ nodes while the other of $100$ nodes. We also added
an additional $30$ or $50$ random edges in $\mathbf{Y}$ as noise. We had $2$
set of experiments: 1) the permutation matrix ($\mathbf{P}$) was set to an
identity matrix ($\mathbf{I}$) referred as ‘NP’ and 2) the first $50$ nodes of
$\mathbf{X}$ were permuted referred as ‘P’. Table 1 shows a comprehensive
comparison of STEPD with state-of-the-art graph matching algorithms w.r.t.
matching error. Supplementary Figure 1 depicts the effect of STEPD model
parameters $\nu$ and $\mu$ on precision, recall and F-score, when comparing
the structured perturbations in the original $\mathbf{Z}$ with the inferred
$\hat{\mathbf{Z}}$ for both SF and ER graphs.

Table 1: Synthetic Data Mean Matching Error  Graph | #noise | UMY | RANK | PATH | MMG | STEPD
---|---|---|---|---|---|---
SF (NP) | 30 | 84.07 | 85.01 | 84.24 | 82.60 | 65.19
SF (NP) | 50 | 84.85 | 85.78 | 82.33 | 82.64 | 66.22
SF (P) | 30 | 84.34 | 85.51 | 84.08 | 81.72 | 60.29
SF (P) | 50 | 85.05 | 85.58 | 77.71 | 82.42 | 60.34
ER (NP) | 30 | 172.18 | 173.33 | 103.33 | 134.55 | 118.27
ER (NP) | 50 | 172.49 | 173.24 | 83.31 | 125.14 | 103.54
ER (P) | 30 | 172.19 | 173.29 | 75.13 | 121.99 | 114.39
ER (P) | 50 | 170.96 | 172.54 | 83.80 | 122.25 | 114.81

### Spectral Recovery Evaluation

We demonstrated the ability of STEPD to recover permuted planted cliques in
the presence of bernoulli noise with Synthetic data. We generated Erodos-Reyni
(ER) graph $\mathbf{X}$ of size $100$ and probability $0.5$. We then generated
permuted graphs $\mathbf{P}$ where five nodes were randomly permuted and the
others were kept fixed. We then planted cliques of sizes varying as multiples
of three up to a maximum clique of size $90$. Furthermore, we added noisy
random edges again, using an ER graph of size 100 and probability 0.1. Thus,
the permuted and perturbed matrix was formed as
$\mathbf{Y}=\mathbf{P}\mathbf{X}\mathbf{P}^{\prime}+\mathbf{Z}+\mathbf{E}$.
STEPD was given input $\mathbf{X}$ and $\mathbf{Y}$ and output two matrices
$\hat{\mathbf{P}}$ and $\hat{\mathbf{Z}}$, the inferred permutation and the
structured perturbation matrices. We then generated the matrix
$\hat{\mathbf{Y}}=\hat{\mathbf{P}}\mathbf{X}\hat{\mathbf{P}^{\prime}}+\hat{\mathbf{Z}}$
and computed the relative change in the maximum eigenvalues of
$\hat{\mathbf{Y}}$ and $\mathbf{Y}$ defined as
$\frac{\lambda_{max}(\hat{\mathbf{Y}})-\lambda_{max}(\mathbf{Y})}{\lambda_{max}(\mathbf{Y})}$.
The result of the relative change as a function of the planted clique size are
shown in Figure 2. It is clear that as the size of the clique increases, the
ability of STEPD to recover matrices which are spectrally equivalent to the
original $\mathbf{Y}$ increases. Moreover, even for small planted cliques the
relative error is small.

Figure 2: Relative difference between the maximum eigenvalues between $Y$ and
$\hat{Y}$

### Real Data

We compare the matching errors of different algorithms on two real world data.

#### C. Elegans network

An experiment was performed on the neuronal connectivity networks of two
roundworms belong to the C. Elegans family (White et al. 1986). The nodes of
the network are a common set of $202$ neurons with different number of edges,
$2870$ and $3090$ respectively. The adjacency matrices $\mathbf{X}$ and
$\mathbf{Y}$ are shown in Figure 3.

As shown in (Fiori et al. 2013), we permuted the first $20$ nodes of
$\mathbf{X}$ to obtain $\mathbf{X}_{p}$ s.t.
$\mathbf{Y}=\mathbf{P}\mathbf{X}_{p}\mathbf{P}^{\prime}+\mathbf{Z}$ and our
goal is to infer both $\mathbf{P}$ and $\hat{\mathbf{Z}}\sim\mathbf{Z}$.

Figure 3: An example of two roundworm (C. Elegans) networks demonstrating the
PSPI problem and solution. Each subplot shows the adjacency matrix of the
respective graph. Here the ‘yellow’ color represents the non-zero edges in the
adjacency matrix. STEPD identifies $4$ structured perturbations as block
diagonal components in the inferred $\hat{\mathbf{Z}}$ with a precision of
$0.703$ and recall of $0.626$ w.r.t. the structured perturbations in the
original $\mathbf{Z}$.

Figure 3 shows a node-cluster membership based re-ordered version of
$\mathbf{Z}$ to reflect the structured perturbations in $\mathbf{Z}$ as
approximate block diagonals. The corresponding structured perturbations in the
inferred $\hat{\mathbf{Z}}$ obtained by the proposed STEPD methodology is
shown in Figure 3F. Table 2 compares the matching error of STEPD algorithm
with several graph matching approaches. The STEPD algorithm has the least
matching error.

Moreover, we perform our unsupervised F-measure based analysis to identify
structured perturbations (i.e. clusters) in both $\mathbf{Z}$ and inferred
$\hat{\mathbf{Z}}$. The optimal $k$ for which F-measure was maximum for
$\mathbf{Z}$ was $k=4$ with F-measure$=$0.33. Similarly, the optimal $k$ for
which F-measure was maximum for $\hat{\mathbf{Z}}$ was $k=4$ with
F-measure$=$0.35. This suggests that our STEPD algorithm correctly identifies
the structured perturbation clusters in $\mathbf{Z}$. To have a comprehensive
analysis, we further compared the location of non-zero weight edges i.e. local
structure of the perturbation clusters in $\mathbf{Z}$ and inferred
$\hat{\mathbf{Z}}$ to attain a precision of $0.703$ and recall of $0.626$.

Table 2: Real Data Matching Error  Graph | UMY | RANK | PATH | MMG | STEPD
---|---|---|---|---|---
CEL | 71.57 | 71.36 | 65.18 | 64.25 | 57.98
TCGA | 124.97 | 122.70 | 141.89 | 113.50 | 109.82

#### Glioma Cancer Data

We analyze the structures of two networks extracted from a glioma cancer
dataset obtained from the Cancer Genome Atlas (TCGA)(Johnson, Li, and
Rabinovic 2007). The two networks correspond to IDH-mutant and IDH-wildtype
subtypes of glioma and are labeled $\mathbf{X}$ and $\mathbf{Y}$ respectively.
The networks includes $12,985$ genes comprising $457$ transcription factors
(TF) and $12,895$ target (T) genes, where $\text{TF}\in\text{T}$. The IDH-
mutant network has $13,683$ (TF $\rightarrow$ T) edges while IDH-wildtype
network has $14,158$ (TF $\rightarrow$ T) edges. Our goal is to identify sub-
networks of TFs having a different regulatory program in these glioma
subtypes.

(a) Original sub-network of $286$ TFs which are part of structured
perturbation cluster $C_{1}$.

(b) Original sub-network of $10$ TFs present in structured perturbation
cluster $C_{2}$.

(c) Original sub-network of $10$ TFs present in structured perturbation
cluster $C_{3}$.

Figure 4: Original sub-networks in the IDH-mutant and IDH-wildtype graphs
consisting of the edges belonging to the 3 perturbation clusters in
$\hat{\mathbf{Z}}$ inferred by STEPD algorithm. Here ‘green’ edges represent
edges present in both IDH-mutant and IDH-wildtype. The ‘red’ edges are from
the IDH-mutant graph while ‘blue’ edges come from IDH-wildtype network. The
‘white’ circles represent the TFs and the peripheral nodes are the target (T)
genes. In accordance with our goal of identification of difference in
regulatory networks of IDH-mutant versus IDH-wildtype, we observe that both
perturbation clusters $C_{1}$ and $C_{2}$ have distinct communities
(predominatly ‘red’ or ‘blue’ connections) associated with IDH-mutant and IDH-
wildtype sub-networks respectively. Thus, STEPD algorithm can detect
structured changes between IDH-mutant and IDH-wildtype sub-networks which have
an equivalent biological relevance. We present the detail of the unveiled
relevance in Table.2 of the Supplementary.

Though there exists one-to-one mapping between the nodes in $\mathbf{X}$ and
$\mathbf{Y}$, details about the inferred permutation matrix and genes taking
over role of other genes is provided in Supplementary. Here we focus on
structured perturbation. Using the F-measure criterion (Mall, Langone, and
Suykens 2013), we detected optimal $k=3$ at F-measure$=$0.30 for $\mathbf{Z}$.
Similarly, we identified $3$ structured perturbations in the
$\hat{\mathbf{Z}}$ using STEPD method at F-measure$=$0.30 (for
$\nu=0.5$,$\mu=0.5$). We attained a precision of $1.0$ and recall of $0.952$
when comparing the edges in the structured perturbations of $\mathbf{Z}$ and
$\hat{\mathbf{Z}}$. The largest cluster comprised of $286$ genes (TFs), while
the $2$ smaller clusters consisted of $10$ TFs each. We showcase the original
sub-networks of the IDH-mutant and IDH-wildtype graphs corresponding to the
TFs in each cluster in Figure 4. Additional caveats are provided in
Supplementary.

To carry out a deeper biological analysis, we use the GO Onotology (GO terms)
which comprises a repository of known functions and processes associated with
genes. For example, each gene can be classified as participating in a
Biological Process (BP), Cellular component (CC), Molecular Function (MF) or
Pathways (PW). Each of these high level processes and functions consist of
terms, for example, “defense response to virus” or “protein modification
process.” For each perturbation cluster, we identified the enriched or over-
represented GO terms using ConsensusPathDB (Kamburov et al. 2012).

A detailed breakdown of the enriched GO terms for both the IDH-mutant
($\mathbf{X}$) and IDH-wildtype ($\mathbf{Y}$) case is provided in
Supplementary. It was shown in (Ceccarelli, Barthel et al. 2016) that the main
difference between IDH-mutant and IDH-wildtype gliomas is the characteristic
hyper-methylation phenotype i.e. chromatin modification and histone
acetylation. Furthermore, it was shown in (Frattini et al. 2018) that
PPAR$\alpha$ is recruited for mitochondrial respiration leading to tumor
oncogenesis in fusion gliomas. This is indicated by the over-represented
biological components and pathways for IDH-mutant and IDH-wildtype case
respectively for perturbation cluster $C_{2}$. Details of the enriched
biological components and pathways associated with perturbation clusters
$C_{1}$ and $C_{2}$ are provided in the Supplementary. The STEPD approach is
not only able to identify known but also detects potential novel enrichments
which need further investigation.

## Conclusion

In this paper, we introduced the problem of Permutation and Structured
Perturbation Inference (PSPI), a combinatorial problem with real world
applications in the fields of system biology and computer vision. PSPI is a
generalization of the graph matching problem and takes structured
perturbations into account while computing the permutations between a pair of
graphs. We proposed STEPD, an iterative algorithm to solve the relaxed version
of the problem. Experiments on computational biology benchmark datasets showed
that PSPI can be used to make potentially insightful discoveries from gene
regulatory networks and can become an important prognostic tool for
biologists. Source code and datasets for reproducibility can be found at
https://github.com/code-halo/rgm/.

## References

  * Borgwardt and Kriegel (2005) Borgwardt, K.; and Kriegel, H. 2005.  Shortest-path kernels on graphs.  In _ICML_.
  * Boyd et al. (2011) Boyd, S.; Parikh, N.; Chu, E.; Peleato, B.; and Eckstein, J. 2011.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  _Foundations and Trends® in Machine Learning_ 3\.
  * Candès et al. (2011) Candès, E. J.; Li, X.; Ma, Y.; and Wright, J. 2011.  Robust principal component analysis?  _Journal of the ACM (JACM)_ 58\.
  * Ceccarelli, Barthel et al. (2016) Ceccarelli, M.; Barthel, F. P.; et al. 2016.  Molecular Profiling Reveals Biologically Discrete Subsets and Pathways of Progression in Diffuse Glioma.  _Cell_ 164\.
  * Cho, Lee, and Lee (2010) Cho, M.; Lee, J.; and Lee, K. M. 2010.  Reweighted Random Walks for Graph Matching.  In _ECCV_ , 492–505. Springer.
  * Conte et al. (2004) Conte, D.; Foggia, P.; Sansone, C.; and Vento, M. 2004.  Thirty years of graph matching in pattern recognition.  _International journal of pattern recognition and artificial intelligence_ 18\.
  * Demarin et al. (2016) Demarin, V.; Roje Bedeković, M.; Bosnar Puretić, M.; and Bošnjak Pašić, M. 2016.  Arts, brain and cognition.  _Psychiatria Danubina_ 28(4): 343–348.
  * Dym, Maron, and Lipman (2017) Dym, N.; Maron, H.; and Lipman, Y. 2017.  DS++: A Flexible, Scalable and Provably Tight Relaxation for Matching Problems.  _ACM Transaction on Graph_ 36(6).
  * Fiori et al. (2013) Fiori, M.; Sprechmann, P.; Vogelstein, J.; Musé, P.; and Sapiro, G. 2013.  Robust multimodal graph matching: Sparse coding meets graph matching.  In _NeurIPS_.
  * Frattini et al. (2018) Frattini, V.; Pagnotta, S. M.; Fan, J. J.; Russo, M. V.; Lee, S. B.; Garofano; et al. 2018.  A metabolic function of FGFR3-TACC3 gene fusions in cancer.  _Nature_ .
  * Grover and Leskovec (2016) Grover, A.; and Leskovec, J. 2016.  Node2vec: Scalable Feature Learning for Networks.  In _KDD_ , 855–864. New York, NY, USA.
  * Hamilton, Ying, and Leskovec (2017) Hamilton, W. L.; Ying, R.; and Leskovec, J. 2017.  Inductive Representation Learning on Large Graphs.  In _NeurIPS_ , 1025–1035.
  * Jia et al. (2016) Jia, K.; Chan, T.-H.; Zeng, Z.; Gao, S.; Wang, G.; Zhang, T.; and Ma, Y. 2016.  ROML: A Robust Feature Correspondence Approach for Matching Objects in A Set of Images.  _Interntaional Journal of Computer Vision_ 117: 173–197.
  * Johnson, Li, and Rabinovic (2007) Johnson, W. E.; Li, C.; and Rabinovic, A. 2007.  Adjusting batch effects in microarray expression data using empirical Bayes methods.  _Biostatistics_ 8\.
  * Kamburov et al. (2012) Kamburov, A.; Stelzl, U.; Lehrach, H.; and Herwig, R. 2012.  The ConsensusPathDB interaction database: 2013 update.  _Nucleic acids research_ 41\.
  * Kashima, Tsuda, and Inokuchi (2003) Kashima, H.; Tsuda, K. T.; and Inokuchi, A. 2003.  Marginalized kernels between labeled graphs.  In _ICML_.
  * Kipf and Welling (2016) Kipf, T. N.; and Welling, M. 2016.  Semi-Supervised Classification with Graph Convolutional Networks.
  * Lee, Rossi, and Kong (2018) Lee, J. B.; Rossi, R.; and Kong, X. K. 2018.  Graph Classification using Structural Attention.  In _KDD_.
  * Lin, Liu, and Su (2011) Lin, Z.; Liu, R.; and Su, Z. 2011.  Linearized alternating direction method with adaptive penalty for low-rank representation.  In _NeurIPS_.
  * Luenberger and Ye (2015) Luenberger, D. G.; and Ye, Y. 2015.  _Linear and nonlinear programming_ , volume 228.  Springer.
  * Lê-Huu and Paragios (2017) Lê-Huu, D. K.; and Paragios, N. 2017.  Alternating Direction Graph Matching.  In _CVPR_ , 4914–4922.
  * Mall et al. (2017) Mall, R.; Cerulo, L.; Bensmail, H.; Iavarone, A.; and Ceccarelli, M. 2017.  Detection of statistically significant network changes in complex biological networks.  _BMC Systems Biology_ 11\.
  * Mall et al. (2018) Mall, R.; Cerulo, L.; Garofano, L.; Frattini, V.; Kunji, K.; Bensmail, H.; Sabedot, T. S.; Noushmehr, H.; Lasorella, A.; Iavarone, A.; et al. 2018.  RGBM: regularized gradient boosting machines for identification of the transcriptional regulators of discrete glioma subtypes.  _Nucleic Acids Research_ .
  * Mall, Langone, and Suykens (2013) Mall, R.; Langone, R.; and Suykens, J. A. 2013.  Self-tuned kernel spectral clustering for large scale networks.  In _Big Data, 2013 IEEE International Conference on_. IEEE.
  * Ming et al. (2014) Ming, F.; Wong, F.; Liu, Z.; and Chiang, M. 2014.  Stock market prediction from WSJ: text mining via sparse matrix factorization.  In _ICDM_. IEEE.
  * Perozzi, Al-Rfou, and Skiena (2014) Perozzi, B.; Al-Rfou, R.; and Skiena, S. 2014.  DeepWalk: Online Learning of Social Representations.  In _KDD_ , 701–710.
  * Przulj (2006) Przulj, N. 2006.  Biological network comparison using graphlet degree distribution.  In _In proceedings of European Conference on Computational Biology_.
  * Shervashidze and M.Borgwardt (2006) Shervashidze, N.; and M.Borgwardt, K. M. 2006.  Fast Subtree Kernels on Graphs.  In _In proceedings of European Conference on Computational Biology_.
  * Singh, Xu, and Berger (2007) Singh, R.; Xu, J.; and Berger, B. 2007.  Pairwise global alignment of protein interaction networks by matching neighborhood topology.  In _Annual International Conference on Research in Computational Molecular Biology_. Springer.
  * Sinkhorn and Knopp (1967) Sinkhorn, R.; and Knopp, P. 1967.  Concerning nonnegative matrices and doubly stochastic matrices.  _Pacific Journal of Mathematics_ 21\.
  * Umeyama (1988) Umeyama, S. 1988.  An eigendecomposition approach to weighted graph matching problems.  _IEEE TPAMI_ 10\.
  * Vogelstein et al. (2015) Vogelstein, J. T.; Conroy, J. M.; Lyzinski, V.; Podrazik, L. J.; Kratzer, S. G.; Harley, E. T.; Fishkind, D. E.; Vogelstein, R. J.; and Priebe, C. E. 2015\.  Fast approximate quadratic programming for graph matching.  _PLOS one_ 10\.
  * Wallace, Martin, and Ambs (2011) Wallace, T.; Martin, D.; and Ambs, S. 2011.  Interaction among genes, tumor biology and the environment in cancer health disparities: examining the evidence on a national and global scale.  _Carcinogenesis_ 32\.
  * Wang, Li, and Konig (2010) Wang, F.; Li, P.; and Konig, A. C. 2010.  Learning a bi-stochastic data similarity matrix.  In _ICDM_. IEEE.
  * Wang et al. (2018) Wang, T.; Ling, H.; Lang, C.; and Feng, S. 2018.  Graph Matching with Adaptive and Branching Path Following.  _IEEE TPAMI_ 40(12): 2853–2867.
  * White et al. (1986) White, J. G.; Southgate, E.; Thomson, J. N.; and Brenner, S. 1986.  The structure of the nervous system of the nematode Caenorhabditis elegans.  _Philos Trans R Soc Lond B Biol Sci_ 314\.
  * Yan et al. (2015) Yan, J.; Xu, H.; Zha, H.; Yang, X.; Liu, H.; and Chu, S. 2015.  A Matrix Decomposition Perspective to Multiple Graph Matching.  In _ICCV_ , 199–207.
  * Yanardag and Vishwanathan. (2015) Yanardag, P.; and Vishwanathan., S. V. N. 2015.  Deep Graph Kernels.  In _KDD_.
  * Yuan and Lin (2006) Yuan, M.; and Lin, Y. 2006.  Model selection and estimation in regression with grouped variables.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 68\.
  * Zaslavskiy, Bach, and Vert (2009) Zaslavskiy, M.; Bach, F.; and Vert, J.-P. 2009.  A path following algorithm for the graph matching problem.  _IEEE TPAMI_ 31\.
  * Zhang, Horvath et al. (2005) Zhang, B.; Horvath, S.; et al. 2005.  A general framework for weighted gene co-expression network analysis.  _Statistical applications in genetics and molecular biology_ 4\.