Joint alignment of multiple protein-protein interaction networks via
convex optimization

arXiv:1604.03482v1 [cs.DS] 12 Apr 2016

Hashemifar, Somaye
hashemifar@ttic.edu

Huang, Qixing
huangqx@ttic.edu

Xu, Jinbo
j3xu@ttic.edu

Abstract.
Motivation: High-throughput experimental techniques have been producing more and more protein-protein
interaction (PPI) data. PPI network alignment greatly benefits the understanding of evolutionary relationship
among species, helps identify conserved sub-networks and provides extra information for functional annotations. Although a few methods have been developed for multiple PPI network alignment, the alignment
quality is still far away from perfect and thus, new network alignment methods are needed.
Result: In this paper, we present a novel method, denoted as ConvexAlign, for joint alignment of multiple PPI networks by convex optimization of a scoring function composed of sequence similarity, topological
score and interaction conservation score. In contrast to existing methods that generate multiple alignments in
a greedy or progressive manner, our convex method optimizes alignments globally and enforces consistency
among all pairwise alignments, resulting in much better alignment quality. Tested on both synthetic and
real data, our experimental results show that ConvexAlign outperforms several popular methods in producing functionally coherent alignments. ConvexAlign even has a larger advantage over the others in aligning
real PPI networks. ConvexAlign also finds a few conserved complexes among 5 species which cannot be
detected by the other methods.

1

Introduction

Protein-protein interaction (PPI) networks provide valuable information for understanding of protein functions and system-level cellular processes. The alignment of PPI networks is a useful means for comparing the
networks of different species. This comparison helps identify evolutionarily conserved pathways/complexes
that may be functionally significant. Studying the conserved modules may provide useful information about
the molecular mechanism contributing to their functions.
PPI networks can be aligned either locally or globally. Local network alignment methods such as Mawish [12] and AlignNemo [2] aim to find small isomorphic subnetworks. Global network alignment (GNA)
methods maximize the overall match between input networks. Some GNA methods such as IsoRank [26,27],
MI-GRAAL [13], GHOST [19], MAGNA [23, 32], Prob [30], NETAL [18] and HubAlign [5] are designed
for pairwise alignment, while others such as IsoRankN [15] and NetCoffee [6] for multiple alignment. GNA
can be one-to-one or many-to-many mapping. The latter allows one protein to be aligned to multiple proteins
of a single network while the former does not.
More attention has been paid to pairwise network alignment. With the availability of more PPI networks,
it becomes inevitable to align multiple networks. Existing GNA methods such as NetworkBlast-M [9, 25]
and GraemLin 2.0 [3] are designed for local alignment of multiple networks, whereas others such as IsoRankN [15], SMETANA [22], NetCoffee [6], BEAMS [1] and FUSE [4] for global alignment of multiple
networks. In addition to sequence similarity, all these methods excluding NetworkBlast-M and NetCoffee also employs topological information. Moreover, all the methods except NetCoffee are designed for
1

many-to-many alignments. NetworkBlast-M starts with a set of highly conserved regions and then extends
them greedily. GraemLin2.0 integrates phylogenetic information and network topology and then employs
a hill-climbing algorithm to generate the alignment. IsoRankN applies IsoRank to compute the alignment
scores between each pair of networks and then uses a PageRank-Nibble algorithm to cluster the proteins.
SMETANA employs a semi-Markov random walk model to measure similarity between proteins. BEAMS
constructs a weighted k-partite graph in which edges are assigned weights derived from protein sequence
similarity. NetCoffee applies a triplet approach similar to T-Coffee to compute the edge weights of the kpartite graph. Both BEAMS and NetCoffee apply a heuristic on the k-partite graph to build an alignment.
BEAMS fulfills this by greedily merging a set of disjoint cliques while NetCoffee by applying a simulated
annealing method on a set of candidates. FUSE applies a non-negative matrix tri-factorization method to
compute edge weights of the k-partite graph.
Most of existing GNA methods do not optimize alignment of all proteins simultaneously. Instead, they
start from the best alignment between a subset of proteins and then gradually extend it by adding more
proteins using a greedy strategy. This may impact alignment quality since errors introduced at an earlier
stage cannot be fixed later.
This paper presents a novel one-to-one GNA algorithm, denoted as ConvexAlign, to align multiple
PPI networks using a new scoring scheme that integrates network topology, sequence similarity and interaction conservation score. It is NP-hard to optimize such a scoring function. We formulate this GNA
problem as an integer program and relax it to a convex optimization problem, which enables us to simultaneously align all the PPI networks, without resorting to the widely-used seed-and-extension or progressive
alignment methods. Then we use an ADMM (alternating direction method of multipliers) method (see
http://stanford.edu/˜boyd/admm.html) to solve the relaxed convex optimization problem and
optimize all the protein mappings together. Tested on the PPI networks of five different species, ConvexAlign outperforms several popular methods such as IsoRankN, SMETANA, NetCoffee and BEAMS in
terms of biological alignment quality. ConvexAlign finds a few conserved complexes among these 5 species
which cannot be found by the other methods. ConvexAlign also performs very well in aligning some publicly available synthetic networks.

2

Method

Definition. We represent a protein-protein interaction network by an undirected graph G = (V, E) where
V is the set of vertices (proteins) and E the set of edges (interactions). Let d(u) denote the degree of
vertex u and e = (u, v) ∈ E represent an edge. A one-to-one global alignment A between N networks
Gi = (Vi , Ei ), 1 ≤ i ≤ N , is given by a decomposition of all nodes V = ∪N
i=1 Vi such that V = A1 ∪· · ·∪AK
where each Ai contains at most one protein from each network and any two Ai and Aj are disjoint. We call
each Ai in the alignment a group or a cluster. Proteins in each cluster are mutually aligned to one another.

2.1

Scoring function for network alignment

Our goal is to find an alignment that maximizes the number of preserved edges and the number of matched
orthologous (or functionally conserved) proteins. For this purpose we use a node score for scoring matched
proteins and an edge score for scoring matched interactions, respectively. For a pair of proteins, their node
score is the combination of their topology score and sequence similarity score. We use a minimum-degree
heuristic algorithm to calculate the topological score, which was used by us to develop a pairwise GNA
method HubAlign [5]. A recent third-party evaluation by Prz̆ulj group [16] has shown that this topological
score works very well in pairwise GNA. Please see our paper [5] for more details. We use the normalized
BLAST bit scores for sequence similarity. Let B(vi , vj ) and T (vi , vj ) respectively denote the sequence sim-

2

ilarity and topology score between a pair of proteins vi ∈ Vi and vj ∈ Vj . Then the node score node(vi , vj )
is calculated as follows:
node(vi , vj ) = (1 − λ1 )B(vi , vj ) + λ1 T(vi , vj ),
(1)
where λ1 controls the importance of the topology score relative to the BLAST score. The node score of
multiple alignment A, i.e. fnode (A), sums the scores among all pairs of matched proteins:
X
X
node(vi , vj ).
(2)
fnode (A) =
1≤i<j≤N Ak ∈A,vi ,vj ∈Ak

The edge score finteraction (A) measures interaction-preserving in an alignment A. This score counts the
number of interactions aligned between all pairs of networks:
X
X
finteraction (A) =
δ((vi , vi0 ) ∈ Ei )δ((vj , vj0 ) ∈ Ej ),
(3)
1≤i<j≤N Ak ,Al ∈A,vi ,vj ∈Ak ,vi0 ,vj0 ∈Al

where δ((vi , vi0 ) ∈ Ei ) is an indicator function. We aim to find the multiple alignment A that maximizes a
combination of node and interaction scores as follows.
f = (1 − λ2 )fnode (A) + λ2 finteraction (A),

(4)

where λ2 describes the tradeoff. See Appendix for determination of λ1 and λ2 by cross-validation.

2.2

Integer and Convex Programming Formulation

Definition. A one-to-one multiple network alignment is valid or feasible if the following condition (also
called consistency property) is satisfied: for any three vertices vi , vj , and vk of three different networks, if
vi is aligned to vj and vj aligned to vk , then vi is aligned to vk .
Parameterizing multiple alignments. Let M be the number of proteins in all the input PPI networks i.e.
N
P
M=
|Vi |. We may represent a valid multiple alignment A by a binary matrix Y = (Y1 ; Y2 ; · · · ; YN ) ∈
i=1

{0, 1}M ×K , where each block Yi encodes the association between Vi and A. Each row of Y corresponds
to one vertex and each column to one alignment cluster. That is, ∀vi ∈ Vi , Yi (vi , Aj ) = 1 if and only if vi
is in cluster Aj . Let 1 be a vector of appropriate size with all elements being 1. Since Y is a one-to-one
alignment, it shall satisfy the following constraints:
• Each row of Y has exactly one non-zero entry, i.e., Y 1 = 1 .
• Each column of Y has at most N non-zero entries, i.e., Y T 1 ≤ N 1.
• Each column of Yi has at most one non-zero entry, i.e., YiT 1 ≤ 1.
On the other direction, any binary matrix Y satisfying the above properties encodes a one-to-one alignment.
Although Y is a good representation of an MNA, the objective function with Y as variable is nonlinear
and thus, hard to optimize. Inspired by [7], we introduce another alignment matrix X as follows.
 


I|V1 | X12 · · · X1N
Y1
 X T I|V | · · · X2N   Y2 

2
 12
 

T
T
T
(5)
X= .
..  =  ..  · Y1 Y2 · · · YN ,
.
.
.
 .
.
···
.   . 
T
X1N

···

· · · I|VN |

YN

3

where each block Xij = Yi YjT is a binary matrix encoding the mapping between Vi and Vj . That is,
Xi,j (vi , vj ) = 1 if and only if vi and vj are aligned (i.e., in the same alignment cluster).
It is easy to see that X is positive semi-definite. Since this section considers only one-to-one mapping,
for any two i and j (i 6= j), each row or column of Xij has at most one non-zero element, i.e., Xij 1 ≤ 1
and XijT 1 ≤ 1 where 1 is a vector of appropriate size with all entries 1. On the other direction, we have the
following proposition (see Appendix D for its proof).
Proposition 1. Let X be a binary block matrix with N × N blocks, and Xij be the block in the ith row
and the j th column. If X satisfies the following conditions: (1) X  0, (2)Xii = I|Vi | for 1 ≤ i ≤ N , and
(3)Xij 1 ≤ 1 and XijT 1 ≤ 1 for 1 ≤ i < j ≤ N , then X encodes a feasible global alignment of N networks
admitting one-to-one mapping and satisfying the cycle consistency property.
Therefore, we may encode a one-to-one GNA using X, which leads to a linear formulation of the
objective function. Following Prop. 1, we impose the following constraints on X:
Xij 1 ≤ 1, XijT 1 ≤ 1, Xij ∈ {0, 1}|Vi |×|Vj |
X  0,

Xii = I|Vi |

(1 ≤ i < j ≤ N )
(1 ≤ i ≤ N )

(6)

Objective function. As Xij is the indicator submatrix for Vi and Vj , the node score can be formulated as
follows.
X
X
X
node(v, v 0 )Xij (v, v 0 ) =
hCij , Xij i,
(7)
fnode =
1≤i<j≤N v∈Vi ,v 0 ∈Vj

1≤i<j≤N

where Cij is a matrix composed of the values of node(v, v 0 ).
To formulate finteraction , we introduce indicator variables yij (vi , vj , vi0 , vj0 ) for edge correspondences:
yij (vi , vj , vi0 , vj0 ) = Xij (vi , vj )Xij (vi0 , vj0 ),
finteraction =

X

∀(vi , vi0 ) ∈ Ei , (vj , vj0 ) ∈ Ej , 1 ≤ i < j ≤ N.

X

yij (vi , vj , vi0 , vj0 ) =

1≤i<j≤N (vi ,vi0 )∈Ei ,(vj ,vj0 )∈Ej

X

h1, y ij i,

(8)
(9)

1≤i<j≤N

where y ij stacks the indicator variables between Vi and Vj .
The nonlinear constraint (8) can be replaced by the following linear inequalities (c.f. [8, 14]):
X
∀vj0 ∈ Vj ,
y(vi , vj , vi0 , vj0 ) ≤ Xij (vi , vj )
vi0 :(vi ,vi0 )∈Ei

X

∀vi0 ∈ Vi ,

y(vi , vj , vi0 , vj0 ) ≤ Xij (vi , vj )

vj0 :(vj ,vj0 )∈Ej

∀vj ∈ Vj ,

X

y(vi , vj , vi0 , vj0 ) ≤ Xij (vi0 , vj0 )

vi :(vi ,vi0 )∈Ei

∀vi ∈ Vi ,

X

y(vi , vj , vi0 , vj0 ) ≤ Xij (vi0 , vj0 )

(10)

vj :(vj ,vj0 )∈Ej

It is easy to prove that (8) implies (10). On the other direction, considering that the coefficients of y is
positive and we want to maximize (9), we shall be able to prove that (10) implies (8). We replace (8) by (10)
to obtain linear constraints and summarize (10) in the matrix form as follows.
Bij y ij ≤ Fij (Xij ),
4

(11)

where Bij is coefficient and Fij is a linear operator that picks the corresponding element of Xij for each
constraint. That is, Fij (Xij (vi , vj )) =< Pij , Xij > where Pij is a binary matrix with the same dimension
as Xij and only one element Pij (vi , vj ) is equal to 1.
Finally, by integrating (7), (9), (11) and Prop. 1, we have the following integer program:

X 
maximize
(1 − λ2 )hCij , Xij i + λ2 h1, yij i
1≤i<j≤N

subject to y ij ∈ {0, 1}|Ei |×|Ej | , Bij y ij ≤ Fij (Xij ),
Xij 1 ≤ 1, XijT 1 ≤ 1, Xij ∈ {0, 1}|Vi |×|Vj | ,
X  0,

Xii = I|Vi | ,

1≤i<j≤N
1≤i<j≤N

1≤i≤N

(12)

The key constraint is X  0, which enforces cycle consistency in the alignments. X  0 still holds even
N = 2.

2.3

Optimization via Convex Relaxation

It is NP-hard to directly optimizing (12) since the variables are binary. We may first relax them to obtain a
convex optimization problem that can be solved to global optimum within polynomial time, and then employ
a greedy rounding scheme to convert fractional solution to integral.
Convex relaxation. By relaxing y ij and Xij to real values between 0 and 1, we have the following convex
program:

X 
maximize
(1 − λ2 )hCij , Xij i + λ2 h1, yij i
1≤i<j≤N

subject to y ij ≥ 0, Bij y ij ≤ Fij (X),

1≤i<j≤N

Xij 1 ≤ 1, XijT 1 ≤ 1, Xij ≥ 0,
X  0,

1≤i<j≤N
1≤i≤N

Xii = I|Vi | ,

(13)

Optimization strategy. We use ADMM (alternating direction of multiplier method) to solve the convex
relaxation (13). The basic idea is to augment its Lagrangian dual(see https://en.wikipedia.org/
wiki/Augmented_Lagrangian_method) and iteratively optimize a subset of variables while keeping the others fixed. This allows us to exploit structure patterns in the constraint set for effective optimization. As the derivation is quite technical, we leave the details in Appendix E.
Rounding into an integer solution. The above convex relaxation has a pretty tight fractional solution. We
propose a greedy rounding strategy to convert fractional solution to integral. We collect all the protein pairs
with an indicator value X(u, v) > 0.05 and place them in a decreasing order into a sorted list X . Then we
build an alignment graph starting with an empty edge set by scanning through X . For each scanned protein
pair (u, v) in X , in the alignment graph we add an edge to connect this pair as long as such an addition does
not violate the constraint that no protein in one network is aligned to two proteins in another network. After
all pairs are scanned, we decompose the alignment graph into connected components, each corresponding
to a cluster of mutually-aligned proteins. The set of all the clusters form an alignment. Most components
are cliques. For the very few non-clique components, we just add some edges to make them cliques.

3

RESULTS

We compare our algorithm ConvexAlign with several popular and publicly available methods IsoRankN [15],
SMETANA [22], NetCoffee [6] and BEAMS [1]. We ran SMETANA and NetCofee with their default pa5

rameters. For both BEAMS and IsoRankN, we set three different values for their parameter α = {0.3, 0.5, 0.7}.
We left other parameters of BEAMS at their default values.

3.1

Test data

We use the PPI networks of H.sapiens (human), S.cerevisiae (yeast), Drosophila melanogaster (fly), Caenorhabditis elegans (worm) and Mus musculus (mouse) taken from IntAct [10]. The human network has 9003 proteins and 34935 interactions, the yeast network has 5674 proteins and 49830 interactions, the fly networks
has 8374 nodes and 25611 interactions, the mouse network has 2897 proteins and 4372 interactions and the
worm network has 4305 proteins and 7747 interactions. Only experimentally-validated PPIs are used.
We also use the NAPAbench [21] synthetic PPI networks. NAPAbench is a benchmark that contains
PPI network families generated by three different network models: crystal growth(CG) [11], duplicationmutation-complementation(DMC) [31] and duplication with random mutation(DMR) [28]. We use the 8way alignment dataset of this benchmark, which contains three network families each with 8 networks
of 1000 nodes generated by one of the three network models. The 8-way alignment dataset simulates a
network family of closely-related species, so this benchmark has very different properties as the above 5
real PPI networks. NAPAbench has recently been used to benchmark SMETANA.

3.2

Alignment quality measures

We evaluate multiple network alignment quality using several topological and functional consistency metrics proposed in different studies. Functional consistency measures however, are more important than topological measures since one of the important applications of network alignment is to functional annotation
transfer. For topological analysis of the output clusters we use the following metrics.
c-coverage: It is the number of clusters composed of proteins from exactly c species. Specifically, total
coverage is the number of clusters composed of proteins from at least two species. Clusters with large c
explain a larger amount of data better than clusters with small c.
Conserved Interaction(CI): It is calculated as the ratio of the number of aligned interactions to the total
number of interactions between output clusters.
A multiple alignment with a higher c-coverage (or total coverage) or CI is not necessarily biologically
meaningful since it may align many unrelated proteins together. Therefore, we also employ GO terms to
measure functional consistency or biological quality of an alignment. GO terms describe roles of proteins
in terms of their associated biological process (BP), molecular function (MF) and cellular component (CC).
We exclude root GO terms from analysis, i.e., GO terms on level higher than 5. We also exclude CC because
proteins with matched CC are not usually considered functionally similar. Moreover, CC only annotates a
small percentage of the proteins. The following measures are based on the observation that functionally
related proteins are more likely to have similar GO terms.
Specificity: We say a cluster annotated if at least two of its proteins have GO annotations. An annotated
cluster is consistent if all of its proteins share at least one common GO term. Specificity is defined as the
ratio of consistent clusters to annotated clusters.
Average of functional similarity (AF S). This score is based on the semantic similarity of the GO terms,
which is derived from their distance in the ontology. We use Schlicker's similarity, based on the Resnik
ontological similarity, to calculate the functional similarity in the BP and MF category [24]. Schlicker's
similarity is one of the best performing methods for computing the functional similarity between proteins
(see Appendix B for more details) [20]. Let scat (u, v) denote the GO functional similarity of proteins u and
v in category cat (i.e. BP or MF). AF S of an output cluster A in category cat is defined as follows:
X
1
AF Scat (A) = |A|×(|A|−1)
scat (vi , vj )
vi ,vj ∈A,i6=j

2

6

Finally we define AF S in category cat as the average of AF S over all clusters. Following [6], we take
into consideration all the clusters that contain at least 60% GO-annotated proteins to avoid ignoring many
functionally meaningful clusters. We separately compare the AF S for clusters for c = 3, 4, 5. We also
provide the distribution of the AFS scores for each given c.
Mean normalized
entropy (M N E). The normalized entropy of a cluster A is defined as: N E(A) =
Pd
1
p
×
log(pi ) where d is the number of different GO annotations in A and pi represents the
×
i=1 i
log(d)
fraction of proteins in A with annotation GOi . A cluster with lower entropy is more functionally coherent.
M N E is the mean of normalized entropy over all annotated clusters.
Conserved orthologous interactions(COI): Similar to SMETANA, COI is calculated as the total number
of interactions between all consistent clusters. COI may be a better measure than CI because it detects
whether the conserved interactions are spurious or actually correspond to real conserved interactions between orthologous proteins. An alignment with larger COI may lead to identifying functionally conserved
subnetwroks (i.e clusters) composed of orthologous genes.
Sensitivity: The closest cluster of a given GO term is the cluster that contains the maximum number of
proteins associated with that GO term. Similar to BEAMS [1], we define Sensitivity as the average (over all
GO terms) of the fraction of proteins in the closest cluster, that are associated with that GO term.

3.3

Alignment quality on real data

Topological quality. Table 1 lists the topological evaluation of the alignments produced by different
methods. The first four multi-rows show the results for the clusters consisting of proteins belonging to
c = 2, 3, 4, 5 species, respectively. In each multi-row, the top and bottom rows show c-coverage and the
number of proteins in the clusters, respectively. ConvexAlign has a larger c-coverage when c = 4, 5 than the
other methods except SMETANA and NetCoffee. However, as we show later, many of clusters generated by
these two methods are not functionally conserved. The total coverage of BEAMS and IsoRank is better than
the others because they produce many clusters composed of proteins from 2 or 3 species. These clusters
can not explain the data as well as clusters containing proteins from 4 or 5 species can. ConvexAlign has
a better CI than all other methods except SMETANA. These conserved interactions may be very helpful in
identifying the functional modules conserved among networks of different species. It is worth mentioning
that most of the conserved interaction resulting from SMETANA may be spurious [1].
Table 1: Topological evaluation of output clusters by different alignment methods. IsoRankN and BEAMS
are tested using three different values of their parameters α.
IsoRankN IsoRankN IsoRankN SMETANA NetCoffee BEAMS BEAMS BEAMS ConvexAlign
(0.3)
(0.5)
(0.7)
(0.3)
(0.5)
(0.7)
c=2
4625
4178
4670
1127
1424
5703
5274
5271
2856
11035
8356
11165
2718
2848
11406
11469
11465
5712
c=3
2259
2270
2304
1653
1739
2192
2557
2556
1833
8521
6810
8750
5808
5217
6576
8128
8118
5499
c=4
1023
731
944
2028
1980
1163
1141
1143
1190
5276
2924
4823
9531
7920
4652
4686
4701
4760
c=5
224
112
184
1622
1217
683
600
600
765
1417
560
1182
10342
6075
3915
3046
3044
3825
Total
8131
7291
8102
6430
6360
9741
9572
9570
6644
coverage 26249
18650
25920
28399
22070
26549
27329
27328
19796
CI
0.03
0.02
0.03
0.10
0.03
0.03
0.03
0.03
0.04
CIQ
0.03
0.02
0.02
0.06
0.02
0.02
0.02
0.02
0.03

Biological quality. Table 2 provides the functional consistency measures of the alignments generated by
7

different methods. The first four multi-rows show the quality of the clusters composed of proteins from c =
2, 3, 4, 5 species. In these multi-rows, the top and middle rows show the number of consistent and annotated
clusters, respectively, and the bottom row shows specificity. Regardless of c, ConvexAlign outperforms the
other methods in terms of specificity and the number of consistent clusters. At the same time, ConvexAlign
generates fewer annotated clusters than BEAMS when c = 2, 3, 4. , Although SMETANA and NetCoffee
generate a larger number of clusters for c = 4, 5 than ConvexAlign, their clusters are not very functionally
consistent. The fifth row shows ConvexAlign has much higher specificity than the others when all the
resulting clusters (c ≥ 2) are considered. These results suggest that ConvexAlign finds more functionally
consistent clusters, not only by generating small clusters (i.e. c = 2, 3) but more importantly large clusters
(i.e. c = 4, 5). These clusters (especially when c = 4, 5) are very valuable because they may provide useful
information about the orthology relationship among the proteins of all species. Moreover, these clusters
can be very useful for identifying conserved sub-networks as well as predicting the function of unannotated
proteins. ConvexAlign yields a COI/CI ratio around 60% that is 1.44 times larger than the second best ratio
by BEAMS. This result may indicate that ConvexAlign is able to identify conserved interactions between
othologous proteins. It also suggests that although SMETANA has the largest CI, many of those conserved
interactions are possibly false and formed by non-orthologous proteins. ConvexAlign also outperforms other
methods in terms of M N E and sensitivity.
Table 2: Functional consistency of output clusters. Note that for MNE, the smaller the better; while for the
other measures, the larger the better.

consistent
annotated
specifity
consistent
c=3 annotated
specifity
consistent
c=4 annotated
specifity
consistent
c=5 annotated
specifity
c ≥ 2 specifity
COI
COI/CI
MNE
Sensitivity
c=2

IsoRankN IsoRankN IsoRankN SMETANA NetCoffee BEAMS BEAMS BEAMS ConvexAlign
(0.3)
(0.5)
(0.7)
(0.3)
(0.5)
(0.7)
906
1259
919
295
495
1539
1568
1569
1914
3614
2862
3646
931
931
3486
3456
3452
2326
0.25
0.44
0.25
0.39
0.53
0.44
0.45
0.45
0.82
203
466
231
188
462
1003
1084
1084
1155
2160
2153
2210
1556
1640
2119
2442
2441
1741
0.09
0.22
0.10
0.12
0.28
0.47
0.44
0.44
0.66
41
106
54
170
406
606
624
624
661
1020
723
942
2019
1640
1159
1136
1138
1079
0.04
0.15
0.06
0.08
0.25
0.52
0.55
0.55
0.61
14
19
9
183
406
383
359
359
493
224
112
184
1621
1955
683
600
600
763
0.06
0.17
0.05
0.11
0.21
0.56
0.60
0.60
0.65
0.17
0.32
0.17
0.14
0.29
0.48
0.48
0.48
0.71
88
188
127
480
553
1237
1311
1305
1668
0.02
0.13
0.03
0.04
0.21
0.40
0.41
0.41
0.59
2.15
2.19
2.14
2.44
2.39
1.97
1.95
1.95
1.93
0.45
0.46
0.45
0.36
0.22
0.33
0.31
0.37
0.51

Table 3 shows the AF S separately for clusters composed of proteins in 3, 4, and 5 species in both
categories BP and MF. The AF S obtained by ConvexAlign is 6 − 20% larger than the other methods. These
results indicate that on average the clusters generated by ConvexAlign are functionally more consistent. That
is, ConvexAlign outperforms the other methods in terms of not only the number of consistent clusters, but
also the average GO semantic similarity. The distribution of AFS scores for clusters composed of proteins
in 3, 4, and 5 species is shown in Fig. A.1, in which the middle line in each box shows the median value.
That is, the median AFS obtained by ConvexAlign is larger than the other methods. These results further
confirm that ConvexAlign yields clusters with higher functional similarity in both categories MF and BP.

8

Table 3: AF S comparison between ConvexAlign and the other methods

c=3
AF S BP c=4
c=5
c=3
AF S M F c=4
c=5

3.4

IsoRankN IsoRankN IsoRankN SMETANA NetCoffee BEAMS BEAMS BEAMS ConvexAlign
(0.3)
(0.5)
(0.7)
(0.3)
(0.5)
(0.7)
0.83
1.02
0.86
0.74
1.03
1.60
1.63
1.63
1.74
0.69
0.97
0.72
0.68
0.99
1.63
1.61
1.60
1.79
0.75
1.01
0.72
0.85
1.16
1.66
1.67
1.67
1.71
0.80
0.94
0.80
0.69
0.99
1.40
1.33
1.34
1.54
0.83
1.02
0.86
0.74
1.03
1.60
1.63
1.63
1.74
0.86
1.06
0.86
0.94
1.18
1.68
1.68
1.68
1.74

Alignment quality on synthetic data

This section explains the results on the NAPAbench benchmark. Fig. 1 shows the number of consistent
clusters generated by different methods and their specificity on clusters composed of proteins in c =
2, 3, 4, 5, 6, 7, 8 species, respectively. In terms of the number of consistent clusters, ConvexAlign is slightly
better than the second best method BEAMS regardless of c, but much better than the others. In terms of
specificity, ConvexAlign has a much larger advantage over the other methods when c = 4, 5, 6, 7. When
c = 8, ConveAlign is slightly better than BEAMS, but much better than the others. These results indicate
that ConvexAlign aligns proteins in a functionally consistent way, without generating too many spurious
clusters in which the proteins appear to be unrelated. Fig. 2 shows that ConvexAlign outperforms all the
other methods in terms of both M N E and COI. Due to space limit, the topological evaluation is presented
in Fig. A.2.

Figure 1: Specificity (top) and the number of consistent clusters (bottom) generated by the competing methods for different c on synthetic data. Only the best performance for IsoRankN and BEAMS is shown.

3.5

Finding conserved subnetworks

One of the applications of network alignment is to reveal subnetworks conserved across the species. These
subnetworks may be helpful for extracting biological information that cannot be inferred from sequence similarity alone. Fig. 3 shows one conserved complex detected by ConvexAlign among the 5 species: human,
9

Figure 2: COI and MNE of the clusters generated by the competing methods on synthetic data. Only the
best performance for IsoRankN and BEAMS is shown.
yeast, fly, mouse and worm, but not appearing in the alignments generated by other methods. This complex
is enriched for proteasome (with p-value< 10−7 in all species), which is essential for the degradation of
most proteins including misfolded or damaged proteins. The aligned nodes are shown in Table A.1. In
Fig. 3, the interactions in IntAct are displayed in solid lines. For fly, mouse and worm, some edges (shown
by dotted lines) are missing in IntAct but present in the STRING database [29] with experimental evidence
at the highest confidence. Note that our input networks consist of interactions only from IntAct but not
STRING. This suggests that ConvexAlign is able to predict missing interactions. We use PANTHER [17]
to check if the aligned nodes are orthologous proteins. Most of the aligned proteins are shown to be least
divergent orthologs. As shown in Table A.1, there are some missing proteins from different species in some
of the clusters. This is because either there are no orthologs in those species or there is no alignment for
them. For example, cluster 5 has no proteins from worm and fly. PANTHER could not find any orthologous
proteins in those species either. Cluster 6 misses orthologous proteins in fly and yeast, which are aligned
by ConvexAlign to proteins not in this proteasome complex. In addition, this proteasome complex has different number of nodes in different species, which implies that ConvexAlign is able to deal with inserted
and deleted nodes. Fig. 4 shows another conserved subnetwork detected by ConvexAlign that is related to
DNA replication (with p-value< 6−10 in all species). Again, this subnetwork cannot be detected by the
other methods. PANTHER suggests that the aligned proteins are orthologous and functionally related (see
Table A.2).

yeast

human

mouse

fly

worm

Figure 3: The ConvexAlign-detected proteasome complex in each input PPI network.
MCM4

dpa

MCM2
MCM4

Mcm4

MCM6

Mcm6

MCM2
Mcm6

MCM6

mcm-2

mcm-4

Mcm2

MCM7

MCM7

Mcm2

Mcm7

Mcm7

mcm-6

mcm-7

TOP2A

RFC5

Rfc2

RFC4

Top2
LIG3

human

CDC9

rfc-2

RfC4

TOP2

lig3

yeast

fly

Top2a

Lig3

mouse

top-2

lig-4

worm

Figure 4: The ConvexAlign-detected DNA replication complex in each input PPI network.
10

4

Discussion

This paper has presented a new method ConvexAlign for global alignment of multiple PPI networks. ConvexAlign uses a network alignment scoring function that integrates sequence and topological similarity
between the matched proteins as well as interaction consistency. Then ConvexAlign uses a novel convex
formulation to simultaneously align all the proteins in multiple input networks, resulting in better alignment
quality. Such a formulation allows us to use an ADDM method to find its optimal solution.
We have tested ConvexAlign on both real PPI networks and the synthetic data, evaluated the output
alignments by different performance metrics and compared it to several popular methods. Experimental
results on the real data show that on average ConvexAlign generates more functionally consistent clusters
consisting of proteins from most of the input species. That is, ConvexAlign can explain a larger amount of
data in a more functionally meaningful way. ConvexAlign can also find a few conserved and biologically
important complexes which cannot be detected by the other alignment methods.
In the future we may extend ConvexAlign to produce many-to-many global alignments, which will
require some revision of our formulation. It will also be interesting to study how to revise our convex
formulation for local alignments of multiple PPI networks. Of course we may also apply ConvexAlign to
the alignment of other biological systems such as metabolic networks and protein structures.
Currently the time complexity of our algorithm is O(M 3 K), where K is the number iterations in
ADMM and M of the total number of proteins. Using a single computer, it takes dozens of hours to align
the real PPI networks of the five species and only 1.5 hours to align 8 synthetic networks. The M 3 factor
comes from the eigen-decomposition of a M × M matrix, incurred by the consistency constraint X  0.
We may explore a few strategies to speed up this step. For example, we may place the positive semidefinite
constraint on a submatrix of X, i.e., enforcing the consistency among only a subset of important nodes (i.e.,
hub nodes and/or nodes adjacent to hubs). since X is sparse and block-structured, we may also apply some
block-based or parallel algorithms to speed up eigen-decomposition.

References
[1] Ferhat Alkan and Cesim Erten. Beams: backbone extraction and merge strategy for the global manyto-many alignment of multiple ppi networks. Bioinformatics, 30(4):531–539, 2014. 1, 3, 3.2, 3.3
[2] Giovanni Ciriello, Marco Mina, Pietro H Guzzi, Mario Cannataro, and Concettina Guerra. Alignnemo:
a local network alignment method to integrate homology and topology. PloS one, 7(6):e38107–e38107,
2012. 1
[3] Jason Flannick, Antal Novak, Chuong B Do, Balaji S Srinivasan, and Serafim Batzoglou. Automatic
parameter learning for multiple network alignment. In Research in Computational Molecular Biology,
pages 214–231. Springer, 2008. 1
[4] Vladimir Gligorijevic, Noël Malod-Dognin, and Nata?a Prz̆ulj. Fuse: multiple network alignment via
data fusion. arXiv:1410.7585 [q-bio.MN], 2014. 1
[5] Somaye Hashemifar and Jinbo Xu. Hubalign: an accurate and efficient method for global alignment
of protein–protein interaction networks. Bioinformatics, 30(17):i438–i444, 2014. 1, 2.1
[6] Jialu Hu, Birte Kehr, and Knut Reinert. Netcoffee: a fast and accurate global alignment approach to
identify functionally conserved proteins in multiple networks. Bioinformatics, page btt715, 2013. 1,
3, 3.2

11

[7] Qi-Xing Huang and Leonidas Guibas. Consistent shape maps via semidefinite programming. In Proceedings of the Eleventh Eurographics/ACMSIGGRAPH Symposium on Geometry Processing, SGP
’13, pages 177–186, Aire-la-Ville, Switzerland, Switzerland, 2013. Eurographics Association. 2.2
[8] Qixing Huang, Vladlen Koltun, and Leonidas Guibas. Joint shape segmentation with linear programming. ACM Trans. Graph., 30(6):125:1–125:12, December 2011. 2.2
[9] Maxim Kalaev, Mike Smoot, Trey Ideker, and Roded Sharan. Networkblast: comparative analysis of
protein networks. Bioinformatics, 24(4):594–596, 2008. 1
[10] Samuel Kerrien, Bruno Aranda, Lionel Breuza, Alan Bridge, Fiona Broackes-Carter, Carol Chen,
Margaret Duesbury, Marine Dumousseau, Marc Feuermann, Ursula Hinz, et al. The intact molecular
interaction database in 2012. Nucleic acids research, page gkr1088, 2011. 3.1
[11] Wan Kyu Kim and Edward M Marcotte. Age-dependent evolution of the yeast protein interaction network suggests a limited role of gene duplication and divergence. PLoS Comput Biol, 4(11):e1000232,
2008. 3.1
[12] Mehmet Koyutürk, Yohan Kim, Umut Topkara, Shankar Subramaniam, Wojciech Szpankowski, and
Ananth Grama. Pairwise alignment of protein interaction networks. Journal of Computational Biology,
13(2):182–199, 2006. 1
[13] Oleksii Kuchaiev and Nataša Pržulj. Integrative network alignment reveals large regions of global
network similarity in yeast and human. Bioinformatics, 27(10):1390–1396, 2011. 1
[14] M. Pawan Kumar, Vladimir Kolmogorov, and Philip H. S. Torr. An analysis of convex relaxations
for map estimation of discrete mrfs. JOURNAL OF MACHINE LEARNING RESEARCH, 10:71–106,
2008. 2.2
[15] Chung-Shou Liao, Kanghao Lu, Michael Baym, Rohit Singh, and Bonnie Berger. Isorankn: spectral
methods for global alignment of multiple protein networks. Bioinformatics, 25(12):i253–i258, 2009.
1, 3
[16] Noël Malod-Dognin and Nataša Pržulj. L-graal: Lagrangian graphlet-based network aligner. Bioinformatics, page btv130, 2015. 2.1
[17] Huaiyu Mi, Qing Dong, Anushya Muruganujan, Pascale Gaudet, Suzanna Lewis, and Paul D Thomas.
Panther version 7: improved phylogenetic trees, orthologs and collaboration with the gene ontology
consortium. Nucleic acids research, 38(suppl 1):D204–D210, 2010. 3.5
[18] Behnam Neyshabur, Ahmadreza Khadem, Somaye Hashemifar, and Seyed Shahriar Arab. Netal: a
new graph-based method for global alignment of protein–protein interaction networks. Bioinformatics,
29(13):1654–1662, 2013. 1
[19] Rob Patro and Carl Kingsford. Global network alignment using multiscale spectral signatures. Bioinformatics, 28(23):3105–3114, 2012. 1
[20] Catia Pesquita, Daniel Faria, Andre O Falcao, Phillip Lord, and Francisco M Couto. Semantic similarity in biomedical ontologies. PLoS Comput Biol, 5(7):e1000443, 2009. 3.2
[21] Sayed Mohammad Ebrahim Sahraeian and Byung-Jun Yoon. A network synthesis model for generating
protein interaction network families. PloS one, 7(8):e41474, 2012. 3.1

12

[22] Sayed Mohammad Ebrahim Sahraeian and Byung-Jun Yoon. Smetana: accurate and scalable algorithm
for probabilistic alignment of large-scale biological networks. PLoS One, 8(7):e67995, 2013. 1, 3
[23] Vikram Saraph and Tijana Milenković. Magna: Maximizing accuracy in global network alignment.
Bioinformatics, 30(20):2931–2940, 2014. 1
[24] Andreas Schlicker, Francisco S Domingues, Jörg Rahnenführer, and Thomas Lengauer. A new measure
for functional similarity of gene products based on gene ontology. BMC bioinformatics, 7(1):302,
2006. 3.2
[25] Roded Sharan, Silpa Suthram, Ryan M Kelley, Tanja Kuhn, Scott McCuine, Peter Uetz, Taylor Sittler, Richard M Karp, and Trey Ideker. Conserved patterns of protein interaction in multiple species.
Proceedings of the National Academy of Sciences of the United States of America, 102(6):1974–1979,
2005. 1
[26] Rohit Singh, Jinbo Xu, and Bonnie Berger. Pairwise global alignment of protein interaction networks
by matching neighborhood topology. In Research in computational molecular biology, pages 16–31.
Springer, 2007. 1
[27] Rohit Singh, Jinbo Xu, and Bonnie Berger. Global alignment of multiple protein interaction networks
with application to functional orthology detection. Proceedings of the National Academy of Sciences,
105(35):12763–12768, 2008. 1
[28] Ricard V Solé, Romualdo Pastor-Satorras, Eric Smith, and Thomas B Kepler. A model of large-scale
proteome evolution. Advances in Complex Systems, 5(01):43–54, 2002. 3.1
[29] Damian Szklarczyk, Andrea Franceschini, Michael Kuhn, Milan Simonovic, Alexander Roth, Pablo
Minguez, Tobias Doerks, Manuel Stark, Jean Muller, Peer Bork, et al. The string database in 2011:
functional interaction networks of proteins, globally integrated and scored. Nucleic acids research,
39(suppl 1):D561–D568, 2011. 3.5
[30] Andrei Todor, Alin Dobra, and Tamer Kahveci. Probabilistic biological network alignment. IEEE/ACM
Transactions on Computational Biology and Bioinformatics (TCBB), 10(1):109–121, 2013. 1
[31] Alexei Vázquez, Alessandro Flammini, Amos Maritan, and Alessandro Vespignani. Modeling of
protein interaction networks. Complexus, 1(1):38–44, 2003. 3.1
[32] V Vijayan, V Saraph, and T Milenković. Magna++: Maximizing accuracy in global network alignment
via both node and edge conservation. Bioinformatics, page btv161, 2015. 1

13

APPENDIX

A

Results

Figure A.1: AFS distribution resulting from different methods on real data.

0.95	
  
0.9	
  
0.85	
  

CI	
  

0.8	
  
0.75	
  
0.7	
  

DMC	
  

0.65	
  

DMR	
  

0.6	
  

CG	
  

0.55	
  

Ra

Iso

Iso

Ra

nk
N(
0.
3)
	
  
nk
N
Iso
(0
.5
Ra
)	
  
nk
N(
0.
7
SM
)	
  
ET
AN
NE
A
TC 	
  
OF
BE
EE
AM
	
  
S(
0.
3)
BE
	
  
M
A
BE (0.5
AM )	
  
S
Co (0.7
nv
)
ex 	
  
Al
ign
	
  

0.5	
  

Figure A.2: CI analysis on synthetic data.

B Schlicker's similarity
Given two GO terms g1 , g2 and their least common ancestor gc , the Resnik ontological similarity is defined
as simR esnik(g1 , g2 ) = IC(gc ), where IC(g) is the information content of the term g in the given annotation dataset. Considering two gene products p and q annotated with the sets GOp and GOq of GO terms,
respectively with sizes N and M , a similarity matrix S is calculated such as:
sij = simResnik (GOip , GOjq ), ∀i ∈ {1, · · · , N }, ∀j ∈ {1, · · · , M }
This matrix contains all pairwise similarity values between all GO terms associated to p and q. The average
over the row maxima and the column maxima, respectively, gives similarity values for the comparison of p

14

Table A.1: Clusters of aligned nodes that are enriched for proteasome.

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

mouse
Psmb4
Psmd7
Psma1
Psmd14
Psmd1
Rab1A
Psmc6
Psmc3
Psmd2
Psma7
Psmc2
Psmc1
Psmc5
Psmd13
Psmd11
Psmd3
Psmd4

worm
rpn-8
rpn-11
rpt-4
rpt-5
rpn-1
pas-4
rpt-2
rpn-9
-

yeast
fly
human
PUP1
PSMB7
RPN8
Rpn8
PSMD7
PRE5 Prosalpha6T PSMA6
RPN11
Rpn11
PSMD14
RPN2
PSMD1
RAB1A
RPT4
Rpt4
PSMC6
RPT5
Rpt4R
PSMC3
RPN1
Rpn1
PSMD2
PRE6 Prosalpha6 PSMA7
RPT1
Rpt1
PSMC2
RPT2
Rpt2
PSMC1
RPT6
Rpt6
PSMC5
RPN9
Rpn9
PSMD13
RPN6
PSMD11
RPN3
PSMD3
RPN10
PSMD4

Table A.2: Clusters of aligned nodes that are enriched for DNA replication.

cluster 1
cluster 2
cluster 3
cluster 4
cluster 5
cluster 6
cluster 7

mouse
Top2a
Mcm2
Mcm4
Mcm6
Mcm7
Lig3
Rfc2

worm
yeast
fly
human
top-2 CDC25C Top2 TOP2A
mcm-2 MCM2 Mcm2 MCM2
mcm-4 MCM4
dpa MCM4
mcm-6 MCM6 Mcm6 MCM6
mcm-7 MCM7 Mcm7 MCM7
lig-4
CDC9
lig3
LIG3
rfc-2
RFC4
RfC4 RFC5

to q and the comparison of q to p:
scorerow =

N
M
1 X
1 X
max(sij )(1≤j≤M ) , scorecolumn =
max(sij )(1≤i≤N ) .
N
M
(i=1)

(j=1)

The Schlicker's similarity is then calculated as simf unc (p, q) = max{scorerow , scorecolumn }.

C Parameter selection
For all the experiments in this paper, we set the parameters as λ1 = 0.3, and λ2 = 0.02. These parameters
are chosen via 10-fold cross-validation in optimizing the GO-term scores of the alignment between the
mouse and worm networks. The weight factor for aligned interactions is small because: 1) there are many
more aligned interactions than aligned nodes, so a small λ2 may place the node and interaction scores at
the similar scale; and 2) the topological score used in our scoring function already encodes some interaction
information and thus, may overlap with the interaction score. Of course we may increase λ2 to favor other
performance metrics such as the number of aligned interactions and the number of annotated clusters.

15

D Proof of Proposition 1
Proof. Let Vi be a vertex set corresponding to all rows in the block Xii . Then Xii has size |Vi | × Vi . Each
block Xij has size |Vi | × |Vj |, describing the relationship between Vi and Vj . Let V = ∪1≤i≤N Vi , then the
binary block matrix X has size |V | × |V |. Ignoring the N blocks Xii (1 ≤ i ≤ N ), the binary block matrix
X can be treated as an adjacency matrix of a simple graph M. That is, starting from X, we may construct
a simple graph M for V such that one non-zero entry in X corresponds to one edge in M.
According to the constraints Xij 1 ≤ 1, XijT 1 ≤ 1, it is easy to see that each connected component of
M contains at most one vertex from each Vi . Now we want to prove that each connected component in M
is a clique. This is equivalent to proving the consistency property, i.e., given three vertices vi ∈ Vi , vj ∈
Vj , vk ∈ Vk , if Xij (vi , vj ) = Xik (vi , vk ) = 1, then Xjk (vj , vk ) = 1. This can be induced from that X is
positive semidefinite. Consider the 3 × 3 principal submatrix of X induced by vi , vj , vk . Since the principle
submatrix of a positive semidefinite matrix is also semidefinite positive, we have


1
1
1
 1
1
Xjk (vj , vk )   0.
1 Xjk (vj , vk )
1
It follows that


1
1
1
1
Xjk (vj , vk ) ) = −(1 − Xjk (vj , vk ))2 .
det( 1
1 Xjk (vj , vk )
1


This implies Xjk (vj , vk ) = 1.
Let A = {A1 , A2 , ..., AK } denote all the connected components of M. Then A is a feasible one-to-one
alignment between the N vertex sets V1 ,V2 ,...,VN . For each vertex Vi , we construct a binary matrix Yi of
size |Vi | × K. For any k (1 ≤ k ≤ K) and vi ∈ Vi , if vi appears in Ak , then Yi (vi , k) = 1; otherwise
Yi (vi , k) = 0. Finally we construct a |V | × K binary matrix Y by stacking Y1 , Y2 ,...,YN along the vertical
direction. It is easy to show that X = Y Y T .

E Optimization via ADMM (Alternating Direction Method of Multipliers)
In this section, we describe in detail how to solve the optimization problem in (13) using ADMM. The
basic idea is to augment its Lagrangian dual (see ??) and iteratively optimize a subset of variables while
keeping the others fixed. This allows us to exploit structure patterns in the constraint set for effective
optimization. To maximize the power of ADMM, we introduce a latent variable sij to break each constraint
λ2
Bij y ij ≤ Fij (Xij ) into two sets of constraints Bij y ij ≤ sij and sij = Fij (Xij ). Let λ = 1−λ
and dij be
2

16

the coefficient of y ij . The relaxed convex optimization problem can be written as follows.

X 
maximize
hCij , Xij i + λhdij , yij i
1≤i<j≤N

: z 0ij ≥ 0

subject to y ij ≥ 0, 1 ≤ i ≤ j ≤ N
Bij y ij ≤ sij , 1 ≤ i ≤ j ≤ N

: z ij,1 ≥ 0

sij = Fij (Xij ), 1 ≤ i < j ≤ N

: z ij,2

Xij 1 ≤ 1, 1 ≤ i ≤ j ≤ N

: z ij,3 ≥ 0

XijT 1 ≤ 1, 1 ≤ i ≤ j ≤ N

: z ij,4 ≥ 0

Xij ≥ 0, 1 ≤ i ≤ j ≤ N

: Zij,5 ≥ 0

Xii = I|Vi | , 1 ≤ i ≤ N

: Zii,6

X0

:S0

(14)

Note that the right column shows the dual variables of the corresponding constraints. Using the dual variables, the Lagrangian of the above problem is as follows:
X 
− hCij , Xij i − λhdij , y ij i − hz ij,0 , y ij i + hz ij,1 , Bij y ij − sij i + hz ij,2 , sij − Fij (Xij )i
L=
1≤i<j≤N


+ hz ij,3 , X1 + 1i + hz ij,4 , X T 1 − 1i − hZij,5 , Xij i − 2hSij , Xij i

X 
+
hZii,6 , Xii − I|Vi | i − hSii , Xii i
1≤i≤N

X

=



− h1, z ij,3 i − h1, z ij,4 i − hCij + FijT (z ij,2 ) − z ij,3 1T − 1z Tij,4 + Zij,5 , Xij i

1≤i<j≤N
T
− hλdij + z ij,0 − Bij
z ij,1 , y ij i − hz ij,1 − z ij,2 , sij i

X 
+
− hI|Vi | , Zii,6 i − hSii − Zii,6 , Xii i .



1≤i≤N

(15)
The augmented Lagrangian dual of the above problem is as follows.
X 
L0 =
h1, z ij,3 i − h1, z ij,4 i + hCij + FijT (z ij,2 ) − z ij,3 1T − 1z Tij,4 + Zij,5 + 2Sij , Xij i
1≤i<j≤N
T
+ hλdij + z ij,0 − Bij
z ij,1 , y ij i + hz ij,1 − z ij,2 , sij i

1
T
+
kCij + FijT (z ij,2 ) − z ij,3 1T − 1z Tij,4 + Zij,5 + 2Sij k2F + kλdij + z ij,0 − Bij
z ij,1 k2 + kz ij,1 − z ij,2 k2
2µ

X 
1
+
hI|Vi | , Zii,6 i + hSii − Zii,6 , Xii i +
kSii − Zii,6 k2F .
(16)
2µ
1≤i≤N

L0

shall be maximized with respect to the primal variables but minimized with respect to the dual variables. We initialize all the primal and dual variables to zero. At iteration k + 1, we update the dual and
primal variables as follows.
(k+1)
Step 1: Optimizing z ij,0 . When variable z ij,0 is active with other variables fixed, the optimization
problem is equivalent to computing
(k)

(k)

T
min kλdij + z ij,0 − Bij
z ij,1 + µy ij k2 ,

z ij,0 ≥0

17

1≤i<j≤N

In this case, we have
(k+1)

z ij,0

(k)

(k)

T
= max(0, Bij
z ij,1 − µy ij − λdij ),

1≤i<j≤N

(k+1)

Step 2: Optimizing z ij,1 . When z ij,1 is active while other variables are fixed, each of them can be
optimized independently as


(k+1)
(k)
(k)
T
min kBij
z ij,1 − λdij + z ij,0 + µy ij k2 + kz ij,1 − z ij,2 k2 .
z ij,1 ≥0

Since each y appears in 4 constraints of (10), we permute the rows of Bij to form 4 submatrices
Bij,1 , Bij,2 , Bij,3 , Bij,4 , so that each column of one submatrix contains exactly one non-zero entry. Accordingly, we can also reorder and rewrite z ij,1 as (z ij,1,1 ; z ij,1,2 ; z ij,1,3 ; z ij,1,4 ). Then we alternate the
optimization of z ij,1,l (l = 1, 2, 3, 4) as follows.


X
(k+1)
(k+1)
(k)
(k)
T
T
z ij,1,l − λdij + z ij,0 + µy ij −
Bij,p
z ?ij,1,p k2 + kz ij,1,l − z ij,2,l k2 .
z ij,1,l = arg min kBij,l
z ij,1,l ≥0

p6=l

(17)
where z ?ij,1,p is the latest value of z ij,1,p . Due to the special structure of Bij,l (l = 1, 2, 3, 4), the elements of
z ij,1,l can be optimized independently, leading to explicit expressions of optimal values:
(k+1)

(k+1)

T
+ I)−1 (Bij,l (λdij + z ij,0
z ij,1,l = max 0, (Bij,l Bij,l

(k)

+ µy ij −

X

(k) 
T
Bij,p
z ?ij,1,p ) + z ij,2,l ) .

p6=l
T + I is a diagonal matrix and z ?
where Bij,l Bij,l
ij,1,p is the latest value of z ij,1,p .
(k+1)

Step 3: Optimizing z ij,2 , z ij,3 , z ij,4 , Zij,5 and Zii,6 . The optimization of each z ij,2 is decoupled and its
optimal value is


(k)
(k) T
(k)
(k)
(k)
(k+1)
(k+1)
z ij,2 = arg min kF T (z ij,2 ) − z ij,3 1T + 1z ij,4 − C − Zij,5 − 2Sij − µXij k2F + kz ij,2 − z ij,1 k2
z ij,2


(k)
(k) T
(k)
(k)
(k) 
(k+1)
(18)
= (FF T + I)−1 F z ij,3 1T + 1z ij,4 − C − Zij,5 − 2Sij − µXij + z ij,1 .
where FF T + I is a diagonal matrix. Through a similar derivation, we can compute the optimal value of
other variables at iteration k + 1 as


(k+1)
(k+1)
(k)
(k)
z ij,3 = arg min kz ij,3 1T − Cij + F T (z ij,2 ) + Zij,5 + µXij − 1z ij,4 (k) k2F + 2µh1, z ij,3 i
z ij,3 ≥0



(k+1)
(k)
(k)
= max(0, Cij + F T (z ij,2 ) + Zij,5 + µXij − 1z ij,4 (k) 1 − µ /|Vj |)
(19)

(k+1)

z ij,4



(k+1)
(k)
(k)
= arg min k1z Tij,4 − Cij + F T (z ij,2 ) + Zij,5 + µXij − z ij,3 (k+1) 1T k2F + 2µh1, z ij,4 i
z ij,4 ≥0


T
(k+1)
(k)
(k)
= max(0, Cij + F T (z ij,2 ) + Zij,5 + µXij − 1z ij,4 (k) 1 − µ /|Vi |)
(20)

18

(k+1)

Zij,5



(k+1)
(k)
(k)
(k+1)
(k+1) T
− Cij − F2T (z ij,2 ) − µXij − 2Sij k2F
= arg min kZij,5 − z ij,3 1T + 1z ij,4
Zij,5 ≥0



(k+1)
(k)
(k)
(k+1)
(k+1) T
− Cij − F2T (z ij,2 ) − µXij − 2Sij ,
= max 0, z ij,3 1T + 1z ij,4
(k+1)

Zii,6

= arg min kSii + µXii − Zii,6 − µI|Vi | k2F
Zii,6 ≥0


(k)
(k)
= max 0, Sii + µXii − µI|Vi | .

(21)

Step 4: Optimizing S. Finally, we optimize S. In this case, the optimization problem is reduced to
S (k+1) = arg min
S≥0

X

(k+1)

(k+1) T

(k+1)

+ Zij,5

(k)

+ 2Sij + µXij k2F

1≤i<j≤N

X

+

(k+1)

kCij + FijT (z ij,2 ) − z ij,3 1T − 1z ij,4
(k+1)

kSii − Zii,6

(k)

+ µXii k

(22)

1≤i≤N

= arg min kS − T (k) k2F
S≥0

= U max(Σ, 0)U T ,

(23)

where U ΣU T is eigen-decomposition of T (k) defined below and max(Σ, 0) takes the positive eigenvalues.
 

 z (k+1) 1T + 1z (k+1) T − C − F T (z (k+1) − µX (k) /2
i 6= j
ij
(k+1)
ij
ij,3
ij,4
ij,2
ij
Tij
=
(24)
(k+1)
(k)

Z
− µX
otherwise
ii,6

ii

Step 5: Optimizing primal variables. Finally the primal variables are updated as follows:
(k+1)

y ij

(k+1)

Xij

s(k+1)


1
(k+1)
(k+1)
λ1 + z ij,0 − B T z ij,1 ,
µ

1
(k+1)
(k)
(k+1)
(k+1)
(k+1) T
(k+1)
+ Zij,5 + 2Sij
= Xij +
Cij + FijT (z ij,2 ) − z ij,3 1T − 1z ij,4
µ
1 (k+1)
(k+1)
= s(k) + (z 1
− z2
).
µ
(k)

= y ij +

19

(25)

