Optimizing over trained GNNs via symmetry breaking
Authors: Shiqiang Zhang
Juan S Campos
Christian Wolfgang Feldmann
David Walz
Frederik Sandfort
Miriam Mathea
Calvin Tsay
Ruth Misener
Venue: NeurIPS 2023 poster (2023)
Keywords: Mixed-integer optimization
Graph neural network
Symmetry-breaking
Molecular design

Abstract: Optimization over trained machine learning models has applications including: verification, minimizing neural acquisition functions, and integrating a trained surrogate into a larger decision-making problem. This paper formulates and solves optimization problems constrained by trained graph neural networks (GNNs). To circumvent the symmetry issue caused by graph isomorphism, we propose two types of symmetry-breaking constraints: one indexing a node 0 and one indexing the remaining nodes by lexicographically ordering their neighbor sets. To guarantee that adding these constraints will not remove all symmetric solutions, we construct a graph indexing algorithm and prove that the resulting graph indexing satisfies the proposed symmetry-breaking constraints. For the classical GNN architectures considered in this paper, optimizing over a GNN with a fixed graph is equivalent to optimizing over a dense neural network. Thus, we study the case where the input graph is not fixed, implying that each edge is a decision variable, and develop two mixed-integer optimization formulations. To test our symmetry-breaking strategies and optimization formulations, we consider an application in molecular design.

Optimizing over trained GNNs via symmetry breaking
Shigiang Zhang" Imperial College London London, UK
Juan S. Campos Imperial College London London, UK
Christian Feldmann
BASF SE Ludwigshafen, Germany
David Walz BASF SE
Frederik Sandfort
Miriam Mathea
BASF SE
BASF SE Ludwigshafen, Germany
Ludwigshafen, Germany
Ludwigshafen, Germany
Calvin Tsay Imperial College London London, UK
Ruth Misener
Imperial College London London, UK
Abstract
Optimization over trained machine learning models has applications including: verification, minimizing neural acquisition functions, and integrating a trained surrogate into a larger decision-maie making problem. This paper formulates and solves optimization problems constrained by trained graph neural networks (GNNs). To circumvent the symmetry issue caused by graph isomorphism, we propose two types of symmetry-breakinn the remaining nodes by lexicographically ordering their neighbor sets. To guarantec that adding these constraints will not remove all symmetric solutions, we construct a graph indexing algorithm and prove that the resulting graph indexing satisfies the proposed yymmetry-breaking constraints. For the classical GNN architectures considered in this paper, optimizing over a GNN with a fixed graph is equivalent to optimizing a dense neural network. Thus, study the case where the input graph is not fixed, implying that each edge a decision variable, and develop two mixed-intege optimization formulations. To test symmetry-breakee strategies and optimization formulations, we consider an application in molecular design.
constraints: one indexing node 0 and one indexing
Introduction
Graph neural networks (GNNs) 13 are designed operate on graph-structured data. By passing messages between nodes via edges (or vice versa), GNNs can efficiently capture and aggregate local information within graphs. GNN architectures including spectral approaches [49 and spa- tial approaches [8), are proposed based on various motivations. Due to their ability to learn non-Euclidean data structure, GNNs have recently been applied to many graph-specific tasks, demon- strating incredible performance. In drug discovery, for example, GNNs can predict molecular properties given the graph representation of molecules [16
19 Other applications using GNNs
include graph clustering 24 25.. text classification 266 27 and social recommendation [28]] 29
Although GNNs are powerful tools for these "forward" prediction tasks, few works discuss the "backward" (or inverse) problem defined on trained GNNs. Specifically, many applications motivate the inverse problem of using machine learning (ML) to predict the inputs corresponding to some
output specification(s). For example, computer-aided molecular design (CAMD) 30 32 aims to
"Corresponding author: s.zhang21 @imperial.ac.uk
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Prediction (Forward): What the properties for given molecule?
molecules
GNN
properties
quantum mechanics
physico-chemical bioactivity
(cco) (eco)-toxicityy
Optimization (Backward): What is the optimal molecule with desired properties?
Figure 1: Illustration of forward and backward problems defined on trained GNNs.
design the optimal molecule(s) for a given target based on trained ML model (3 Several GNN-based approaches have been proposed for CAMD 47, however, these works still use GNNs as "forward" functions to make predictions over graph-domain inputs. Therefore, both the graph structures of inputs and the inner structures of the trained GNNs are ignored. Figure I conceptually depicts the difference between forward and backward/inverse problems.
Mathematical optimization over trained ML models is an area of increasing interest. In this paradigm, a trained ML model translated into an optimization formulation, thereby enabling decision decision-making making problems over model predictions. For continuous, differentiable models, these problems can be solved via gradient-basee methods -| Mixed-integer programming (MIP) has also been proposed, mainly to support non-smooth models such as ReLU neural networks 45-50), their ensembles [51], and tree ensembles (52) Optimization over trained ReLU neural networks has been especially prominent [55). finding applications such as verification 65), reinforcement learning (65-67],
compression (68), and black-box optimization [69]. There is no reason to exclude GNNs from this rich field. Wu et al. (70] managed to verify GNN-based job scheduler, where MIP is used obtain tighter bounds in a formulations for GCN 6 and GraphSAGE [13 and tested their performance
abstraction refinement framework. Mc Donald [71] built MIP
CAMD problem, which is the first work to directly apply MIP on GNNs to the best of our knowledge.
Optimization over trained GNNs involves two significant challenges. One is the symmetry issue arising from the permutation invariance of GNNs (i.e., isomorphism graphs share the same GNN output(s)). For forward problems, this invariance is preferred since it admits any graph representation for inputs. From optimization perspective, however, symmetric solutions can significantly enlarge the feasible domain and impede a tranch-and-bous
(B&B) search tree. One way to break symmetry in integer programs is by adding constraints to remove symmetric solutions (2]. The involvement of graphs not only adds complexity to the problem in terms of symmetry, but also in terms of implementation. Due to the complexity and variety of GNN architectures, a general framework is needed. This framework should be compatible with symmetry breaking techniques.
This paper first defines optimization problems on trained GNNs. To handle the innate symmetry, we propose two sets of constraints to remove most isomorphism graphs (i.e., different indexing for any abstract graph). To guarantee that these constraints do not exclude all possible symmetrics (i.c., any abstract graph should have at least one feasible indexing), undirected graph and prove that the resulting indexing is feasible under the proposed constraints. To encode GNNs, we build general framework for MIP GNNs based on the open-source Optimization and Machine Learning Toolkit (OMLT) 79 Many classic GNN architectures are supported after proper transformations. Two formulations are provided and compared on a CAMD problem. Numerical results show the outstanding improvement of symmetry- breaking.
design an algorithm to index any
Paper structure: Section introduces the problem definition, proposes the -breaknn constraints, and gives theoretical guarantees. Section 3 discusses the application, implementation details, and numerical results. Section 4ection[[coneee concludes and discusses future work.

---

2 Methodology
2.1 Definition of optimization trained GNNs
We consider optimization over a trained GNN on a given dataset D = ((X', A*),,A"), M Each sample in the dataset consists of input graph (with features X and adjacency matrix A') and an output property y². The target of training is to approximate properties, that is:
GNN(X' A') = y', VI iisss M
When optimizing over a trained GNN, the goal is to find the input with optimal property:
(X", A²) = arg min GNN(X, A) (X,A)
(OPT)
s.t. fy (X, A) S S.jEET 9s [X, A) = =
where fj. 9s are problem-spp sepeeifi- constraints and J, are index sets, Note that optimality (OPT) is defined on the trained GNN instead of true properties. To simplify our presentation, we will focus on undirected graphs with node features, noting that directed graphs can be naturally transformed into undirected graphs. Likewise, using edge features does not influence our analyses. We discuss practical impacts of such modifications later.
2.2 Symmetry handling
For any input (X, A). assume that there are N nodes, giving a N x N adjaceney matrix A and a feature matrix X (Xo, XI...., ...,,,N 1N= ) Nxf, where X; E Rº contains the features for i-th
node, For any permutation over set [N] == {0,
1], let:
A = An(w),(()) X = X,()). Vu, v E [N]
Then the permuted input (X', A') has isomorphic graph structure, and therefore the same output, due the permutation invariance of GNNs (i.e., GNN A) GNN AA)). In other words, the symmetries result from different indexing for the same graph. In general, there exist N! ways to index N nodes, each of which corresponds to solution of (OPT). Mc Donald [71] proposed set of constraints to force molecules connected, which can also be used to break symmetry in our general setting. Mathematically, they can be written
Ve [N]](}}, Eu < c, s.t. Au,e =1
(S1)
Constraints (S1) require each node (except for node 0) to be linked with node with smaller index. Even though constraints (ST) help break symmetry, there can still exist many symmetric solutions. To resolve, at least relieve, this issue, need to construct more commetry--rreking constraints.
Breaking symmetry the feature level. We can first design some constraints on the features to define a starting point for the graph (i.e., assigning index () to a chosen node). Note that multiple nodes may share the same features. Specifically, we define an arbitrary function h RRF. - RR to assign a hierarchy to each node and force that node 0 has the minimal function value:
h(Xo) S h(X,), Yu E [N]][O}
(S2)
The choice of l can be application-speeificc for example, as in Section 3.1 The key is to define such that not too many nodes share the same minimal function value.
Breaking symmetry at the graph level. is unlikely that constraints on features will break much of the symmetry. After adding constraints (S1) and (S2), any neighbor of node could be indexed 1, then any neighbor of node 0 or could be indexed 2, and so on. We want to limit the number of possible indexing. natural idea is account for neighbors of each node: the neighbor set of node with smaller index should also have smaller lexicographical order.
Before further discussion, need to define the lexicographical order. Let S(M,, be the set of all non-decreasing sequences with L integer elements [0, MM (i.c., S(M, L) = {(11, .2.....a........................ at) 0 S a1 I 02 < a1 < M)). For any S(M, L), denote its lexicographical order by LO(a). Then for any a, E S(M, L), we have:
(a) < T.O(b) + B1E (,,2,.... D.,, s.l. [a b. 15i<l
(L.O)
aa b
3

---

For any multiset A with no more than L integer elements in [0, M 11, we first sort all elements in A in a non non-decreasing order, then add elements with value M until the length equals L. In that way, we build an injective mapping from A to a sequence E SMM, T.). Define the lexicographical order of A by (a) For simplicity, we use (A) to denote the lexicographical order of A.
Algorithm 1 Indexing algorithm
Input: G (V. F) with node set V (Ub, 01, ................................ UN-1) (N [V). Denote the neighbor set of node e as (u), Vu E V.
I(vb) +0 0
Assume that vo is indexed with O
--1 19 + {vo}
Index for next node
Initialize set of indexed nodes
while N do
V - V
b Set of unindexed nodes
< - (I(4) (EE(()) Vo V V
D Obtain all indexed neighbors
(v) - (LO((\(())) LO(N(())) Yu E Vi Ye E V D Assign a rank to each unindexed node
[v) - [(v),
VEE (1) + E V
b Assign temporary indexes
NT(v) + (I*(w) u Mu (())), Vu e V
> Define temporary neighbor sets based on I > Neighbors of ve has minimal order multiple nodes share the same minimal order, arbitrarily choose one > Index to node "*
vd + arg min LO(N(())
NEVY
I(v²)= =
vet1 t VF
> Add De to set of indexed nodes
5 ++1
D Next index is + 1
end while
Output: u(v).
Result indexing
The definition of I.O(-) and Section 2.3 reuse A to denote multiset. This will not cause ambiguity since multisets are only introduced the theoretical part, where the adjacency matrix is not involved. For brevity, use capital letters A, B.... to denote multisets and lowercase letters a, b.... to denote the corresponding sequences in L),
Remark:
Since we only need lexicographical orders for all neighbor sets of the nodes of the graph, let M N, N 1. With the definition of LO(.), we can represent the aforementioned idea by: (OVN()) + 1}) S LO(N(v + 1){{u)) Vu e [N - -]}{(}
(S3) where ((v) is the neighbor set of node v. Note that the possible edge between nodes v and v + is ignored in (S3) to include the cases that they are linked and share all neighbors with indexes less than In this case, LOON((v + 1)) is necessarily smaller than (w))(w)) since E N + 1).
Constraints S exclude many ways index the rest of the nodes, for example they reduce the possible ways to index the graph in Figure from 120 to 4 ways. But, we still need to ensure that there exists at least one feasible indexing for any graph after applying constraints (S3). Algorithm constructs an indexing and Section 2.3 provides theoretical guarantee of its feasibility.
2.3 Theoretical guarantee
This section proves that Algorithm always provides feasible indexing. We first give some properties that are used to derive an important lemma.
< Ve Vi
Property 1. For any s 1,2,..., N 1,
-1, 2s, Vu V/2

---

V0
U5
Figure 2: Consider a graph having 6 nodes with different features. Without any constraints, there are 6! 720 ways to index Utilizing (ST) force connectivity results in 636 ways. Using (S2) to choose vo to be indexed 0, then there are 5! = 120 ways. Finally, applying (S3) to order the rest of the nodes, there are only ways, which can be derived from Algorithm For details on using Algorithm to index this graph, see Appendix[A.2
Property 2. For any $1, 82 NI1,
$1 $1 = (v) N Vu V22
Property 3. Given any multisets A, B with more than integer elements in |0, M 1]. have: LO(A) & LO(B) = LO(An [m]) s LO(B n [m]), Vm = 1, ..,
Using these properties, smaller than v, then at each iteration, the temporary index assigned to u not greater than v.
prove Lemma D which shows that if the final index of node is
Lemma 1. For any two nodes and v.
I(u) < I(v) Iu) [v),, vs = =,2=,2,, .. N
Now can prove Theorem
Theorem 1. Given any undirected graph G (V, with one node indexed 0. The indexing yielded from Algorithm 7 satisfies SS3
Proof. Assume that constraints SS3 not satisfied and the minimal index that violates S3 is
Denote nodes with index s and s by u, respectively (i.c., I(u) = s, s,I(v) s + 1).
Let [(u)(() == {u1, u2. ................................ Mm um be all neighbors of u except for where:
I(u:) < I(u++1), Vi 1, mm
Similarly, let (())) {u,, 12 ... )n be all neighbors of except for where: I(0;) < I(1441), Vj = 1,2, .2
Denote the sequences in S(N,N 1) corresponding to sets (I(u,) 11 ii m] and (I(e;) 11s j S n) by (0), aN 12,...,aa-11. = (b,, b,.,,.,,,,,,,,,. By definition of LO(.):
[I(u;), N.
1
i N by = [(us), N,
1 S
n n< < N
Since nodes u and v violate constraints
S3 there exists a position
1 < N --1 satisfying:
0; bi:
v1 sick
by
from where we know that nodes u and v share the first k - 1 neighbors (i.c. = U. V1 < < k). From by ax < N know that node definitely has its k-th neighbor node Ug. Also, note that Uk is not a neighbor of node u. Otherwise, we have u4 and then ag = I(ux) = I(vx) by-
Case 1: If az N, that is, node has neighbors.
In this case, node has neighbors of node u as well as node v.. Therefore, we have:
TO(NF( (O()((e))
5

---

which violates the fact that node u is chosen to be indexed s at s-th iteration of Algorithm 1
Case 2: If ag. N, that is, node has nodes uy its k-tth neighbor.
Since I(ux) ok > by I(mx), we can apply Lemma on node ug and uy (s + 1)-th iteration, to obtain
I (ux) 2 (vx)
Similarly, if we apply Lemma 1 to all the neighbors of node and node at s-th iteration, we have:
I' (w:) S (u++1) Vi 1,2....... m I'())) S (1)+1)) Vj = .
Given that O = -(ux) is the k-th smallest number in a, we conclude that I (ux) ) is equal to the k-th smallest number in Mf()). Likewise, I (va) equals to the k-th smallest number in (v).).. Meanwhile, I"(u;) = (u;) since u, = 2/ A1 S < k. After comparing the lexicographical orders between of (u) and W/ (with the same k smallest elements, (ux) and Ie (vx) as the k-th smallest element, respectively), node u is chosen. Therefore, we have:
(*k)) S (va))
from which we know that:
LO(N (ux)) LO(N((v)))
At (s + 1)-th iteration, node uy has one more indexed neighbor (i.e., node a with index s), while node U/ has no new indexed neighbor. Thus we have:
LO(V+t((u))
(ss)) < LO(N((w)) LO(N"(e)) LO(W++((o))
which yields:
I (wx) < I
(<)
The contradiction between and completes this proof.
Given any undirected graph with node features, after using S2 to choose node 0, Theorem guarantees that there exists at least one indexing satisfying (S3). However, when applying (ST
force a connected graph, we need to show the compatibility between SI and S3 as shown in Lemma 2 [A proofs of properties and lemmas.
Lemma 2. For any undirected, connected graph G = (V, E), if one indexing of G sastifies S3 then it satisfies (II..
Note that generic graph-based search problems, as long as there exists a symmetry issue caused by graph isomorphism. GNNs can be generalized to all permutation invariant functions defined over graphs. Although the next sections apply MIP formulations to demonstrate the symmetry-breakinnn constraints, could alternatively use constraint programming [80] satisfiability [8]] paradigm. For example, constraints in Reluplex [3].
S3 are not limited to optimizing over trained GNNs. In fact, they can be employed
we could have encoded the Elgabou and Frisch
2.4 Connection Difference to the symmetry-breaking literature
Typically, MIP solvers detect symmetry using graph automorphism, for example SCIP 84 uses BLISS 85 86), and both Gurobi (8 and SCIP break symmetry using orbital fixing [8 and orbital pruning [89]. When MIP solver detects symmetry, the only graph available is the graph formed by the variables and constraints in the MIP.
Our symmetries come from alternative indexing of abstract nodes. Each indexing results in an isomorphic graph. In MIP, however, each indexing corresponds an element of much larger symmetry group defined on all variables, including node features (O(NF)), adjacency matrix (O(N2)), model (e.g.. GNN), and problem-spectie
features. For instance, the first row of Table 6 the input graph has N = 4 nodes, but the MIP has 616 + 411 variables. We only need to consider permutation group with N! = 24 elements. However, because the MIP solver does not have access to the input graph structure, it needs to consider all possible automorphic graphs with 616 + 411 nodes. By adding constraints (ST) (S3). there no need consider the symmetry group of all variables to find a much smaller subgroup corresponding to the permutation group defined abstract nodes.

---

The closest setting to ours is distributing m different jobs to identical machines and then minimizing the total cost. Binary variable ; denotes if job assigned machine j. The requirement that each job can only be assigned to one machine (but each machine can be assigned to multiple jobs). Symmetries come from all permutations of machines. This setting appears in noise dosage problems 190, 19], packing and partitioning orbitopes [92] 931, and scheduling problems [94]. However, requiring that the sum of each row in Ai,; equals to 1 simplifies the problem. By forcing decreasing lexicographical orders for all columns, the symmetry issue is handled well. Constraints (S3) can be regarded as a non-trivial generalization
these constraints from a bipartite graph to an arbitrary undirected graph: following Algorithm I will produce the same indexing documented in [90 94
2.5 Mixed-integer formulations for optimizing over trained GNNs
As mentioned before, there are many variants of GNNs 4418 with different theoretical and practical motivations. Although it is possible to build MIP for a specific GNN from the scratch (e.g., (7]).
general definition that supports multiple architectures is preferable.
2.5.1 Definition of GNNs
We define a GNN with L layers as follows:
GNN: Rd R
Rd ... 2 Re
[V| times
[V]
where V is the set of nodes of the input graph.
Let x iRdo be the input features for node v. Then, the 1-th layer (1 ...,.... .......,)....... is defined by:
(
w 11) 10 11-1) (0-1)
+ b00 AVEV
where ((v) is the set of all neighbors of e, a could be identity any activation function.
With linear aggregate functions such
and mean, many classic GNN architectures could be
rewritten in form (x). for example: Spectral Network [4], ChebNet (5), GCN 6), Neural FPs [10
DCNN [1I, PATCHY-SAN [12, GraphSAGE 13] and MPNN [66
2.5.2 Mixed-integer optimization formulations for non-fixed graphs
the graph structure for inputs is given and fixed, then ++) is equivalent to a fully connected layer,
whose MIP formulations already well-established 46 49 But if the graph structure is non non-fixed (i.e., the elements adjacency matrix are also decision variables). two issues arise in : (1)((v) is not well-defined; (2) !o. and 6! may be not fixed and contain the graph's 's information. Assuming that the weights and biases are constant, we can build two MIP formulations to handle the first issue.
The first formulation comes from observing that the existence of edge from node to node determines the contribution link from x!!-1) (!-1) x.!) Adding binary variables cu-sv for all u, V. we can then formulate the GNNs with bilinear constraints:
x a Σ NEV
Cu-ty
(1-1)) + 6(D) VVEV
(bilinear)
Introducing this nonlinearity results in mixed-integer quadratically constrained optimization problem (MIQCP), which can be handled by state-of-the-art solvers such as Gurobi 87
The second formulation generalizes the big-M formulation for GraphSAGE in to all GNN architectures satisfying (*) and the assumption of constant weights and biases. Instead of using binary variables to directly control the existence of contributions between nodes, auxiliary variables 20-1e are introduced to represent the contribution from node u to node in the (-th layer:
(1-1)
x!!!
Σ Ewi
(L! (1-1) + b000
Ve E
(big-M)
7

---

where:
(111
0,
Ex-tv =
(1-1) ca-t1
atu
Assuming that each feature is bounded, the definition of =u1-1' can be reformulated using big-M:
x 1) (1 Cu-w) S=(=a S n/ =x
M
11 eu-ro)
(1)))) ( cu-o 5x((-1) & M (1-11)e--r ex-
where aaa 1 < (
Cu-+v {0, 1). By adding extra continuous variables and constraints, as well as utilizing the bounds for all features, the big-M formulation replaces the bilinear constraints with linear constraints. Section numerically compares these two formulations.
3 Computational experiments
We performed all experiments on a 3.2 GHz Intel Core 17-8700 CPU with 16 GB memory. GNNs are implemented and trained in PyG [95]. MIP formulations for GNNs and CAMD implemented based on OMLT [79), and are optimized by Gurobi 10.0.1 87 with default relative MIP optimality gap (i.e., 10 10-4). +). The code is available at https://githubb.. www//ooo-mmeriiallllllllooooo
MIP CAMD These are all engineering choices and we could have, for example, extended software other than OMLT to translate the trained GNNs into an optimization solver 996 97].
3.1 Mixed-integer optimization formulation for molecular design
MIP formulations are well-established in the CAMD literature 98-110 The basic idea is represent- ing molecule as a graph, creating variables for each atom (or group of atoms), using constraints to preserve graph structure and satisfy chemical requirements. A score function is usually given as the optimization target. Instead of using knowledge from experts statistics to build the score function, GNNs (or other ML models) could be trained from datasets to replace score functions. Here we provide general and flexible framework following the MIP formulation for CAMD in (71) Moreover, by adding meaningful bounds and breaking symmetry, our formulation could generate more reasonable molecules with less redundancy. Due to space limitation, briefly introduce the formulation here (see Appendix B for the full MIP formulation).
To design a molecule with N atoms, we define N x F binary variables Xe,f, v E [N], E [F] to represent F features for N atoms. Hydrogen atoms are not counted in since they can implicitly considered as node features. Features include types of atom, number of neighbors, number of hydrogen atoms associated, and types of adjacent bonds. For graph structure of molecules, we add three sets of binary variables An,v; DB...,, TB.,, to denote the type of bond (i.e., any bond, double bond, and triple bond) between atom u and atom v, where A is the adjacency matrix.
C2 handle structural feasibility following (71). Additionally, we propose new constraints bound the number of each type of atoms (C22). double bonds (C23). triple bonds (C24), and rings (C25). In our experiments, we calculate these bounds based on the each dataset itself so that each molecule in the dataset will satisfy these bounds (see Appendix C for details). By setting proper bounds, we can control the composition of the molecule, and avoid extreme cases such as all atoms being set to oxygen, or molecule with too many rings or double/triple bounds. In short, constraints (C22). (C25) provide space for chemical expertise and practical requirements. Furthermore, our formulation could be easily applied to datasets with different types of atoms by only changing the parameters in Appendix B.I Moreover, group representation of molecules [98] 102 is also compatible with this framework. The advantage is that all constraints can be reused without any modification.
To design reasonable molecules, constraints CC
constraints (C5) are the realization of (SI). Except for (C5), these structural constraints are independent of the graph indexing. Therefore, we can compatibly implement constraints (S2) and (S3) to break symmetry.
Among constraints (CI)
C25) (as shown in Appendix B.3
Corresponding to (2), we add the following constraints over features:
F-- No.5 V
Σ
2- X Ve E [N]((0}
(C26)
SE(F]
SE(F]
8

---

where 2 [F] are coefficients to help build a bijective h between all possible features and
all integers [0, 20 11. These coefficients are also called "universal ordering vector" in 72 78).
On graph level, constraints (S3) can be equivalently rewritten as:
2 Au,
2~--1 Au,tt, VU E IN 1JM(0)
(C27)
a
Similarly, the coefficients 2 u E [N] are used to build bijective mapping between all possible sets of neighbors and all integers in [0, 1].
For illustration purposes, we can view CAMD as two separate challenges, where the first one uses structural constraints to design reasonable molecules (including all symmetric solutions), and the second one uses symmetry-breaking
constraints to remove symmetric solutions. Note that the diversity of solutions will not change after breaking symmetry, because each molecule corresponds to least one solution (guaranteed by Section 2.3)
3.2 Counting feasible structures: performance of symmetry-breakinn
We choose two datasets QM7 104 105 and QM9 106 107 from CAMD literature to test the proposed methods. See Appendix C. for more information about QM7 and QM9. To test the efficiency of our sammetry-breakinne techniques, we build MIP formulation for CAMD and count all feasible structures for different N. By setting PoolScarchMode=2,, PoolSolutions=100, Gurobi can find many (up to 10º) feasible solutions to fill in the solution pool.
Table shows the performance of our symmetry -breaking constraints (S2) and comparing the baseline of (ST). Without adding (ST), we need to count the number of any graph with compatible
MiN-12 This number will be much larger features. Even ignoring features, the baseline would be 2 --
after introducing features, which loses the meaning as a baseline, so we use (SI) as the baseline.
Table 1:
Numbers of feasible solutions. The time limit is 48 hours. least 2.5 10 solutions are found for each time out (t.o.). For each dataset, the last column reports the percentage of removed symmetric solutions after adding (S2) and (S3) to the baseline of (ST). Higher percentage means breaking more symmetries.
QM7
QM9
N
SI
SI
S2
ST
SS3
(%)
ST
ST
(22)
SI
S3
(%)
17
10
10
41
15
9
40
112
37
175
54
69
3, 323
726
416
87
4, 536
1, 077
631
86
67, 020
11, 747
3, 003
117, 188
21, 441
5, 860
95
6
t.o.
443, 757
50, 951
2 98
t.o.
527, 816
59, 492
2 98
t.o.
t.o.
504, 952
t.o.
t.o.
776, 567
*
t.o. t.o.
t.o.
t.o.
t.o.
t.o.
3.3 Optimizing over trained GNNs for molecular design
For each dataset, we train a GNN with two GraphSAGE layers followed by an pooling layer, and three dense layers. Details about their structures and training process are shown in Appendix C.2 For statistical consideration, we train 10 models with different random seeds and use the models with smaller losses for optimization. Given N (E (4,5, 6, 7, 8}), and formulation (E (bilinear, bilincar+BrS. big-M, big-M+BrSJ), where "+BrS" means adding
brreetry-breaking constraints (C26) and (C27). we solve the corresponding optimization problem 10 times with different random seeds Gurobi. This means that there are 50 runs for each N and formulation.
Figure 3 shows the significant improvement of symmetry-breakinne time, big-M performs better. For the bilinear constraints, Gurobi 10.0.1 appears transform bilinear constraints into linear constraints. Appendix C shows that, after Gurobi's presolve stage, the big-M formulation has more continuous variables but fewer binary variables compared to the bilinear formulation. The fewer binary variables after presolve may explain the better big-M M performance.
Considering the average solving

---

47
12
bilinear big-M
bilinear+Br5 big BrS
20
time (x10 5
49
time (x102 S)
37
15
29
47
49
5
10
0 N =
50
N
= N = 6
N 4
= ==
N
Formulations
Formulations
Figure 3: These graphs report the time limeral achieve relative MIP optimality gap 10-4 averaged over the number of successful runs among 50 runs. We plot a bar for each formulation at different values of N. The left- and right-hand graphs correspond to datasets QM7 and QM9, respectively. The number of runs that achieve optimality is shown above each bar. Results with fewer than 20 successful runs do not have bar because so few runs terminated. The time limit is 10 hours.
In addition to the average solving time, we compare the performance of two formulations with breaking symmetry each run. In 341 runs out of 500 runs (i.e., 50 runs for each dataset {QM7, }M9} with each N {4, 5, 6, 7, 8))) the big-M formulation achieves optimality faster than the bilinear formulation. Additionally, Gurobi uses much of the solving time to improve the bounds. Table 6 and in Appendix C.3 report timeopt that denotes the first time to find the optimal solution. In 292 runs, the big-M formulation finds the optimal solution earlier than the bilinear formulation.
Discussion & Conclusion
We introduce optimization over trained GNNs and propose constraints break symmetry. We prove that there exists at least one indexing (resulting from Algorithm [ satisfying these constraints. Numerical results show the significant improvement after breaking symmetry. These constraints not limited to the problem (i.e., optimizing trained GNNs), technique (i.e., MIP), and application (i.e., CAMD) used in this work. For example, one can incorporate them into genetic algorithms instead of MIP, or replace the GNN by artificially-desinned score functions or other ML models. In other graph-based decision-making problems, as long as the symmetry issue caused by graph isomorphism exists, these constraints could be used to break symmetry. Moreover, the proposed frameworks for building MIP for GNNs as well as CAMD provide generality and flexibility for more problems.
Limitations. Assuming constant weights and biases excludes some GNN types. This limitation is an artifact of our framework and is possible to build MIP formulations for some other architectures. Using MIP may limit GNN size, for instance, edge features may enlarge the optimization problem.
Future work. One direction is make the MIP-GNN framework more general, such as adding edge features, supporting more GNN architectures, and developing more efficient formulations. Another direction is using the symmetry -beekkiggb constraints in graph-based applications beyond trained GNNs because we can consider any function that invariant to graph isomorphism. To facilitate further CAMD applications, more features such aromacity and formal charge could be incorporated. Also, optimizing lead structures a defined chemical space potential to design large molecules.
Acknowledgements
This work was supported by the Engineering and Physical Sciences Research Council [grant numbers EP/W003317/1 and EP/T001577/11 Imperial College Hans Rausing PhD Scholarship to SZ, BASF/RAEng Research Chair Data-Drives Optimisation RM, and Imperial College Research Fellowship to CT.
We really appreciate the comments and suggestions from reviewers and meta-reviewers during the peer review process, which were very helpful to clarify and improve this paper. We also thank Guoliang Wang for the discussion about graph indexing.
10

---

References
[1] Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun. Graph neural networks: A review of methods and applications. Al Open, 1:57-81, 2020.
[2] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S. Yu Philip. A comprehensive survey graph neural networks. IEEE Transactions on Neural Networks and Learning Systems, 32(1):4-22) 2020.
[3] Ines Chami, Sami Abu-El-Haijp Bryan Perozzi, Christopher Ré, and Kevin Murphy. Machine learning on graphs: A model and comprehensive taxonomy. Journal of Machine Learning Research, 23(89): 1-64, 2022.
[4] Joan Bruna, Wojeiech Zaremba, Arthur Szlam, and Yann Lecun. Spectral networks and locally connected networks on graphs. In ICLR, 2014.
[5] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. In NIPS, 2016.
[6] Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In ICLR, 2017.
[7] Ruoyu Li, Sheng Wang. Feiyun Zhu, and Junzhou Huang. Adaptive graph convolutional neural networks, In AAAI, 2018.
[8] Chenyi Zhuang and Qiang Ma. Dual graph convolutional networks for graph-based semi- supervised classification. In www, 2018.
[9] Bingbing Xu, Huawei Shen, Qi Cao, Yunqi Qiu, and Xueqi Cheng. Graph wavelet neural network. In ICLR, 2019.
[10] David K. Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P. Adams. Convolutional networks on graphs for learning molecular fingerprints. In NIPS, 2015.
[11] James Atwood and Don Towsley. crffusionn-cnvllutionnn neural networks. In NIPS, 2016.
[12] Mathias Niepert, Mohamed Ahmed, and Konstantin Kutzkov. Learning convolutional neural networks for graphs. In ICML, 2016.
[13] Will Hamilton, Zhitao Ying, and Jure Leskovee. Inductive representation learning on large graphs. In NIPS, 2017.
[14] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. In ICLR, 2017.
[15] Federico Monti, Davide Boscaini, Jonathan Masci, Emanuele Rodola, Jan Svoboda, and Michael M. Bronstein. Geometric deep learning on graphs and manifolds using mixture model CNNs. In CVPR, 2017.
[16] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. In ICML, 2017.
[17] Hongyang Gao, Zhengyang Wang, and Shuiwang Ji. Large-scale learnable graph convolutional networks. In Knowledge Discovery & Data Mining, 2018.
[18] Jiani Zhang, Xingjian Shi, Junyuan Xie, Hao Ma, Irwin King, and Bit-Yan Yeung. GaAN: Gated attention networks for learning on large and spatiotemporal graphs. In UAI, 2018.
[19] Hiroyuki Shindo and Yuji Matsumoto. Gated graph recursive neural networks for molecular property prediction. arXiv preprint A9Xiv:I909.00 00259, 2019.
[20] Xiaofeng Wang, Zhen Li, Mingjian Jiang, Shuang Wang, Shugang Zhang, and Zhiqiang Wei. Molecule property prediction based on spatial graph embedding. Journal of Chemical Information and Modeling, 3979983877-3888, 2019.

---

[21] Youjun Xu, Jianfeng Pei, and Luhua Lai. Deep learning based regression and multiclass models for acute oral toxicity prediction with automatic chemical feature extraction. Journal of Chemical Information and Modeling, 57(1 57(11):6672-2683 2017.
[22] Michael Withnall, Edvard Lindelöf, Ola Engkvist, and Hongming Chen. Building attention and edge message passing neural networks for bioactivity and physical-chemical property prediction. Journal of Cheminformatics, 12 2020.
[23] Artur M. Schweidmmann, Jan G. Rittig, Andrea König, Martin Grohe, Alexander Mitsos, and Manuel Dahmen. Graph neural networks for prediction of fuel ignition quality. Energy & Fuels, 34(9): 100111195--144400 2020.
[24] Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton, and Jure Leskovec. Hierarchical graph representation learning with differentiable pooling, In NIP, 2018.
[25] Xiaotong Zhang, Han Liu, Qimai Li, and Xiao-Ming Wu. Attributed graph clustering via adaptive graph convolution. In IJCAI, 2019.
[26] Hao Peng. Jianxin Li, Yu He, Yaopeng Liu, Mengjiao Bao, Lihong Wang, Yangqiu Song, and Qiang Yang. Large-scale hierarchical text classification with recursively regularized deep graph-CNN. In WWW, 2018.
[27] Liang Yao, Chengsheng Mao, and Yuan Luo. Graph convolutional networks for text classifica- tion. In AAAI, 2019.
[28] Qitian Wu, Hengrui Zhang, Xiaofeng Gao, Peng He, Paul Weng, Han Gao, and Guihai Chen. Dual graph attention networks for deep latent representation of multifaceted social effects in recommender systems. In WwW, 2019.
[29] Wenqi Fan, Yao Ma, Qing Li, Yuan He, Eric Zhao, Jiliang Tang, and Dawei Yin. Graph neural networks for social recommendation. In WWW, 2019.
[30] Rafiqul Gani. Computer-aided methods and tools for chemical product design. Chemical Engineering Research and Design, 82(11): 22(11):1494 1504, 2004.
[31] Lik Yin Ng, Fah Keen Chong, and Nishanth G. Chemmangattuvaalppii. Challenges and opportunities computer aided molecular design. Computer Aided Chemical Engineering. 34:25-34, 2014.
[32] Nick D. Austin, Nikolaos V. Sahinidis, and Daniel W. Trahan. Computer- aided molecular design: An introduction and review of tools, applications, and solution techniques. Chemical Engineering Research and Design, 116:2-26, 2016.
[33] Daniel C. Elton, Zois Boukouvalas, Mark D. Fuge, and Peter W. Chung. Deep learning for molecular design- - a review of the state of the art. Molecular Systems Design & Engineering. 4(4):828 849, 2019.
[34] Abdulelah S. Alshehri, Rafiqul Gani, and Fengqi You. Deep learning and knowledge- based methods for computer-aided molecular design toward a unified approach: State-of-the-an and future directions. Computers & Chemical Engineering, 141:107005, 2020.
[35] Faezeh Faez, Yassaman Ommi, Mahdieh Soleymani Baghshah, and Hamid R. Rabiee. Deep graph generators: survey. IEEE Access, 9: 9:106675-10600,, 2021.
[36] Wenhao Gao, Tianfan Fu, Jimeng Sun, and Connor W. Coley. Sample efficiency matters: benchmark for practical molecular optimization. In NIPS Truck Datasets and Benchmarks, 2022.
[37] Xiaolin Xia, Jianxing Hu, Yanxing Wang, Liangren Zhang, and Zhenming Liu. Graph-based generative models for de Novo drug design. Drug Discovery Today: Technologies, 32:45-53, 2019.
12

---

[38] Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez Timothy Hopper, Brian Kelley, Miriam Mathea, al. Analyzing learned molecular representations for property prediction. Journal of Chemical Information and Modeling, 59(8)3337033888 2019.
[39] Jiacheng Xiong, Zhaoping Xiong, Kaixian Chen, Hualiang Jiang, and Mingyue Zheng. Graph neural networks for automated de novo drug design. Drug Discovery Today, 2021.
[40] Jan G. Rittig. Martin Ritzert, Artur M. Schweidtmann, Stefanie Winkler, Jana M. Weber, Philipp Morsch, Karl Alexander Heufer, Martin Grohe, Alexander Mitsos, and Manuel Dah- men. Graph machine learning for design of high-octane fuels. AIChE Journal, page e17971, 2022.
[41] Christian Szegedy, Wojeiech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus. Intriguing properties of neural networks. In ICLR, 2014.
[42] Ga Wu, Buser Say, and Scott Sanner. Scalable planning with deep neural network learned transition models. Journal of Artificial Intelligence Research, 68:571-606, 2020.
[43] Rudy R. Bunel, Oliver Hinder, Srinadh Bhojanapalli, and Krishnamurthy Dvijotham. An efficient nonconvex reformulation of stagewise convex optimization problems. NIPS, 2020.
[44] Blanka Horvath, Aitor Muguruza, and Mehdi Tomas. Deep learning volatility: deep neural network perspective pricing and calibration in (rough) volatility models. Quantitative Finance, 21(1):11-27, 2021.
[45] Francesco Croce, Maksym Andriushchenko, and Matthias Hein. Provable robustness of relu networks via maximization of linear regions. In AISTATS, 2019.
[46] Ross Anderson, Joey Huchette, Will Ma, Christian Tjandraatmadja, and Juan Pablo Vielma. Strong mixed-integer programming formulations for trained neural networks. Mathematical Programmingg 183(1-))3-39, 2020.
[47] Matteo Fischetti and Jason Jo. Deep neural networks and mixed integer linear optimization. Constraints, 23(3):296-309, 2018.
[48] Alessio Lomuscio and Lalit Maganti. An approach to reachability analysis for feed-forward relu neural networks. arXiv preprint arXiv: 170.:07366.7751,.0 2017.
[49] Calvin Tsay, Jan Kronqvist, Alexander Thebelt, and Ruth Misener. Partition- formula- tions for mixed-integer optimization of trained ReLU neural networks. In NIPS, 2021.
[50] Alessandro De Palma, Harkirat S. Behl, Rudy Bunel, Philip Torr, and M. Pawan Kumar. Scaling the convex barrier with active sets. In ICLR, 2021.
[51] Keliang Wang, Leonardo Lozano, Carlos Cardonha, and David Bergman. Optimizing over an ensemble of trained neural networks. INFORMS Journal on Computing, 2023.
[52] Velibor V. MiKić. Optimization of tree ensembles. Operations Research, 68(5):1605-1624 2020.
[53] Miten Mistry, Dimitrios Letsios, Gerhard Krennrich, Robert M. Lee, and Ruth Misener. Mixed- integer convex nonlinear optimization with gradient-booste trees embedded. INFORMS Journal on Computing, 33(3): 103-11 1119, 19, 2021.
[54] Alexander Thebelt, Jan Kronqvist, Miten Mistry, Robert M. Lee, Nathan Sudermann-Merx, and Ruth Misener. Entmoot: framework for optimization over ensemble tree models. Computers & Chemical Engineering, 151:107343, 2021.
[55] Joey Huchette, Gonzalo Muñoz, Thiago Serra, and Calvin Tsay. When deep learning meets polyhedral theory: A survey, arXiv preprint arXiv:2305.002. 2023.
[56] Chih-Hong Cheng, Georg Nührenberg, and Harald Ruess. Maximum resilience of artificial neural networks. In Automated Technology for Verification and Analysis, 2017.
13

---

[57] Rudy R. Bunel, Ilker Turkaslan, Philip Torr, Pushmeet Kohli, and Pawan K. Mudigonda. A unified view of piecewise linear neural network verification. NIPS, 2018.
[58] Huan Zhang, Tsui-Wei Weng, Pin- Yu Chen, Cho-Jui Hsieh, and Luca Daniel. Efficient neural network robustness certification with general activation functions. In NIPS, 2018.
[59] Vincent Tjeng, Kai Y. Xiao, and Russ Tedrake. Evaluating robustness of neural networks with mixed integer programming. In ICLR, 2019.
[60] Elena Botoeva, Panagiotis Kouvaros, Jan Kronqvist, Alessio Lomuscio, and Ruth Misener. Efficient verification of ReLU-based neural networks via dependency analysis. In AAAI, 2020.
[61] Rudy Bunel, P Mudigonda, Ilker Turkaslan, Philip Torr, Jingyue Lu, and Pushmeet Kohli. Branch and bound for piecewise linear neural network verification. Journal of Machine Learning Research, 21(2020). 2020.
[62] Rudy Bunel, Alessandro De Palma, Alban Desmaison, Krishnamurthy Dvijotham, Pushmeet Kohli, Philip Torr, and M. Pawan Kumar. Lagrangian decomposition for neural network verification. In UAI, 2020.
[63] Kaidi Xu, Huan Zhang, Shiqi Wang. Yihan Wang, Suman Jana, Xue Lin, and Cho-Jui Hsieh. Fast and complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers. In ICLR, 2021.
[64] Shiqi Wang, Huan Zhang, Kaidi Xu, Xue Lin, Suman Jana, Cho-Jui Hsieh, and Zico Kolter. Beta-crown: Efficient bound propagation with per-neuron split constraints for neural network robustness verification. In NIPS, 2021.
[65] Buser Say, Ga Wu, Yu Qing Zhou, and Scott Sanner. Nonlinear hybrid planning with deep net learned transition models and mixed-integer linear programming. In IJCAI, 2017.
[66] Arthur Delarue, Ross Anderson, and Christian Tjandraatmadja. Reinforcement learning with combinatorial actions: An application to vehicle routing. NIPS, 2020.
[67] Moonkyung Ryu, Yinlam Chow, Ross Anderson, Christian Tjandraatmadja, and Craig Boutilier. CAQL: Continuous action Q-learning. In ICLR, 2020.
[68] Thiago Serra, Xin Yu, Abhinav Kumar, and Srikumar Ramalingam. Scaling up exact neural network compression by ReLU stability. In NIPS, 2021.
[69] Theodore P. Papalexopoulos, Christian Tjandraatmadja, Ross Anderson, Juan Pablo Vielma, and David Belanger. Constrained discrete black-box optimization using mixed-intege pro- gramming. In ICML, 2022.
[70] Haoze Wu, Clark Barrett, Mahmood Sharif, Nina Narodytska, and Gagandeep Singh. Scal- able verification of GNN-based job schedulers. Proceedings of the ACM on Programming Languages, 6(0OPSLA2): 1065-1065, 2022.
[71] Tom Mc Donald. Mixed integer (non-) linear programming formulations of graph neural networks. Master Thesis, 2022.
[72] Eric J. Friedman. Fundamental domains for integer programs with symmetries. In Combinato- rial Optimization and Applications, 2007.
[73] Leo Liberti. Automatic generation of symmetry -reeekin constraints. Lecture Notes in Computer Science, 5165:328-333 2008
[74] François Margot. Symmetry in integer linear programming. 50 Years of Integer Programming 1958-2008: From the Early Years to the State-of-the-Arr, 2009.
[75] Leo Liberti. Reformulations in mathematical programming: Automatic symmetry detection and exploitation. Mathematical Programming. 131 -044 2012.
[76] Leo Liberti and James Ostrowski. Stabilizer-hased symmetry breaking constraints for mathe- matical programs. Journal of Global Optimization, 60:183-194, 2014.
14

---

[77] Gustavo Dias and Leo Liberti. Orbital independence in symmetric mathematical programs. In Combinatorial Optimization and Applications, 2015.
[78] Christopher Hojny and Mare E. Pfetsch. Polytopes associated with symmetry handling. Mathematical Programming, 175(1-2):197-240, 2019.
[79] Francesco Ceccon, Jordan Jalving, Joshua Haddad, Alexander Thebelt, Calvin Tsay, Carl D. Laird, and Ruth Misener. OMLT: Optimization & machine learning toolkit. Journal of Machine Learning Research, 23(1):15822-15836 2022.
[80] Alan Frisch, Brahim Hnich, Zeynep Kiziltan, Ian Miguel, and Toby Walsh. Global constraints for lexicographic orderings. In Principles and Practice of Constraint Programming- CP 2002, pages 93-108, 2002.
[81] Joao Parques-Silvall Josep Argelich, Ana Graça, and Inés Lynce. Boolean lexicographic optimization: algorithms & applications. Annals of Mathematics and Artificial Intelligence, 62:317-343, 2011.
[82] Hani A. Elgabou and Alan M. Frisch. Encoding the lexicographic ordering constraint in sat modulo theories. In Thirteenth International Workshop on Constraint Modelling and Reformulation, 2014.
[83] Guy Katz, Clark Barrett, David Dill, Kyle Julian, and Mykel J. Kochenderfer. Reluplex: An efficient SMT solver for verifying deep neural networks. In Computer Aided Verification, 2017.
[84] Tobias Achterberg. SCIP: solving constraint integer programs. Mathematical Programming Computation, 1:1 1:1-41, 2009.
[85] Tommi Junttila and Petteri Kaski. Engineering an efficient canonical labeling tool for large and sparse graphs. In ALENEX, 2007.
[86] Tommi Junttila and Petteri Kaski. Conflict propagation and component recursion for canonical labeling. In TAPAS, 2011.
[87] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2023. URL https:tp://// htttss:/wwtt /www gurobi.com
[88] Volker Kaibel, Matthias Peinhardt, and Marc E. Pfetsch. Orbitopal fixing. Discrete Optimiza- tion, 8(4):595-610, 2011.
[89] François Margot. Pruning by isomorphism in branch-and-cut. Mathematical Programming, 94:71-90, 2002.
[90] Hanif D. Sherali and Cole Smith. Improving discrete model representations via symmetry considerations. Management Science, (1110)::396-140)) 2001.
[91] Ahmed Ghoniem and Hanif D. Sherali. Defeating symmetry in combinatorial optimization via objective perturbations and hierarchical constraints. IIE Transactions, 43(8):575-588, 2011.
[92] Volker Kaibel and Mare Pfetsch. Packing and partitioning orbitopes. Mathematical Program- ming. 114(1):1-36 2008
[93] Yuri Faenza and Volker Kaibel. Extended formulations for packing and partitioning orbitopes. Mathematics of Operations Research, 34(3):686-69) 697, 2009.
[94] James Ostrowski, Miguel F. Anjos, and Anthony Vannelli. Symmetry in scheduling problems. Citeseer, 2010.
[95] Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR 2019 Workshop on Representation Learning Graphs and Manifolds, 2019.
[96] David Bergman, Teng Huang, Philip Brooks, Andrea Lodi, and Arvind U. Raghunathan. Janos: An integrated predictive and prescriptive modeling framework. INFORMS Journal on Computing, 34(2):807-816, 2022.
15

---

[97] Laurens Lueg, Bjarne Grimstad, Alexander Mitsos, and Artur M. Schweidtmann. reluMIP: Open source tool for MILP optimization ReLU neural networks, 2021. URL https Igittuu..o com/ ChweehngAI/ReeLUuuaC. ANN MILP
[98] O. Odele and Sandro Macchietto. Computer aided molecular design: A novel method for optimal solvent selection. Fluid Phase Equilibria, 82:47-54, 1993.
[99] Nachiket Churi and Luke E.K. Achenie. Novel mathematical programming model for computer aided molecular design. Industrial & Engineering Chemistry Research, 35(10):3788-3777 3794, 1996.
[100] Kyle V. Camarda and Costas D. Maranas. Optimization in polymer design using connectivity indices. Industrial & Engineering Chemistry Research, 38(5):1884-189,, 1999.
[101] Manish Sinha, Luke E.K. Achenie, and Gennadi M. Ostrovsky. Environmentally benign solvent design by global optimization. Computers & Chemical Engineering, 23(10):1381-133 1999.
[102] Nikolaos V. Sahinidis, Mohit Tawarmalani, and Minrui Yu. Design of alternative refrigerants via global optimization. AIChE Journal, 4977):1711-11775 2003.
[103] Lei Zhang, Stefano Cignitti, and Rafiqul Gani. Generic mathematical programming formulation and solution for computer-aided molecular design. Computers & Chemical Engineering, 78:
79-84, 2015.
[104] Lorenz C. Blum and Jean-Louis Reymond. 970 million druglike small molecules for virtual screening in the chemical universe database gdb-13. 13. Journal of the American Chemical Society, 131(25):8732 8733, 2009.
[105] Matthias Rupp. Alexandre Tkatchenko, Klaus-Robert Müller, and O. Anatole Von Lilienfeld. Fast and accurate modeling molecular atomization energies with machine learning. Physical Review Letters, 108(5):058301, 2012.
[106] Lars Ruddigkeit, Ruud Van Deursen, Lorenz C. Blum, and Jean-Louis Reymond. Enumeration 166 billion organic small molecules in the chemical universe database gdb- 17. Journal of Chemical Information and Modeling, 52(1 1):2864 2875, 2012.
[107) Raghunathan Ramakrishnan, Pavlo O. Dral, Matthias Rupp, and O. Anatole Von Lilienfeld. Quantum chemistry structures and properties of 134 kilo molecules. Scientific Data, 1(1):1-7, 2014.
16

---

A Theoretical guarantee: supplementary
A.1 Proofs properties and lemma
< 8, Yu E Vi S, Vu E V/2
Property 1. For any s 1,2, ..... 1, (v)
Proof. At the s-th iteration, nodes in Vi" have been indexed by 0, 0.1....,5 8 - 1. Therefore, if v E Vi". then (v) I(v) s. I V22, since (*) is the sum of and rank (v) (which is non-negative), then we have (u) 2 s.
Property 2. For any $1 82
1
S - No (v) (v) [s1]. Ve V/22
Proof. First, N* is well-defined on V since V C V By definitions of N* and Vf, for any V222 we have:
(v) [[[(w) [ n }
[I(u) VV"
[I(u) uE (N(
[[[()) [E) {1
[I((u) (EN(()))
(v)
where the second equation the fact that V C V13, and the fourth equation holds since " E V + I(u)= < 81 (using Property D.
Property 3. Given any multisets A, B with no more than L integer elements in |0, M 1], we have:
LO(A) & LO(B)= = LO(An [m]) s LO(B n [m]), Vm = 1,2, M
Proof. By the definition lexicographical order for multisets, denote the corresponding sequence to A, An [m], B n [m]] by a, am 6tdh (SM, L). respectively. Then it equivalent show that:
LO(a) LO(b) LO(a")) z LO(b"), Vm 1.2. M
Let a = (a1, .2.............................. a1) and b (bı, by...... b.). If (a) = (b))). then b and at = b" Thus LO(a")) Lb Otherwise, if LO(a) < LO(b), then there exists < 1 S L such that:
a: bi
l1 vi<i
a1 bi
m S a,, then at = 1 which means LO(a")) LO(b"). Otherwise, if m >> dr, then am and fm share the same first 1 - elements. But the 1-th element of am is at, while the 1-th element of b., either by or M. In both cases, we have LO(a")) < LO(D").
Lemma 1. For any two nodes u and v,
I(u) < I(v) II (T) Vs =,2,,
Proof. Case 1: IfI(u) < I(v) < s, then:
(") = I(u) I(u) (v)
Case 2: I()) < I(p), then [²( I[u)) and:
I2( 2 > I(u)= = (w)
where Property D is used.
17

---

Case 3: If S I(u) < I(e), at I(u) ([u iteration, vi is chosen to be indexed Z(u). Thus:
LO ( (u)) IO (v))
According to the definition of (2 and Ne (2 (.), we have:
(T()) N = =, (u) [[u)]] N N (v) N₂ (v) (v) [T()]
Using Property (with A = N1 ( (u), B = Ne 2 (v), m = I(u)) yields:
10( & LO
Apply Property 2 (with $1 S I(u)) for and e, we have:
= (u) [s]. (²)()) (n
Using Property aggain (with A NZ(") B= N m 8) gives:
LO(N(()) S lo()
Recall the definition of I*)),, we have:
(u) S (u)
Lemma 2. For any undirected, connected graph G = (V, if one indexing of G sastifies S3 then it satisfies SI.
Proof. Node itself is a connected graph. Assume that the subgraph induced by nodes {0, 10.1...... m} connected, it suffices to show that the subgraph induced by nodes {0, I,..., 1} is connected. Equivalently, we need to prove that there exists u < v + such that A,, =
Assume that Au = 0, Vu < x + 1. Since G is connected, there exists v' > + 1 such that:
Eu 1, s.l. Au,u' =1
Then we know:
N(c + - N(')) IV + +1
Recall the definition of LO(+), we obtain that:
(0(N(++1)nB) + 1)) LO(N(U') n 1I)
Since the indexing satisfies S3 then have:
LO(Nww) + +1)) s LO(N(* +))[w]), <<w<d
Applying Property 3
LO(NNw)){W ()) < LO( 1J), W < <w<<<< <d 1))
Note that w v. Therefore,
LO((w)) + 1I) L) Ve < w<d
Choosing - gives:
LO(N(+1))lo+ 1)) < LO((('))++1))
The contradiction between and completes the proof.
18

---

A.2 Example for applying Algorithm
Given a graph with N nodes as shown in Figure 4 where up already indexed 0. We provide
details of Algorithm [ by indexing the rest 5 nodes step by step for further illustration.
Before indexing, we first calculate the neighbor sets for each node:
((v)) {ul, U2, vs, 04, (((1) {vo, v2, U3, v4),N((1))) {vo, vi, 05
((3)) = {vo, U1, 24., {vo, v1, ca)). = feo, 02)
=1 1/ {vo). V {ui 22 23, u5
00
11
22
25
Figure 4: V1 {oo}
Obtain indexed neighbors each unindexed node:
(v1)) (0), {0). = (v)) = (0), (0)
Rank all unindexed nodes:
rank(w1) rank(v2) rank(us) rank(va) rank(us)
Assign a temporary index to each node based on previous indexes (for indexed nodes) and ranks (for unindexed nodes):
(vo) 0, (.)
1,,
After having indexes for all nodes, define temporary neighbor sets:
(V (m1) (0. 1, ()) = (o.) (0.1,1) (0,1,1),
Ny (v4) = (0, (o, = (0
Based on the temporary neighbor sets, 01 is chosen to be indexed 1 (i.e., I(mi = 1).
2: V (vo, J. V2 (12, 23, 04a Us}
0/
14
Us
Figure 5: V2 (uo,
19

---

Obtain indexed neighbors to each unindexed node:
(12) (0.1).A = (0. (0.1). = (0,1): (0)
Rank all unindexed nodes:
rank(up) = ran)(va) = rank(ug) = 0, rank(us) =
Assign temporary index to each node based on previous indexes (for indexed nodes) and ranks (for unindexed nodes):
(2,
=
After having indexes for all nodes, define temporary neighbor sets:
(12) (0,1,3). (0,1,)), = (0.1.)).
{0,2)
Based the temporary neighhor sets, both uy and 04 be chosen to be indexed 2. Without loss of generality, let I(vz) = 2.
Note:
This step explains why temporary indexes and neighbor sets should be added into Algorithm Otherwise, V/2 is also valid to be index 2, following which vg 14, Us will be indexed 3, 5. Then
I
the neighbor set for I(w2) 2 is {0, 1, 1,5) while the neighbor set for I(vs) = 3 is {0, 1, 4], which violates constraints
S3
3 V/²
{vo, v1, v3), V22
{u2, 14, v5
0
010
2
Figure 6: Vi = {wo, U1 03}
Obtain indexed neighbors to each unindexed node:
No (22) = 10. (w. {0)
Rank all unindexed nodes:
rank(u2) = 1, rank(us) = 0, rank(us) 2
Assign temporary index to each node based previous indexes (for indexed nodes) and ranks (for unindexed nodes):
= (,. = =5
After having indexes for all nodes, define temporary neighbor sets:
N (v2) = 10, (0.1.)), = {0, 4]
Based on the temporary neighbor sets, va is chosen to be indexed 3 (i.e., I(va) = 3).
4 V {vo, U1, 23, 24 Vg!
{12, v5}
20

---

2
3
F3
1/
16
Figure Vi (vo, v1, 03; v24
Obtain indexed neighbors to each unindexed node:
(12) {0, [0.1).. {0}
Rank all unindexed nodes:
rank(v2) = 0, rank(us)=
Assign temporary index to each node based on previous indexes (for indexed nodes) and ranks (for unindexed nodes);
I' = 4, 4,
After having indexes for all nodes, define temporary neighbor sets:
(12) {0, 1, LV {0, Based on the temporary neighbor sets, 12 is chosen to be indexed 4 (i.e., I(v2) = 4).
5 V
(vo, U1, 22, va). V22
{us)
up
2
03
Figure 8: V/$ {vo, U1, 82, U3 U1]
Since there only node U9 unindexed, without running the algorithm we still know that vy chosen to be indexed 5 (i.e., I(us) = 5).
Output:
U1
Ua
U/
Figure 9: V/
{vo, U1, U2, V3, U4, U5
21

---

B MIP formulation for CAMD: details
In this part, we provide details of the MIP formulation for CAMD, including parameters, variables and constraints. Especially, we use atom set NC 0, as an example to help illuminate. The number nodes N exclude hydrogen atoms since they are implicitly represented by N features.
B.1 Parameters
Table 2: List of parameters
Parameter
Description
Value
N
number of nodes
F
number of features
16
N
number of atom types
4
N
number of neighbors
5
N
number of hydrogen
5
I
index for N index for N index for N
{0, 1,,1,233
I"
{4,5, 6, 7, 8}
Th
{9, 10, 11, 12, 13)
Idb
index for double bond
14
index for triple bond
15
Atom
atom types covalence of atom
{C, N,O,S} 24,
Cov
B.2
Variables
Table 3: List variables for atom features
X
Type
X
#Neighbors
X
#Hydrogen
0
C
4
9
0
5
10
2
0
6
11
2
7
12
8
13
For the existence and types of bonds, we add the following extra features and variables:
X E [N]: if atom is included in at least one double bond.
N v E [N]: if atom v included at least one triple bond.
An,os * [N]: if there is one bond between atom and
Av.,, E [N]: if node exists.
.DBa,v, # [N]: if there is one double bond between atom u and e.
. .B..,, u # v E [N]: if there is one triple bond between atom and
design molecules with most N atoms, which brings flexibility in real-world applications. In our experiments, however, we fix them to consider possible structures with exact N atoms. Therefore, constraints
Remark: Variables A,, admits
replaced by:
A
1, Yu E [N]
Without this replacement, we can add an extra term on C26) to exclude non-existing nodes and only index existing node:
Σ
2F F-5-1 f-1.
Xo.5 5 Σ F-1-11
X + 2 (1 An,,), Ve [N]V(O)
SE(F]
FE(P]
22

---

B.3 Constraints
There are at least two atoms and one bond between them:
Ao,o = A1,1 = Au,1
(CI)
Atoms with smaller indexes exist:
Ae 2 Aut1,+t1 Vu E IN 1]
(C2)
The adjacency matrix A is symmetric:
Au, = Auju Vu, v E [N], u <
(C3)
Atom v linked to other atoms and only if exists:
(N 1) Au, 2 EAu,. Ve E [N]
(C4)
Note: N - here is big-M constant.
If atom v exists, then it has to be linked with at least one atom with smaller index:
A S
Σ
Augs Vo E [N]
(C5)
u<P
Atom u cannot have double bond with it self:
DB 0, Ve E [N]
(C6)
DB is also symmetric:
DB = DB, Yu, E [N],
(C7)
Atom v cannot have triple bond with it self:
TB O, [NE
(C8)
TB is also symmetric:
TB TBaus Vu, v [N], <
(C9)
There at most one bond type between two atoms:
DB + TBu,e S Aa,es Vu, E [N]. u<v
(C10)
If atom v exists, then it has only one type:
Au,
Exes Ve E [N]
(C11)
fEl
If atom e exists, then it has only number of neighbors:
A
Σ
Xus, Vo E [N]
(C12)
SEL
If atom e exists, then it has only one number of hydrogen:
A
Σ
X
Vu [N]
(C13)
FEI
The number of neighbors of atom u equals to its degree:
A
Σ i N.,,, Ve (N)
(C14)
SE[Na]
If there is possibly one double bond between atom u and v, then the double bond feature of both atoms should be 1 and these two atoms are linked:
3 DB S SXn,jan + Xu + Au,. VU,UE [N], u <
(C15)
23

---

If there is possibly one triple bond between atom and v, then the triple bond feature of both atoms should be 1 and these two atoms are linked:
3 T T N,,, + x,, + Au,, Vu, E [N],
(C16)
The maximal number of double bonds for atom is limited by its covalence:
Cou;
Σ
DB
S
Σ
Xx, Vo [N]
(C17)
2
WE[N]
(E[N9]
The maximal number of triple bonds for atom v is limited by its covalence:
Con
Σ TB S
Σ
Xx,t Vo E [N]
(C18)
3
ME[N]
&E|N"
If atom is not linked with other atoms with double bond, then its double bond feature should be
Σ
DBu,es Ve [N]
(C19)
WE[N]
If atom e is not linked with other atoms with triple bond, then its triple bond feature should be 0:
N
Σ TBu,, Vo E [N]
(C20)
WE[N]
The covalence of atom equals the sum of its neighbors, hydrogen, double bonds, and triple bonds:
Cou.
X
Σ
X
X
+
Σ
xE[N"
iE(N")
vE(NA]
(C21)
Σ
DB
+
Σ
T.
2.
Ve E [N]
NE[N]
NE[N]
Practically, there could be bounds for each type of atom:
Atom. S Σ
X S UUBttom Vi E [N9]
(C22)
DE[N]
Practically, there could be bounds for number of double bonds:
LB S Σ 5 DB. S UB
(C23)
VEIN] u5C
Practically, there could be bounds for number of triple bonds:
LB S Σ S BB
(C24)
VE(N
Practically, there could be bounds for number of rings. Since force the connectivity, bonds are needed to form a tree and then each extra bond forms ring. Therefore, can bound the number of rings by:
LBring S
Σ N
Σ
An,o - (N 1) S UUB Bring
(C25)
u<w
24

---

Implementation details and more results
C.1 Dataset preparation
Both datasets used this paper are available in python package chemprop Molecular features (133 atom features and 14 bond features) are extracted from their SMILES representations. Based on these features and our definition for variables Appendix [.2] we calculate F 16 (since both datasets have types of atoms) features for each atom, construct adjacency matrix A as well as double and triple bond matrices DB and B. In this part, we first briefly introduce these datasets and then provide details about how preprocess them.
QM7 dataset: Consist of 6, 830 molecules with most heavy atoms C, N, O, One quantum mechanics property is provided for each molecule. We filter the dataset and remove three types of incompatible molecules: (1) molecules with aromaticity (i.e., with float electron(s) inside ring); (2) molecules constraints S atom with covalence 6 (since one type of atom can only have one covalence in our formulation); (3) molecules with non-zero radical electrons. After preprocessing, there are 5, 822 molecules left with properties re-scaled to 0 ~ 1.
To calculate the bounds needed in constraints (C22) The maximal ratios between the number of N.O and S atoms and the total atoms are 3/7, 1/3, 1/7,, respectively, while the minimal ratio between the number of C atoms and the total atoms is 1/2. The maximal ratios between the number of double bonds, triple bonds, rings and the number of atoms are 3/7,3/7, 3/7, respectively. Therefore, we set the following bounds:
we check all molecules in QM7 dataset.
LBc
N
UBN max
1,
3N )
UB max (1.
(113))
N AT32 ) .UBs max
N 7
2
7
N
N 2
N
UBa
UB
U Bring
2
2
QM9 dataset: Consist of 133, 885 molecules with at most 9 heavy atoms C, N, O. 12 quantum mechanics properties are provided for each molecule. We filter the dataset and remove two types of incompatible molecules: (1) molecules with aromaticity; (2) molecules with non-zero radical electrons. After preprocessing, there are 1 108, 723 molecules left. Among these features, we choose one property Ug. which denotes the internal energy OK. These properties are re-scaled 0
Similarly, we set the following bounds for the experiments on the QM9 dataset:
N
3N 5
4N
4N UBr
LBc
UBx
UBo
5
7
5
N
N
2N
UB
.BB
.ring
2
2
C.2 Training of GNNS
For QM7 dataset, the GNN's structure is shown in Table 4 For each run, we first randomly shuffle the dataset, then use the first 5000 elements to train and the last 822 elements to test. Each model is trained for 100 iterations with learning rate 0.01 and batch size 64. The average training (41 loss is 0.0241 and test loss is 0.0285.
Table Structure of GNN for QM7 dataset
index
layer
in features
out features
activation
SAGEConv
16
16
ReLU
SAGEConv
16
32
ReLU
Pooling
32NV
32
None
Linear
32
16
ReLU
Linear
16
4
ReLU
Linear
4
1
None
25

---

For QM9 dataset, the GNN's structure is shown in Table We apply the same experimental set-up before, but with larger number of channels to this much larger dataset. For each run, after randomly shuffled, the first 80000 elements are used to train and the last 28723 elements are used to test. The average training loss is 0.0036 and test loss is 0.0036.
Table 5: Structure of GNN for QM9 dataset
index
layer
in features
out features
activation
SAGEConv
16
32
ReLU
SAGEConv
32
64
ReLU
Pooling
64.V
64
None
Linear
64
16
ReLU
Linear
16
4
ReLU
6
Linear
4
1
None
For both datasets, the trained GNNs achieve comparable accuracy to the state-of-the-art methods 38 Note that this comparison is unfair since we preprocess each dataset and remove some incompatible molecules. Also, we only choose one property train and then optimize. However, low losses mean that the trained GNNs approximate true properties well.
C.3 Full experimental results
We report full experimental results in this section. For each dataset with given N and formulation, numbers of variables (#car:: continuous, #vary: binary) and constraints (#con) after presolve stage in Gurobi are first reported. Then count the number of successful runs (i.e., achieving optimality) #run. For all successful runs, we provide the mean, 25th percentile (Q1), and 75th percentilie (Qs) of the running time timezot well as the first time to find the optimal solution timeopt- The time limit for each run is 10 hours. Except for the time limit and random seed, we use the default setting in Gurobi for other parameters such the tolerance.
Table 6: Numerical results for QM7 dataset.
limetol (s)
Limeopt (s)
method
#vare
#varn
#com
#run
mean
Q
Q3
mean
Q
Q3
bilinear
616
411
1812
50
477
314
579
101
13
147
bilinear+BrS
616
383
1755
50
149
115
167
94
66
131
big-M
792
281
1838
50
514
206
722
130
150
big-M+BrS
760
269
1748
50
114
91
128
53
13
94
bilinear
840
611
2776
37
9365
3551
9029
471
169
339
bilinear+BrS
840
565
2679
50
659
282
652
198
136
257
big-M
1142
359
2795
50
6369
1424
9703
729
74
616
big-M+BrS
1097
346
2660
50
501
266
660
205
91
281
bilinear
1096
828
3893
23853
22081
25627
389
373
396
bilinear+Brr
1064
747
3602
48
6400
1748
6225
1751
334
2154
6
big-M
1628
436
4018
10
21484
16095
24848
14368
9966
19930
big-M+BrS
1505
420
3665
47
4700
1656
3759
1156
515
1440
bilinear
1384
1076
5193
*
*
*
bilinear+BrS
1320
950
4684
47
12989
4651
19229
3366
497
4587
7
big-M
2134
516
5364
*
*
*
big-M+BrS
1943
495
4806
10802
4159
16101
4008
64
5251
bilinear
1704
1355
6676
bilinear+BrS
1608
1181
5941
12
16219
7623
22951
3082
90
223
8
big-M
2676
599
6869
*
*
*
big- M+BBS
2390
572
6040
12
6606
2132
3120
2854
55
87
26

---

Table 7: Numerical results for QM9 dataset.
time (s)
time (s)
N method
# uare
#uarb
#oo
#rum
mean
Q1
Q3
mean
Q1
Q3
bilinear
1178
599
3050
50
368
252
353
240
186
279
bilinear+BrS
1181
580
3020
50
256
211
288
184
60
258
4
big-M
1340
453
2933
50
258
193
311
190
141
235
big-M+BrS
1316
441
2863
50
167
89
227
111
43
171
bilinear
1617
846
4563
48
5439
1108
7573
2604
585
1494
bilinear+BrS
1627
823
4545
50
1039
704
1248
628
579
744
5
big-M
1962
584
4480
50
3242
1111
3535
1642
757
1762
Big-M+BrS
1877
564
4260
50
1127
672
1336
640
397
868
bilinear
2118
1110
6319
21
16767
7598
25642
2969
479
1343
bilinear+BrS
2068
1049
6017
48
6907
3595
8691
2394
491
1937
6
big-M
2768
718
6562
29
12625
5573
18020
3626
314
2181
big-Mee big-M+BrS
2500
684
5829
50
1985
660
2613
1132
302
1243
bilinear
2681
1405
8348
3
22759
19543
24667
1417
995
1722
bilinear+BrS
2571
1308
7768
34
21649
16723
27748
2983
813
3808
7
big-M
3549
846
8649
6210
6210
6210
6108
6108
6108
big-M+BrS
3180
810
7613
49
6586
2946
9655
1427
564
1255
bilinear
3306
1731
10650
Filinear+BrS
3136
1597
9790
0
8
big-M
4420
977
11013
big-M+BrS
3922
935
9595
29
17864
11113
21431
3405
705
4308
27


Reviews:
[{"review_id": "WPFTeZn3eA", "reviewer": "Reviewer_eveB", "summary": "This paper presents novel symmetry-breaking constraints for optimizing trained graph neural networks (GNNs) and addresses the graph isomorphism issue. The authors develop two mixed-integer optimization formulations and evaluate their methods in the context of molecular design. ", "questions": "1. Do the indexing produced by Algorithm 1 for a connected graph satisfy S1? ", "limitations": "n/a ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. **Relevance and Application**: The paper addresses the important issue of GNN optimization and demonstrates practical utility through molecular design application.\n\n2. **Empirical Evidence**: The authors provide convincing experimental results supporting their claims, underscoring the proposed method's effectiveness.\n\n3. **Clarity**: The manuscript is well-written, presenting complex ideas with notable clarity.\n", "weaknesses": "1. **Insufficient exploration of S1 and S3 compatibility**: The paper falls short in exploring whether S1 and S3 requirements can coexist for certain indexings of a connected graph with features.  \n\n2. **Insufficient explanation of structure feasibility constraints**: There is a lack of clarity concerning the interplay between structure feasibility constraints C1 - C21 and the indexing process. Specifically, the paper does not elucidate the number of these constraints that are related to indexing and if they are compatible with S3.  \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Sincere thanks for these comments which are fair and help us clarify and improve our work.\n\n**Weakness 1 \\& Question 1 [Compatibility of S1 and S3]**\\\nThanks for mentioning this. Our original paper ignored the compatibility of S1 and S3, but the reviewer is correct that we should add this to the paper. S3 is tighter than S1 for connected graphs (as shown in Lemma 2). The reason we use both S1 and S3 is that S1 can force a connected graph.\n\n**Lemma 2** For any undirected, connected graph $G=(V,E)$, if an indexing satisfies (S3), then it satisfies (S1).\n\nWe briefly show their compatibility here and will add Lemma 2 with proof in the final paper. Given any connected graph with an indexing satisfying (S3), we can prove that this indexing satisfies S1 using mathematical induction: node $1$ has to be a neighbor of node $0$, node $2$ has to be a neighbor of node $0$ or $1$, and so on.\n\n**Weakness 2 [Compatibility between constraints]**\\\nWe appreciate this suggestion. We did not (but should) clarify the interplay between structural constraints and symmetry-breaking constraints. We clarify the compatibility here and will add this explanation in the final paper.\n\nStructural constraints C1 - C25 force graph structures and feasible features that are compatible with molecular design. To present a general formulation, C1 - C25 support controlling the existence of each node by adding variables $A_{v,v},v\\in [N]$, which means that we can design molecules with fewer than $N$ atoms. In our analysis and numerical experiments, however, we set $A_{v,v}=1,v\\in [N]$ for simplicity. Under such setting, only C5, which corresponds to S1, is related to graph indexing (C1 - C4 and C6 - 25 are irrelevant to indexing). Without assuming $A_{v,v}=1,v\\in [N]$, we can add an extra term on C26 ($2^F\\cdot (1-A_{v,v})$) to exclude non-existing nodes:\n$$\n\\begin{equation*}\n    \\begin{aligned}\n        \\sum\\limits_{f\\in [F]}2^{F-f-1}\\cdot X_{0,f}\\le \\sum\\limits_{f\\in [F]} 2^{F-f-1}\\cdot X_{v,f} + 2^F\\cdot (1-A_{v,v}),~\\forall v\\in [N]\\backslash\\{0\\}.\n    \\end{aligned}\n\\end{equation*}\n$$\nBy doing that, we are still considering the indexing for all existing nodes without being influenced by non-existing nodes.\n\nFor illustration purposes, we can view CAMD as two separate challenges, where the first one uses structural constraints to design reasonable molecules (including all symmetric solutions), and the second one removes symmetric solutions."}, {"review_id": "hyvHsqZCZY", "reviewer": "Reviewer_ni19", "summary": "This paper suggests using a symmetry breaking indexing for nodes in a graph. The argument is that for inverse problems, such as molecule design, graph isomorphism results in finding many equivalent solutions. The symmetry breaking indices are supposed to partially alleviate this redundancy. They use pre-trained GNN (based on GraphSage) for computer-aided molecule design (CAMD). During the design, they use Mixed-Integer Programming (MIP) (solved using Gurobi) to solve a set of constraints to make the molecules chemically realistic as well as constraints arising from their indexing scheme. \n\nThey conduct experiments on QM7 and QM9 and find that the symmetry breaking methods perform much better. They show more successful runs (i.e. producing molecules satisfying the constraints) and achieve this in a fraction of the time, compared to models without symmetry breaking. ", "questions": "1. Table 1: In each row, is higher necessarily better? If we want less symmetry copies, should sometimes less be better?  \n2. How do we know that the solutions found with the symmetry broken model are as diverse as the one with symmetries? These are small molecules. Can you use a fast Weisfeiler-Lehman type hash to check isomorphism of found molecules at least in small N experiments? Maybe also use a measure of diversity, some kind of entropy or something. \n \n### Minor questions:  \nline 75 “Directed to undirected” how?   \nline 106 A is now a multiset, not the adjacency matrix anymore?   \nline 183 (*) is not permutation invariant, as written. $w_{u\\to v}$ is like a generalized Graph Attention (GAT).   \nline 187 Dense NN: do you mean a fully connected layer (no graph)?  ", "limitations": "The authors discuss technical limitations, but I don't think there are immediate societal impacts, as the paper pertains to a technical method to speed up generative design in graph problems. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper makes a good case for how adding an ordering (indexing nodes in a specific manner) can dramatically improve inverse problems on graphs, e.g. designing molecules, solved using MIP. \nThe paper is well structured and discusses the literature in good depth. It explains the algorithm and experimental setup. The results are also very strong for the experiments.  ", "weaknesses": "The experiments show that with the indexing, MIP converges much faster. Table 1 states that MIP find many feasible solutions for each setting and set of constraints. One thing I hope the authors can assess (maybe using a Weisfeiler-Lehman hash or similar; see questions below) is how much the symmetry breaking *actually* removes isomorphic solutions, compared to the runs without symmetry breaking. they could do it for a small setting, maybe QM9, N=3 or 4. The goal would be two fold: \n1. To show that the symmetry breaking indeed removed many isomorphic molecules\n2. to check how the diversity of solutions in the symmetry broken case compares with the diversity of molecule found in the vanilla algorithm. \n\nI was also hoping some kind of baseline could be discussed, but I understand that to have a fair comparison they need to use the same pretrained GNNs with and without symmetry breaking.  \n\nAlso, the claims about generality of the approach and applicability to other graph problems, though reasonable, lack any substantive evidence or hint as to how to do it concretely. I suggest wither softening the claims or giving some hints as to why that is. \n\nIn terms of structure, I think the proof of theorem 1 is not necessary to be in the main text and occupies too much space. Instead some version of the last two tables in the appendix could be more useful in the main text. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for these helpful comments and suggestions.\n\n**Weakness [How many symmetric solutions are removed]**\\\nWe acknowledge that the limited information in Section 3.2 and Table 1 is insufficient for our purposes. Table 1 is not used to show that MIP can find *many* feasible solutions under different settings. These experiments numerically count *all* feasible solutions based on different levels of symmetry-breaking constraints.\n\nIn Table 1, the *(S1)* column is the baseline since each molecule should at least be a connected graph. Constraints S1 (which force connectivity) come from the literature and we wanted to show the improvement after adding our proposed symmetry-breaking constraints S2 and S3 with respect to S1. Therefore, the differences between the *(S1)* column and the *(S1)-(S3)* column are the number of removed symmetric solutions. Note that we did not use *symmetric molecules* as the reviewer mentioned. Please see our comments to your question about diversity of solutions.\n\nThe *(S1)* column is our baseline because, without adding S1, the baseline is the number of any graph with compatible features. If we ignore all features and only consider all possible graph structures, the baseline is $2^{\\frac{N(N-1)}{2}}$. This number will be much larger after introducing features, which loses meaning as a baseline.\n\nThe final paper will provide more information about the purpose of Section 3.2 and clarify the baseline in Table 1. We will add one more column for each dataset in Table 1 to show the percentage of removed symmetric solutions after adding S2 and S3. \n\n**Weakness [Generality consideration]**\\\nThanks for this great suggestion. Due to space limitation, we only briefly discussed the generality in Section 2.4 and Section 4. We provide some ideas here and will enrich the relevant discussion in the final paper to include these ideas.\n\n- The realization of these constraints is not limited to MIP: they also can exclude symmetric solutions in other optimization methods, such as genetic algorithms.\n\n- The constraints are independent of GNN: GNN in our work is just a permutation-invariant function. Our constraints directly work on graphs and features. In traditional computer-aided molecular design (CAMD), for example, the score functions are artificially-designed with closed forms. Since these functions are commonly-defined on features of molecules, they are also permutation-invariant. Similarly, other ML models could replace GNNs as long as they are permutation-invariant. Therefore, CAMD *without* a GNN is another possible application of our work.\n\nAs recommended by the reviewer, we will weaken the statements about the generality of the applications. The applications that we know of for this work are (i) GNN, (ii) CAMD, and (iii) the scheduling instances mentioned in Section 2.4. However, the statement of applications beyond these three is an aspirational one and we do not concretely know of other applications.\n\n**Weakness [Paper structure]**\\\nWe appreciate the concern about the paper structure. Since there will be one more page in final paper, we will try to add a simplified version of last two tables to the main text. We’re hoping that we can get both the proof and the tables into the final version (while also addressing the other promised changes).\n\n**Question [Diversity of solutions]**\\\nIn Table 1, fewer is always better, because it means that more symmetries are removed. Adding symmetry-breaking constraints will not influence the diversity of solutions. That is why we need to prove that there is at least one feasible indexing for each molecule. For example, if there are $n$ possible molecules in the design space, there are far more than $n$ feasible solutions due to symmetry, that is, multiple solutions correspond to the same molecule. The target of symmetry-breaking is to make the number of solutions as close as possible to $n$. No matter how many symmetries are removed, the diversity of solutions will not change because each molecule corresponds to at least one solution.\n\n**Minor questions**\n- line 75: We can just ignore the direction of each edge to transform a directed graph to an undirected graph.\n\n- line 106: Apologies for the ambiguity. $A$ is used in the definition of $LO(\\cdot)$ and the theoretical part to denote a multiset. After introducing constraints S1, there is no need of the adjacency matrix. Except for lines 106 - 126 and the proof of Property 3 (Section A.1), $A$ is always the adjacency matrix. In final paper, we will either clarify this point before line 106, or change the notation for multisets.\n\n- line 183: For a mixed-integer optimization formulation of GNNs, we do not require that the GNN is permutation invariant, but only need that the GNN could be rewritten in form (*). To optimize over a trained GNN with a non-fixed graph structure, we require permutation invariance. If we wanted to support GAT, then the symmetry-breaking constraints would be incorrect since different indexing would correspond different GAT outputs.\n\n- line 187: Yes, Dense NN is used to denote a neural network consist of fully connected layers (or to say, dense layers, linear layers). We will clarify this in the final paper."}, {"review_id": "Br025XdxCa", "reviewer": "Reviewer_5u3v", "summary": "This paper investigates the optimization of trained GNNs. This is a permutation-invariant problem and all points on an orbit of the symmetric group have the same performance. The authors proposed some symmetry-breaking approaches to have smaller search space and hence less redundancy and more efficient optimization. The correctness is proved theoretically and the efficiency is examined numerically.", "questions": "If you have Gurobi to detect and break the symmetry, will Gurobi add any constraints and if so, are those constraints the same as or different from (C26) and (C27)?", "limitations": "The authors clearly discuss the limitations in Section 4 and I agree with their discussion.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Optimizing a trained GNN is an important topic in practice and breaking symmetry can reduce the redundancy in many optimization problems.\n2. The authors provide the theoretical guarantee of the proposed algorithm.\n3. The numerical results look nice. It is reasonable to expect such improvement since the MIP formulation does have large symmetry group.", "weaknesses": "1. I do not think this paper has enough literature review on detecting and breaking MIP symmetry -- how do previous work detect symmetry? what is the cost for symmetry detection? what is the connection/difference between previous symmetry-breaking approaches and the proposed approach? I think Gurobi can detect symmetry but the computational cost might depend on the size of the symmetry group. For optimizing GNNs, the symmetry group is obvious (just the permutation group). So one does not need to detect the symmetry and can design problem-specific symmetry-breaking approaches. Those ideas/background should be discussed more in the paper.\n2. As the authors discuss in Section 4, the numerical experiments are not general enough, but I think the current experiments have already shown that the approach is promising.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for these considerations about the relationship and difference between our work and the literature.\n\n**Weakness [Literature review]**\\\nSincere thanks for this useful comment on the literature review: indeed we will add more details in the final paper. MIP solvers typically detect symmetry using graph automorphism, for example the open-source SCIP solver uses BLISS [i, ii], and then both Gurobi and SCIP break symmetry using orbital fixing [iii] and orbital pruning [iv].\n\nThe question about the cost of symmetry-breaking is great because it helps explain why our symmetry-breaking strategy is working so well. When a MIP solver detects symmetry, the only graph it has available is the graph formed by the variables and constraints in the MIP. So, for instance when $N = 4$, we only need to consider a permutation group with $N!=24$ elements. However, because the MIP solver does not have access to the graph structure, it needs to consider all possible automorphic graphs with $616+411$ nodes (these are the number of continuous and binary variables in the first line of Table 6 in our paper). So while Gurobi and the other MIP solvers *can* detect and break the symmetry, the computational cost of doing so is very high. We are reducing the computational cost for the solver by eliminating the need to look for these symmetries.\n\n`[i] Junttila, Tommi and Kaski, Petteri. Conflict propagation and component recursion for canonical labeling. TAPAS, 2011.`\\\n`[ii] Junttila, Tommi and Kaski, Petteri. Engineering an efficient canonical labeling tool for large and sparse graphs. ALENEX, 2007.`\\\n`[iii] Kaibel, Volker and Peinhardt, Matthias and Pfetsch, Marc E. Orbitopal fixing. Discrete Optimization, 2011.`\\\n`[iv] Margot, François. Pruning by isomorphism in branch-and-cut. Mathematical Programming, 2002.`\n\n**Weakness [Connection \\& Difference to the literature]**\\\nWe have (an admittedly short) answer to this question in lines 169 - 173 of Section 2.4. However, we will make this clearer by adding some details.\n\nThe most similar setting to ours is an application that distributes $m$ different jobs to $n$ identical machines and minimizes the total cost. The requirement is that each job can only be assigned to one machine (but each machine can be assigned to more than one job). Symmetry comes from all permutations of machines while variables $A_{i,j}$ denote if job $i$ is assigned to machine $j$. This setting appears in noise dosage problems [83, 84], packing and partitioning orbitopes [85, 86], and scheduling problems [87]. However, requiring that the sum of each row in $A_{i,j}$ equals to $1$ simplifies the problem. By forcing decreasing lexicographical orders for all columns, the symmetry issue is handled well. Our proposed constraints can be regarded as a non-trivial generalization of these constraints from a bipartite graph to an arbitrary undirected graph.\n\nThe main difference between the symmetry issue in this paper versus in the literature is that our variables are features and an adjacency matrix while the symmetry comes from alternative indexing of abstract nodes (e.g., atoms). As the reviewer mentioned, the symmetry group is the permutation group. Any element in the symmetry group acts on $N$ nodes to produce a isomorphic graph in graph level. In MIP, however, it corresponds to an element of symmetry group for all variables, including node features ($O(NF)$), adjacency matrix ($O(N^2)$), model (e.g., GNN), and problem-specific features (e.g., double bonds $DB$). By adding the proposed symmetry-breaking constraints, we don't need to consider a much larger symmetry group defined on all variables to find a subgroup which corresponds to the true symmetric group: the permutation group defined on abstract nodes."}, {"review_id": "XuDilO7CKy", "reviewer": "Reviewer_eLVB", "summary": "To overcome the symmetry issue when solving inverse problems on trained GNNs, this paper proposes two types of symmetry-breaking constraints to break symmetry. The authors construct an indexing algorithm, and prove that the resulting graph indexing satisfies the proposed symmetry-breaking constraints.  They also develop two mixed-integer optimization formulations in the case where the input\ngraph is not fixed. ", "questions": "1. I suggest that the authors provide complexity analysis to the proposed algorithm, since in Table 1 there are many time out cases. \n\n2. I suggest that the authors discuss the relationship between the proposed algorithms and related works (please see Weaknesses (2)).\n\n3. I suggest that the authors provide necessary experiments (please see Weaknesses (1)).", "limitations": "Yes.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The symmetry issue which this paper aims to overcome is important, and the proposed indexing is a general approach that applies in many kinds of situations. \n2. The theoretical analysis is solid, and the proof seems to be correct. ", "weaknesses": "1. The experiments are not persuasive enough. The reasons are as follows. \n\n(1) For Sec 3.1 (Mixed-integer optimization formulation for molecular design): In my opinion, a formulation itself WITHOUT experimental verification should not be thought of as a contribution. \n\n(2) For Sec 3.2 (Counting feasible structures: performance of symmetry breaking): The experimental results (Table 1) do not have baselines. \n\n2. There is no related work analysis (except a little description in Introduction). That makes it a little hard to justify the originality. For example, “The second formulation generalizes the big-M formulation for GraphSAGE in [70]” (L197). However, the authors did not explain how the big-M formulation looks like, in which way the proposed formulation generalizes it, and are there overlaps between the constraints of the big-M formulation and the proposed ones?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Sincerely appreciate these comments and questions.\n\n**Weakness 1(1) [Experimental verification of molecular constraints]**\\\nSection 3.1 is not a contribution of our paper but a description of how we followed the literature. This section just provides background for our numerical experiments. \n\nWe don't experimentally verify the Section 3.1 constraints because these are not new in the CAMD literature. On line 218 of Section 3.1, we explain the source of constraints C1 - C21 [70]. These constraints C1 - C21 are very similar to those which have been in the literature since 1993 [92].\n\nThe new constraints that we add are C22 - C25. The reviewer makes a good point that we should justify these constraints and provide an explanation where these constraints appear. These constraints bound the number of atoms, bonds, and rings. We calculated these bounds based on the QM7 and QM9 data sets, so every molecule in the QM7 and QM9 database will satisfy these bounds, respectively. Appendix C.1 gives more details as to how these bounds were calculated, so we will also provide a link to Appendix C.1 for C22 - C25. The reason that we do not experimentally verify these new molecular constraints is that removing the constraints would allow the GNN to extrapolate to molecules dissimilar to its dataset, for example removing the bound on the number of oxygen atoms could allow a molecule with 7 or 9 atoms to have 4 oxygen atoms when there is no example in the database that has 4 oxygen atoms.\n\n**Weakness 1(2) \\& Question 3 [Baseline of Table 1]**\\\nThanks for pointing this out. The *(S1)* column in Table 1 is the baseline, but we will clarify this in both Section 3.2 and the caption of Table 1. The reason that the *(S1)* column is the baseline is that these constraints come from the experimentally-verified literature and force the connectivity of molecules.\n\n**Weakness 2 \\& Question 2 [Originality of big-M formulation]**\n\nIt's a good point that we should more clearly explain how we generalize this big-M formulation. Currently (and we will make this more clear), the formulation in lines 197 - 204 is exactly the same as [70] except that we use our notation to keep consistency.\n\nThe way we generalize this formulation is not by changing these mathematical equations but rather by showing that this formulation is not limited to GraphSAGE. The mathematical formulation in lines 197 - 204 for all GNN architectures satisfying equation (*) in line 181 and the assumptions in line 191. We agree that the link from lines 181 \\& 191 to lines 197 - 204 are not clear, so we will make the connection explicit in Section 2.5.2.\n\n**Question 1 [Complexity analysis]**\\\nTable 1 counts all feasible solutions after adding more and more symmetry-breaking constraints. Since it's very hard to directly calculate these numbers analytically, we use MIP solver Gurobi to search the branch-and-bound tree to find every feasible solution. This process is super time-consuming. Our MIP for CAMD only involves binary variables, the maximum number of branches is $2^B$, where $B$ is the number of binary variables. That's why there are many time outs in Table 1: Gurobi can not visit all branches within the time limit (i.e., 48 hours)."}]

Author rebuttal:
To all reviewers,

Thank you very much for taking your precious time to review our paper and give great comments! Your feedback is really helpful for us to clarify and improve our work. More importantly, we are very lucky to have four reviewers focus on different aspects of our paper:

- Reviewer eveB identified the lack of an explanation of the compatibility between constraints (e.g., between S1 and S3, between C1 - C25 and S2 - S3).
    
- Reviewer ni19 commented on the improvement (e.g., how many symmetric solutions are removed, diversity of solutions) and generality of our proposed approach.
    
- Reviewer 5u3v asks about the relationship between our work and the symmetry-breaking literature (e.g., how does the literature detect and break symmetry, what are the connections and differences between our techniques and the literature). 
    
- Reviewer eLVB questioned  the experiments (e.g., experimental verification of molecular constraints, baseline, complexity analysis) and the originality of the big-M formulation.

As stated in each separate response, we will add the following clarifications to improve the paper based on the feedback of reviewers:

- A lemma after Theorem 1 will show the compatibility of S1 and S3 (as shown in the response to Reviewer eveB). A paragraph at the end of Section 3.1 will clarify that C1 - C25 will not affect the graph indexing.
    
- More information about the experiments in Section 3.2, including the purpose, the baseline, the percentage of removed symmetric solution, and the explanation of time out (as shown in the response to Reviewer ni19 and eLVB).
    
- A literature review about symmetry-breaking in Section 1, including the references we use in the response to Reviewer 5u3v. More discussion about the connections and differences between our approach and the literature in Section 2.4.
    
- Details about the generality of big-M formulation in Section 2.5.2, and the proposed symmetry-breaking constraints in Section 2.4 and Section 4 (as shown in the response to Reviewer ni19 and eLVB).

Thanks again for your comments, which are more carefully addressed in each separate response.

Meta-review: The paper addresses the issue of GNN optimization and demonstrates practical utility through molecular design application. The authors provide convincing experimental results supporting their claims, underscoring the proposed method's effectiveness. The manuscript is well-written, presenting complex ideas with notable clarity.
Decision: Accept (poster)