Equivariant Quantum Graph Neural Network for Mixed-Integer Linear Programming

Q1: What are the connections between this work and previous work on permutation-invariant quantum machine learning models (cited above)?

A4: Thanks for your question. If we understand correctly, "(cited above)" refers to arXiv:2211.16998 you mentioned. Our proposed model is quite different from the mentioned work. The work (arXiv:2211.16998) focuses on calculating ground states and time-evolved expectation values for permutation-invariant Hamiltonians. Some simple combinatorial optimization problems, such as the Maxcut problem, can correspond to a problem Hamiltonian, and we can find the solution by solving the ground state of the Hamiltonian. However, the way we solve the MILP problem is quite different from this class of problems. The MILP problem is abstracted as a bipartite graph, and predicting feasibility/optimal value/optimal solution is abstracted as a supervised graph-level/node-level classification/regression learning task. This can also be seen in the loss function used in section 3.5 of our latest submitted PDF.

Q2: What explains the terrible generalization performance of the classical GNNs? The authors should perform additional numerics fixing the parameter counts of the classical and quantum models to control for this behavior.

A5: Thanks for your question. As reported in Response to W1, (1) the number of parameters of the GNN is limited by its structure and is related to the embedding sizes. We cannot guarantee the parameters of GNN are exactly the same as our EQGNN. But we can compare their performance on a similar number of parameters. For example, we find two sets of results from the table above.

As we can see, even though GNN has a similar number of parameters to our EQGNN, it still doesn't perform very well. (2) The table of Response to W1 shows that GNN does not have good generalization error bounds in the MILP dataset, resulting in poor generalization performance. In fact, MILP is a very difficult problem, and Chen et al., 2023b are the first to apply GNN directly to approximate the solution of the MILP. The previous works only use GNN to help one part of their MILP algorithm, such as approximating Strong Branching or parameterizing cutting strategies. Moreover, Chen et al., 2023b only show how well GNN fits MILP training samples in their paper but do not show the performance of generalization on test sets. According to them, the generalization for MILPs on the foldable dataset is good, but the generalization for MILPs on the unfoldable dataset is not good enough. This is indeed the case in our tests. Figure 6 in our latest PDF (Figure 8 of the original PDF) shows the generalization error of GNN is small in the foldable dataset, but the generalization error on the unfoldable dataset is not good, as shown in the above tables. This may be caused by the difference between the MILP test set and the training set. In contrast, the good generalization performance shown by the proposed EQGNN is able to benefit the MILP community.

We hope that the above answers can address your concern. If you have any further questions, please contact me again. I am looking forward to your reply.

W2: It is also unclear to me whether there is actually any theoretical quantum advantage when the quantum model has no auxiliary qubits as arXiv:2211.16998 gives efficient classical simulation algorithms for permutation-invariant quantum systems and machine learning models. This might limit the utility of the introduced quantum model to the case where there are many auxiliary qubits, which then runs into problematic training behavior as the authors point out.

A2: Thanks for your comment. One point we would like to clarify is that the contribution of our model may not lie in the quantum advantage you mentioned, i.e., the proposed quantum model cannot be simulated by efficient simulation algorithms. In other words, even if our model can be classically simulated, it still has an advantage. As stated by Reviewer pyPd, our proposed EQGNN can "overcome a fundamental limitation of traditional GNNs (i.e., GNNs can not distinguish pairs of foldable MILP instances)". More specifically, according to statistics, around 1/4 of the problems in MIPLIB 2017 (Gleixner et al., 2021) [1] involved foldable MILPs. It means that practitioners using GNNs cannot benefit from that if there are foldable MILPs in their dataset of interest (Chen et al., 2023b) [2]. In contrast, our proposed EQGNN is able to distinguish graphs that GNN cannot distinguish while maintaining permutation equivariance, and thus can be applied to general MILP instances, which is important and meaningful for the MILP community.

Moreover, it may be good news if our model can be simulated by efficient classical algorithms. This is because classical GNNs have been adopted to represent mappings/strategies for MILP, for example, approximating Strong Branching (Gupta et al., 2020)[3], and parameterizing cutting strategies (Paulus et al., 2022)[4]. Based on the fundamental limitation of traditional GNNs, GNNs may struggle to predict a meaningful neural diving or branching strategy [1]. If our model can be simulated by efficient classical algorithms, we can replace the GNN part of these methods to better integrate with the classical framework. Thank you for letting us realize that our article is not clear enough in this part. We have added the above description to our revised paper.

References:

[1] Ambros Gleixner, Gregor Hendel, Gerald Gamrath, Tobias Achterberg, Michael Bastubbe, Timo Berthold, Philipp M. Christophel, Kati Jarck, Thorsten Koch, Jeff Linderoth, Marco Lubbecke, Hans D. Mittelmann, Derya Ozyurt, Ted K. Ralphs, Domenico Salvagnin, and Yuji Shinano. MIPLIB 2017: Data-Driven Compilation of the 6th Mixed-Integer Programming Library. Mathematical Programming Computation, 2021

[2] Ziang Chen, Jialin Liu, Xinshang Wang, and Wotao Yin. On representing mixed-integer linear programs by graph neural networks. In Proceedings of the International Conference on Learning Representations, 2023b

[3] Prateek Gupta, Maxime Gasse, Elias Khalil, Pawan Mudigonda, Andrea Lodi, and Yoshua Bengio. Hybrid models for learning to branch. Advances in neural information processing systems, 33: 18087–18097, 2020.

[4] Max B Paulus, Giulia Zarpellon, Andreas Krause, Laurent Charlin, and Chris Maddison. Learning to cut by looking ahead: Cutting plane selection via imitation learning. In International Conference on Machine Learning, pp. 17584–17600. PMLR, 2022.

W3: A much more minor point, but there are also many typos: "and yielding" at the bottom of page 3, "expressivity power" in the abstract, "guarantee to distinguish" in the abstract, and "TorchQauntum" at the bottom of page 7 are some that I found.

A3: Thanks for pointing out our grammatical/typographical mistakes.

• "and yielding" should be modified to "yield" to keep the form consistent with the previous verb "learn".
• "expressivity power" should be modified to "expressive power".
• "can guarantee to distinguish" can be modified to "can distinguish" because 'guarantee' is not appropriate here.
• "TorchQauntum" should be modified to "TorchQuantum". We have corrected these errors in the latest uploaded PDF, and further checked the full text for grammatical/typographical mistakes.

However, in order to further study the performance of GNN and further reduce the number of GNN parameters, we set the layer of GNN $L=1$, and the number of parameters is calculated by $16d^2+20d$. Note that the number of parameters when $L=1$ and $d=1$ is already the minimum number of parameters for this model. Then, we vary $d$ from 1 to 12 to get the following result:

As we can see, with the increase in embedding size, the training error of the model decreases continuously, but the generalization errors increase. Moreover, the test errors of GNN stay stable at around $0.09$, which indicates that GNN does not generalize well on the MILP dataset. As a comparison, we vary the number of blocks $T$ to test the performance of EQGNN, i.e.,

In contrast, the test and train errors of EQGNN can decrease as the parameter counts increase, and the generalization errors are always small, such that EQGNN can achieve better results than GNN on the unfoldable dataset. In addition, we can compare the GNN with $d=1$ with EQGNN with $T=3$, because they have the same number of parameters. In this case, EQGNN can achieve better results than GNN in both training and testing. Similarly, comparing GNN with $d=2$ with EQGNN with $T=9$, we can obtain similar results. Furthermore, as we can see, only when GNN has more than 696 parameters does its train performance exceed EQGNN with 120 parameters. However, GNN with 696 parameters has poor generalization performance, such that the test performance of GNN is much lower than that of EQGNN. We have added the above to Appendix H in our revised paper.

To further study the generalization performance of GNN, we test GNN on a larger dataset, i.e., the number of variables is 20, and the number of constraints is 6. The results of the GNN are:

From the result, we can see that the training error can decrease continuously, but the test error does not decrease starting from embedding size 3. In other words, when embedding size=3, it is the best GNN model on the test set, but the number of GNN parameters at this time is not enough to fit the training set well. As the number of parameters increases, GNN gets better at fitting the training set, but GNN does not generalize to the test set. This experiment shows the problem of GNN generalization on the MILP dataset. The above performance results of the GNN are obtained by using the code provided by the original paper (Chen et al., 2023b)(https://github.com/liujl11git/GNN-MILP).

W1: I am more skeptical of the numerical results demonstrating a separation in expressive power between the EQGNN and the GNN the authors cite from (Chen et al., 2023b) (though admittedly I am not an expert on GNNs). The divergence of testing performance in Fig. 10 seems to me that the classical GNN is overfitting, potentially due to the order of magnitude difference in the numbers of parameters between the quantum and classical models. I highly recommend the authors perform supplemental numerics where these parameter counts are brought in line to control for this behavior.

A1: Thanks for your suggestion. Next, we add more numerical experiments to explain this problem. But first, we need to clarify that the number of parameters of the GNN is limited by its structure. That is, we cannot guarantee that the number of parameters in the GNN is exactly linear. More specifically, the GNNs used in (Chen et al., 2023b) are constructed by the following structure.

(1) The initial embedded features are $s_i^0 = f_{in}^V(h_i^V)$ and $t_j^0 = f_{in}^W(h_j^W)$, where $V$ represents the the set of variable nodes, and $W$ represents the the set of constraint nodes. $h_i^V$ indicates the input features of variable nodes, and $h_j^W$ indicates the input features of constraint nodes. $f_{in}^V: \mathbb R^{d_{in}} \rightarrow \mathbb R^{d}$, where $d_{input}$ is the dimension of input features and $d$ is defined as the embedding size. $f_{in}^V$ and $f_{in}^W$ are the learnable functions, which are parameterized with multilayer perceptrons (MLPs), and have two hidden layers with a ReLU activation function. The number of parameters of $f_{in}^V$ is calcuated by $2*(d_{in}*d+d)$.

(2) Message-passing layers for $l = 1, 2, . . . , L$. Given learnable mappings $\\{f_l^V , f_l^W , f_l^V , f_l^W, f_{l_{out}}^V , f_{l_{out}}^W \\}_{l=0}^L$, one can update vertex features by the following formulas:

$s_i^l = f_{l_{out}}^V(f_l^V(s_i^{l-1}), \sum_{j=1}^n E_{i,j} f_l^W(t_j^{l-1})), \quad t_j^l = f_{l_{out}}^W(f_l^W(t_j^{l-1}),\sum_{i=1}^m E_{i,j} f_l^V(s_i^{l-1})),$ where all learnable mappings $f$ are also parameterized with multilayer perceptrons (MLPs) and have two hidden layers with a ReLU activation function.

(3) For predicting feasibility, the output is $f_{out} (\sum_{i=0}^n{s_i^L}, \sum_{j=0}^n{t_i^L})$, $f_{out}: \mathbb R^{d} \rightarrow \mathbb R$, and the number of parameters is $2*d^2 + 2d$.

Therefore, the number of parameters of the GNN is $(14L+2)d^2 + (12L+8)d$ in the case where $h_i^V \in \mathbb R^4$ and $h_j^W \in \mathbb R^2$. In the paper (Chen et al., 2023b), they set $L = 2$ for the GNN. Note that due to the over-smoothing of GNNs, most GNNs only stack two to four layers, i.e., $L = 2 \sim 4.$ Therefore, the best way to control performance is to vary the embedding size $d$. When $L=2$, the number of parameters of the GNN is $30d^2+30d$. Since $d \in \mathbb Z^+$, the number of parameters of GNN can only take $60, 180, 360, 600, 900,...$.

Our model uses the number of blocks $T$ to control the number of parameters, and each block has 12 parameters, so the number of our parameters is $12*T$. We converted Figure 7(a) of our paper into a table form and attached the embedding size of the GNN and the number of blocks of our EQGNN, as shown below.

The results of the GNN:

The results of our EQGNN:

As we can see, the purpose of this experiment is to control the embedding size of the GNN and the number of blocks of EQGNN to be consistent. In all our experiments, we first gradually increase the embedding size/blocks to test the performance of the models and then find the embedding size/blocks $d^*$ or $T^*$ corresponding to the best performance. Then, we select the values near $d^*$ or $T^*$ and show their corresponding results in the experiments of the paper.

The MILP graph is represented by $G = (V \cup S, A)$ and $\mathcal{G}_{q,p}$ as the collection of all such weighted bipartite graphs whose two vertex groups (variables and constraints) have size $q$ and $p$, respectively. All the vertex features are stacked together as $H = (\boldsymbol{f}^V_1,...,\boldsymbol{f}^V_q,\boldsymbol{f}^S_1,...,\boldsymbol{f}^S_p) \in \mathcal{H}^V_q \times \mathcal{H}^S_p$.

$hash(·)$ is a function that maps its input feature to a color in $\mathcal C$. The algorithm flow can be seen as follows. First, all nodes in $V \cup S$ are assigned to an initial color $c_v^0$ and $c_s^0$ according to their node features. Then, for each $v_i \in V$, $hash$ function maps the previous color $c_{v_i}^{l-1}$ and aggregates the color of the neighbors of $\\{c_{s_j}^{l-1}\\}_{{s_j}\in \mathcal N(v_i)}$.

Similarly, $hash$ function maps the previous color $c_{s_j}^{l-1}$ and aggregates the color of the neighbors of $\\{c_{v_i}^{l-1}\\}_{{v_i}\in \mathcal N(s_j)}$. This process is repeated until $L$ reaches the maximum iteration number. Finally, a histogram $h_G$ of the node colors can be obtained according to $\\{\\{c_v^l: v\in V, c_s^l: s\in S \\}\\}$, which can be used as a canonical graph representation. The notation $\\{\\{·\\}\\}$ denotes a multiset, which is a generalization of the concept of a set in which elements may appear multiple times: an unordered sequence of elements. That is, the 1-WL test transforms a graph into a canonical representation, and then if the canonical representation of two graphs is equal, the 1-WL test will consider them isomorphic.

Next, we show a case where the 1-WL test fails. Taking Figure 2 in our paper as an example, we feed (a) and (b) from Figure 2 into this algorithm separately and compare the differences at each step. In the setting of Figure 2, $\{\boldsymbol{f}^V_i\}$ are equal, $\{\boldsymbol{f}^S_j\}$ are equal, are equal, and the weights of all edges are equal. In these two graphs, each node has two neighbors, but the only difference lies in their topological structure. More importantly, these two graphs have different feasibility, i.e., (a) is feasible while (b) is infeasible. Initially, each variable vertex has the same color $c_v^0$, and each constraint vertex has the same color $c_s^0$, because their features are the same. This step is the same for Figure 2(a) and Figure 2 (b). Then, each vertex iteratively updates its color based on its own color and information from its neighbors. Note that, since each $v_i$ is connected to two $s_{j_1}$, $s_{j_2}$ in both Figure 2(a) and 2(b), they will also obtain the consistent $c_v^l$ and $c_s^l$ at this step. That is, Figure 2(a) and 2(b) will get the same representation by the algorithm, and they cannot be distinguished by 1-WL test GNNs. Moreover, Xu et al. [7] show that GNNs are at most as powerful as the 1-WL test in distinguishing graph structures. Therefore, a pair of 1-WL indistinguishable graphs are indistinguishable by GNNs.

For EQGNN, the subfigure (a) and (b) of Figure 2 can be clearly distinguished. Since the two graphs differ only in the connectivity of edges, we only consider the modules associated with edges in EQGNN, i.e., $U_{g1}$ and $U_{g2}$. Specifically, $U_{g_1}(G,\beta) = \text{exp}(-**i**(\sum_{(i,j)\in \mathcal E} (A_{i,j}+\beta) Z_{2i}Z_{2q+j}))$, where $i\in\\{0,1,...,q\\}$ indicates the index of variable nodes and $j\in\\{0,1,...,p\\}$ indicates the index of constraint nodes. Since two qubits represent one variable node, $2q$ qubits are used to represent all variables, and the index of the qubit representing the constraint node is numbered from $2q$. In subfigure (a) node $v_3$ is connected to $s_3$ and $s_4$, while in subfigure (b) node $v_3$ is connected to $s_3$ and $s_1$. Therefore, the difference between $U_{g_1}(G_a,\beta)$ and $U_{g_1}(G_b,\beta)$ lies in the two-qubit gate acted on between qubits representing $v_3$ and $s_4$ and the two-qubit gate acted on between qubits representing $v_3$ and $s_1$. After the nodes correspond to their qubit index, we can obtain the formulas $\text{exp}(-**i**((A_{3,4}+\beta) Z_{6}Z_{16}))$ and $\text{exp}(-**i**((A_{3,1}+\beta) Z_{6}Z_{13}))$. Since the edge weights and the parameter $\beta$ are equal, we can just compare $\text{exp}(-**i**(Z_{6}Z_{13}))$ and $\text{exp}(-**i**(Z_{6}Z_{16}))$. These two values are obviously not equal, so $U_{g_1}(G_a,\beta) \neq U_{g_1}(G_b,\beta)$ and the overall unitary of our model $U(G_a,H) \neq U(G_b,H)$. Thereby, EQGNN can distinguish such graphs that the 1-WL test and GNNs cannot distinguish. We have added the above to the revised paper, see Appendix B. The experiments in Figure 6 show that the performance of EQGNN can achieve an error rate close to zero in these foldable datasets. In contrast, the error rate of GNNs stabilizes at 0.5.

(To be Continued)

Thanks for your insightful comments. We appreciate that some strengths of our paper are recognized, such as fairly clear presentation, and the interesting and unique way to solve MILPs. However, there seems to be some misunderstanding concerning our article, and we will do our best to provide a detailed explanation to address your concerns and improve our manuscript according to your suggestions. This is followed by point-to-point responses.

Question 1: It is not clear to me why there are instances of MILP that cannot be distinguished by GNNs, as in Figure 2. While the vertex degrees are the same, the connectivity is different, so shouldn't GNNs treat them differently? Perhaps I am missing something about how standard GNNs deal with this problem in the context of MILP.

A1: This inability to distinguish is related to the essential limitation of the neighbor-aggregation mechanism of GNNs, and this limitation does not only exist in the context of MILP, which we will elaborate on next.

To better understand this limitation, let's first take a more fundamental example, which is also not distinguishable by GNNs. The problem of distinguishing two graphs can be abstracted as a graph isomorphism testing problem, i.e., to determine whether two graphs are isomorphic.

Definition 1: Two graphs G1 and G2 are isomorphic if there exists a matching between their vertices so that two vertices are connected by an edge in G1 if and only if corresponding vertices are connected by an edge in G2.

That is, two graphs that are isomorphic means that they will have the same topological structure or connectivity. Let's look at a simple example of non-isomorphic graphs, which you can see by clicking on this Figure R1. Both graphs are 2-regular graphs, and the regular graph is a graph where each vertex has the same number of neighbors, i.e., every vertex has the same degree. The only difference is that Figure R1 (a) has two cycles, and Figure R1 (b) has only one cycle.

The 1-dimensional Weisfeiler-Lehman (1-WL) algorithm is a well-known heuristic for deciding whether two graphs are isomorphic. It is known to work well in general, with a few exceptions, e.g., regular graphs [1,2]. Moreover, Xu et al. [3] show that GNNs are at most as powerful as the 1-WL test in distinguishing graph structures. Therefore, a pair of 1-WL indistinguishable graphs are indistinguishable by GNNs. Next, we will briefly introduce the 1-WL test algorithm and why the Figure R1 cannot be distinguished by 1-WL or by GNNs.

Algorithm 1: 1-WL (color refinement)

Input: $G = (V, E, X_V)$

1. $c_v^0 = hash(X_v)$ for all $v\in V$
2. repeat
3. $\quad c_v^l = hash(c_v^{l-1}, \\{\\{c_w^{l-1}: w \in \mathcal N_G(v)\\}\\})$, for all $v \in V$
4. until $c_v^l= c_v^{l-1}$ for all ${v\\in V}$
5. return $\\{\\{c_v^l: v\in V \\}\\}$

In the above algorithm, we consider an undirected graph $G = (V,E,X_V)$, where $V$ is the set of vertices; $E$ is the set of edges; and $X_V$ is the set of node features. The notation $\\{\\{·\\}\\}$ denotes a multiset, which is a generalization of the concept of a set in which elements may appear multiple times: an unordered sequence of elements. $hash(·)$ is a function that maps its input to a color in $\mathcal C$. The algorithm flow can be seen as follows. First, all nodes in $G$ are assigned to an initial color $c_v^0$ according to their node features. Then, $hash$ function maps $v$’s previous color $c_v^{l-1}$ and the multiset of the previous colors of $v$’s neighbors $\\{\\{c_w^{l-1}: w \in \mathcal N_G(v)\\}\\}$ to a color $c_v^l$ in $\mathcal C$, for each $v \in V$. This process is repeated until the colors of all nodes do not change, i.e., $c_v^l= c_v^{l-1}$ for all ${v\\in V}$. Finally, a histogram $h_G$ of the node colors can be obtained according to $\\{\\{c_v^l: v\in V \\}\\}$, which can be used as a canonical graph representation. That is, the 1-WL test transforms a graph into a canonical representation, and then if the canonical representation of two graphs is equal, 1-WL will consider them isomorphic.

In Figure R1, we assume that all node features are equal, and thereby all nodes are assigned to the same initial color. Then, the hash function aggregates the color of $v$ and the colors of the two neighbors of $v$. It should be noted that Figure R1 (a) and Figure R1 (b) have the same result in this step because their vertex degrees and node features are the same. This also directly leads to the fact that Figure R1 (a) and Figure R1 (b) will certainly have the same canonical graph representation after the 1-WL algorithm. However, they are clearly non-isomorphic graphs, and 1-WL fails to distinguish them.

(to be continued)

A7: Thanks for your comment. The $R_{zz}$ gates with different control qubits indeed have the same unitary matrix. However, Figure R3 shows an example of using different control qubits. As we can see, although the input quantum state and output quantum state of the whole circuit are kept consistent for Figure (a) and (b), the position where the global phase appears on a single qubit will be different, such as Figure (a) outputs $|\psi_0^a\rangle = |1\rangle$ and Figure (b) outputs $|\psi_0^b\rangle = e^{i\frac{\theta}{2}}|1\rangle$. However, the global phase $e^{i\frac{\theta}{2}}$ has no effect on the observations resulting from a measurement of the state, and it is only the artefact of the mathematical framework and has no physical meaning. We specify the control qubits solely for the certainty of the mathematical formulation. Nevertheless, our specification of control qubits and the use of cross symbols may cause misunderstanding, so we modified Figure 6 in our paper, please see Figure 5 of the latest submitted PDF. Thank you for your suggestion.

Comment 3: Typo in the sentence "the identical parametric gates are acted when the order of input nodes."

A8: Thanks for pointing out the problem. The sentence is incomplete and should be modified to "the identical parametric gates are acted on the circuit no matter how the order of input nodes changes."

Comment 4: Typo in the sentence "We now study the effect of width of the circuit increased"

A8: We apologize for not writing this sentence clearly. It should be modified to "We now study the influence on the performance with the increase of circuit width".

Weakness 2: Unless I am misunderstanding something, I feel that this permutation equivariance is not particularly insightful. For instance, the equivalent circuits in Figure 7 seem obvious, just that the circuit wires are drawn to either have different input order or output order. Perhaps it would be more useful to show a construction that fails to satisfy permutation equivariance.

A9: First of all, we need to clarify that Figure 7 is just a very simplified example of showing permutation equivariance, and it does not represent our circuit, and designing permutation-equivariant circuits is not easy. To address your concerns, we have rewritten Section 3.5 and 3.6 in the latest uploaded PDF. Section 3.5 formulates the proposed model. In particular, Section 3.6 presents four definitions and three theorems to prove the permutation equivariance of our proposed model. Since this part is about two pages long, please refer to the highlighted section of our latest uploaded PDF. In addition, we also provide a shorter version to answer the question.

One of the most intuitive examples of failure to satisfy permutation equivariance is HEA (Hardware efficient ansatz). Figure R4 shows an example of an HEA with three qubits, which is used to learn the representation of graphs with three nodes. Figure R4 uses $R_x$ gate to encode the feature of three nodes, and uses $RY(\theta_0), RY(\theta_1)$, $RY(\theta_2)$ and CNOT gate as the learning module. Suppose that the order of node is $(v_0, v_1, v_2)$, we measure each qubit and get the probability of $|0\rangle$ at each qubit as the output $(y_{v_0}, y_{v_1}, y_{v_2})$. Then, the order of node changes to $(v_1, v_0, v_2)$ and the output is $(y_{v_1}^{\\prime}, y_{v_0}^{\\prime}, y_{v_2}^{\\prime})$. Obviously, $y_{v_0} \neq y'_{v_0}$. The permutation equivariance of the model is vital for graph learning, and the output of the model should be reordered consistently with the permutation on the input nodes, otherwise the model may overfit to the order of nodes in the training data. In the field of classic machine learning, GNN can use the neighborhood-aggregation mechanism to achieve the permutation equivariance. However, in the field of quantum computing, the permutation-equivariant model requires fine designs due to the constraints of quantum systems.

Weakness 1: I feel that many parts of the construction of the ansatz is not well motivated. However, this seems to be a typical issue for quantum neural networks, perhaps even moreso than classical neural networks.

A10: Thanks for your comment. In fact, because our network is required to satisfy permutation equivariance, each layer of our model must be designed for encoding or learning while ensuring permutation equivariance, which is much more motivated than other QNNs. We believe we have addressed some of your concerns in our previous answers A2, A3, and A9. Moreover, the new content added in section 3.6 in our current submitted PDF can further address your concerns.

We hope that the above answers can address your concern. If you have any further questions, please contact me again. I am looking forward to your reply.

• The constraints for the $2k$-th problem (infeasible) is $x_{j_1}+x_{j_2}=1$, $x_{j_2}+x_{j_3}=1$, $x_{j_3}+x_{j_1}=1$, $x_{j_4} + x_{j_5} = 1$, $x_{j_5} + x_{j_6} = 1$, $x_{j_6} + x_{j_4} = 1$. An explanation for why it is infeasible: we add the first three formulas together to get $2(x_{j_1}+x_{j_2}+x_{j_3})=3$, but x_j \in \\{0, 1\\\} for $j \in J$, so there will be no solution that satisfies this equation and the MILP problem must be infeasible.

For unfoldable MILPs, we first set the number of variables and constraints to $m$ and $n$.

• For each variable, $c_j ∼ \mathcal N(0, 0.01)$, $l_j, u_j ∼ \mathcal N(0, \pi)$. If $l_j > u_j$, then switch $l_j$ and $u_j$. The probability that $x_j$ is an integer variable is 0.5.
• For each constraint, $\circ_i ∼ \mathcal U({≤, =, ≥})$ and $b_i ∼ \mathcal U(-1, 1)$.
• After randomly generating all the MILP samples, we use the 1-WL test algorithm to calculate their graph representation for each instance, ensuring that there are no duplicate graph representations in the dataset so that we can determine that this dataset does not contain 1-WL indistinguishable pairs of MILP instances.

We have added the above description of dataset generation to our revised paper, and please see Appendix D in the latest PDF submission. Thanks for your suggestion.

References:

[1] Ziang Chen, Jialin Liu, Xinshang Wang, and Wotao Yin. On representing mixed-integer linear programs by graph neural networks. In Proceedings of the International Conference on Learning Representations, 2023b.

Question 5: When predicting the solution vector, what is the actual representation of the classical information after reading out the qubits? For example, does a 2-qubit state for $v_1$ correspond to integers 0,1,2,3? Also, what happens when reading out a result that does not satisfy the constraints?

A5: Thanks for your question. As mentioned earlier, the hyperparameter $\omega$ can control the number of qubits representing each node. If the number of qubits is more than 1, like the node $v_1$ in Figure 6, we will add the control gate at the end of the node update layer, such as the CRY gate in Figure 6. Then, at the end of the circuit, we only measure the first qubit representing the node. As mentioned in section 3.5 of our paper, we perform the Pauli-Z measurement on the qubit, which can obtain the $P(|0\rangle)-P(|1\rangle)$, which ranges from -1 to 1. That is, each variable node can obtain one value $F_{sol}^i(G)$ to predict the values of variables. In addition, the maximum range $\lambda$ of variables of the training sample is known. The solution vector of the final output of the model is $\\{\lambda F_{sol}^i(G)\\}_{i=0}^{q-1}$.

Assuming the groundtruth of the solution vector is $\\{y_{sol}^i(G)\\}_{i=0}^{q-1}$, we evaluate the performance of the predicted solution vector by:

$\frac{1}{Mq}\sum_{m=0}^{M-1} \sum_{i=0}^{q-1}(y_{sol}^i(G_m)-\lambda F_{sol}^i(G_m))^2$

where $G_m$ indicates the $m$-th tested MILP graph, and the number of tested instances is $M$. As we can see, the values of the solution vector are continuous. For integer variables, a post-processing can be applied to make a rounding approximation. In addition, we need to make it clear that the proposed method is similar to the previous methods [1], only providing an approximate solution rather than providing an exact solution. We have added the above explanation to our latest uploaded PDF, and please see our highlights in Section 3.5.

Comment 1: In the definition of the feasible region $X_{fea}$ in Section 2, the constraint $l \leq x \leq u$ is missing.

A6: Thank you for pointing this out, and we have corrected the definition of the feasible region $X_{fea}$ in our latest uploaded PDF.

Comment 2: In Figure 6, I feel the circuit for $R_{zz}$ gates is somewhat misleading. The cross symbol is typically used for the SWAP gate. Also, the $R_{zz}$ is symmetric with respect to the two qubits it acts on, so there is no difference between choosing which is the target and which is the control qubit.

(to be continued)

For a more intuitive understanding, let us give a simple example of why such a property is needed in our circuit and why $R_{zz}$ are chosen over other two-qubit gates. As shown in Figure R2 (this is a hyperconnection), suppose there are three qubits representing three nodes, and there is an edge $e_{12}$ between $v_1$ and $v_2$, and there is an edge $e_{23}$ between $v_2$ and $v_3$. Each edge is associated with a two-qubit gate. For learning graphs, the model needs to ensure that the circuit is the same unitary regardless of the arrangement of the edges. For example, the order in which the two-qubit gate associated with the edges acts should not affect the model results. $R_{zz}$ gate are carefully selected to ensure that the edge arrangement invariance, while $R_{xx}$ and $R_{yy}$ do not have this property. For example, $R_{zz}$ is a diagonal matrix, then $I \otimes R_{zz}(\theta_{23})$ or $R_{zz}(\theta_{12}) \otimes I$ are still diagonal matrix, and that diagonal matrices commute, i.e., $[R_{zz}(\theta_{12})\otimes I][I\otimes R_{zz}(\theta_{23})] = [I\otimes R_{zz}(\theta_{23})][R_{zz}(\theta_{12})\otimes I]$. Hence, we can arbitrarily choose the order of action of the two-qubit gates associated with the edge. However, $[R_{xx}(\theta_{12})\otimes I][I\otimes R_{xx}(\theta_{23})] \neq [I\otimes R_{xx}(\theta_{23})][R_{xx}(\theta_{12})\otimes I]$ and $[R_{yy}(\theta_{12})\otimes I][I\otimes R_{yy}(\theta_{23})] \neq [I\otimes R_{yy}(\theta_{23})][R_{yy}(\theta_{12})\otimes I]$, which they do not meet the design requirements of QGNNs. Thank you for letting us realize that our article is not clear enough in this part. We have added the above description to our revised paper, please refer to our latest PDF submission, especially Theorem 3.

Question 3: In the feature encoding layer, is there any reason why $R_x(c_i)$ and $R_x(u_i)$ as well as $R_x(l_i)$ and $R_x(\epsilon_i)$, are applied to different qubits?

A3: As mentioned in Section 3.2, we set a hyperparameter $\omega$ to control the number of features on each qubit in our method. Figure 5 of our paper only shows an example of $w = 2$. Each variable has four features, and two qubits represent it. Of course, we can also set $\omega$ to 4 so that $(c_i, u_i, l_i, \epsilon_i)$ is encoded on the same qubit. As we can see, the smaller the value of $w$, the larger the number of qubits representing one node. This means that the higher the dimension of the feature representation we can learn for each node, the more expressive the model will be. The choice of $\omega$ is a trade-off between the available quantum resources (the number of qubits) and the expressive power of the model. In addition, since each feature dimension is associated with a learnable parameter, the order in which these four features are used is not important.

Question 4: Out of curiosity, how the MILP instances are generated to be foldable or unfoldable? I suppose this requires something like solving the graph isomorphism problem.

A4: Thanks for your questions. The generation of foldable or unfoldable MILP instances is consistent with that of the reference (Chen et al., ICLR 2023b) [1], which does not involve solving the graph isomorphism problem. Next, we will introduce the generation process in detail and add this part to our revised paper.

Note that the foldable and unfoldable datasets are generated separately. Foldable dataset is constructed by many pairs of 1-WL indistinguishable graphs, and figure 2 in our paper is a foldable example, which is a pair of non-isomorphic graphs that cannot be distinguished by WL-test or by GNNs, as mentioned in the response to Question 1. In contrast, the unfoldable dataset refers to other MILP instances that do not have 1-WL indistinguishable graphs.

In section 4.1, we randomly generate 2000 foldable MILPs with 12 variables and $6$ constraints, and there are 1000 feasible MILPs with attachable optimal solutions while the others are infeasible. We construct the $(2k − 1)$-th and $2k$-th problems via following approach $(1 ≤ k ≤ 500)$.

• Sample $J = \\{j_1, j_2, . . . , j_6\\}$ as a random subset of $\\{1, 2, . . . , 12\\}$ with 6 elements.

1. For $j \in J$, $x_j \in \\{0, 1\\}$, i.e., $x_j$ is a binary integer variable.
2. For $j \notin J$, $x_j$ is a continuous variable with bounds $l_j∼ \mathcal{U}(0, \pi), u_j∼ \mathcal{U}(0, \pi)$. If $l_j > u_j$, then switch $l_j$ and $u_j$.

• $c_1 =. . . = c_{12} = 0.01$.
• The constraints for the $(2k−1)$-th problem (feasible) is $x_{j_1}+x_{j_2}=1$, $x_{j_2}+x_{j_3}=1$, $x_{j_3}+x_{j_4}=1$, $x_{j_4} + x_{j_5} = 1$, $x_{j_5} + x_{j_6} = 1$, $x_{j_6} + x_{j_1} = 1$. For example, $x = (0, 1, 0, 1, 0, 1)$ is a feasible solution.

(to be continued)

As mentioned earlier, one well-known limitation of GNNs is their separation power, which is upper bounded by the 1-WL algorithm for graph isomorphism testing [3] because they have similar neighbor aggregation mechanisms. Therefore, for a pair of 1-WL indistinguishable graphs, as shown in Figure R1 or Figure 2, any GNNs will learn the exact same representations for these graphs, yielding the same prediction for both.

As we can see, Figure 2 differs from Figure R1 in that it is a bipartite graph for the MILP task. To prove that Figure 2 is also a pair of 1-WL indistinguishable graphs, we present a variant of 1-WL test specially modified for MILP that follows the same lines as in (Chen et al. 2023b, Algorithm 1)[4].

Algorithm 2: WL test for MILP-Graphs

Input: A graph instance $(G, H) \in \mathcal{G}_{q,p} \times \mathcal{H}^V_q \times \mathcal{H}^S_p$, and iteration limit $L > 0$.

1. Initialize with $c_v^0 = hash(\boldsymbol f^V)$ for all $v\in V$, $c_s^0 = hash(\boldsymbol f^S)$ for all $s\in S$.
2. for $l = 1,2,...L$ do
3. $\quad c_{v_i}^l = hash(c_{v_i}^{l-1}, \sum_{j=0}^{p-1} A_{i,j} hash(c_{s_j}^{l-1}))$, for all $v \in V$
4. $\quad c_{s_j}^l = hash(c_{s_j}^{l-1}, \sum_{i=0}^{q-1} A_{i,j} hash(c_{v_i}^{l-1}))$, for all $s \in S$
5. end for
6. return the multisets containing all colors $\\{\\{c_v^l: v\in V, c_s^l: s\in S \\}\\}$.

We can represent the MILP graph as $G = (V \cup S, A)$ and $\mathcal{G}_{q,p}$ as the collection of all such weighted bipartite graphs whose two vertex groups have size $q$ and $p$, respectively. All the vertex features are stacked together as $H = (\boldsymbol{f}^V_1,...,\boldsymbol{f}^V_q,\boldsymbol{f}^S_1,...,\boldsymbol{f}^S_p) \in \mathcal{H}^V_q \times \mathcal{H}^S_p$. Thereby, the weighted bipartite graph with vertex features

$(G,H) \in \mathcal{G}_{q,p} \times \mathcal{H}^V_q \times \mathcal{H}^S_p$

contains all information on the MILP problem. Next, we feed (a) and (b) from Figure 2 into this algorithm separately and compare the differences at each step. In the setting of Figure 2, $\{\boldsymbol{f}^V_i\}$ are equal, $\{\boldsymbol{f}^S_j\}$ are equal, and the weights of all edges are equal. Therefore, initially, each variable vertex has the same color $c_v^0$, and each constraint vertex has the same color $c_s^0$. This step is the same for Figure 2(a) and Figure 2 (b). Then, each vertex iteratively updates its color, based on its own color and information from its neighbors. Note that, since each $v_i$ is connected to two $s_{j_1}$, $s_{j_2}$ in both Figure 2(a) and 2(a), they will also obtain the consistent $c_v^l$ and $c_s^l$ at this step. That is, Figure 2(a) and 2(b) will get the same representation by Algorithm 2, and they cannot be distinguished by 1-WL test and GNNs. We have added the above to our latest uploaded PDF, please see Appendix B.

References:

[1] Jin-Yi Cai, Martin Fürer, and Neil Immerman. An optimal lower bound on the number of variables for graph identification. Combinatorica, 12(4):389–410, 1992.

[2] Brendan L Douglas. The weisfeiler-lehman method and graph isomorphism testing. arXiv preprint arXiv:1101.5211, 2011.

[3] Xu, K., Hu, W., Leskovec, J., and Jegelka, S. How powerful are graph neural networks? In ICLR, 2019.

[4] Ziang Chen, Jialin Liu, Xinshang Wang, and Wotao Yin. On representing mixed-integer linear programs by graph neural networks. In Proceedings of the International Conference on Learning Representations, 2023b.

Question 2: Why choose $R_{zz}$ over other two-qubit gates that commute with itself such as $R_{xx}$ or $R_{yy}$? This goes back to the weakness of the circuit construction not being well motivated.

A2: First of all, we need to clarify that the commute in this article does not refer to $R_{zz}(\theta_{12})R_{zz}(\theta_{23}) = R_{zz}(\theta_{23})R_{zz}(\theta_{12})$, but to $[R_{zz}(\theta_{12})\otimes I][I\otimes R_{zz}(\theta_{23})] = [I\otimes R_{zz}(\theta_{23})][R_{zz}(\theta_{12})\otimes I]$. More specifically, it refers to $\text{exp}(-\text{i}(\theta_{i_1,j_1} Z_{i_1} Z_{j_1} + \theta_{i_2,j_2} Z_{i_2} Z_{j_2})) = \text{exp}(-\text{i}(\theta_{i_2,j_2} Z_{i_2} Z_{j_2} + \theta_{i_1,j_1} Z_{i_1} Z_{j_1}))$, where $Z_{i_1} = I \otimes ... Z_{i_1} \otimes ... \otimes I$. Since $\text{exp}(-\text{i} A_{i,j} Z_i Z_j)$ is a diagonal matrix, and all diagonal matrices are commute, the equation is proven. However, other two-qubit gates, e.g., $\text{exp}(-**i** A_{i,j} X_i X_j)$ does not satisfy the equation. In other words, the $\text{R}_\text{zz}(\theta) = \text{exp}(-**i** \theta Z_i Z_j)$ gate is carefully chosen. In addition, we have added the above to section 3.6 of the updated PDF and shown why this equation is required to hold.

Dear Reviewer AMEi,

Thank you very much for your reply. We would like to express our sincere gratitude for your valuable feedback. You are right. It is true that $X \otimes X \otimes I$ commute $I \otimes X \otimes X$. We apologize that the previous statement in the rebuttal is inappropriate. $R_{XX}$ gate and $R_{YY}$ gate are similar to $R_{ZZ}$ gate, and they can be used to preserve the permutation invariant of edges. $R_{ZZ}$ gate is not the only option to preserve the permutation invariant of edges, and the two-qubit gates that can commute in the circuit are able to be used to learn the graph information interaction. We have added the above to our latest uploaded PDF, please see the second to last paragraph of section 3.6.

In addition, we want to point out that our proposed model provides a general framework for designing equivariant QNN on learning graph structure data besides the specific construction of the circuits. Specifically, as mentioned in section 3.6, an equivariant GQNN can be built by the layers independent of input orders and the layers related to input orders. The former can be constructed by the proposed parameter-sharing mechanism and the latter can be established by choosing the two-qubit gates that can commute in the circuit to preserve the permutation of edges. Moreover, the proposed multi-qubit encoding mechanism and optional auxiliary layer can also provide new ideas for the construction of other EQNNs to handle high-dimensional data features. In other words, we provide a design rule to build an equivariant graph QNN.

We are grateful for your feedback that helped improve our submission.

Best regards

Dear Reviewer AMEi,

Thanks again for your review. We have replied to the problem about the choice of two-qubit gates in the last response (we have further updated it just now). We wanted to reach out to see if our response has addressed your concerns.

In addition, we would appreciate it if you could consider updating our score.

We are more than happy to discuss further if you have any further concerns and issues, please kindly let us know your feedback.

May you have a blessed and happy Thanksgiving!

With gratitude

We control the changes in $\mathcal g_1$ and $\mathcal g_2$ individually to demonstrate the difference to our model by that change, and thus to demonstrate that different graphs correspond to different unitary of quantum circuits. We have added the above to our latest uploaded PDF, and please see our highlights in Appendix C.

Question 2: Node permutation equivariance or invariance comes naturally with graph representations. The authors need to clearly state whether the claimed generalizability, expressiveness, and efficiency should be attributed to EQGNN or VQC? What are the main technical contributions? How shall they be positioned in the context of Table 1? Also, why the mentioned reference Schatzki et al. (2022) was not included in Table 1?

A7: On your first question, I may not have fully understood what you meant. Firstly, it needs to be stated that not every VQC is permutation-equivariant. Figure R1 is a (Hardware efficient ansatz) HEA, which is an intuitive example of failure to satisfy permutation equivariance. As we can see, the output will not be the same when the order of nodes changes. In contrast, we prove the permutation equivariance of our approach in section 3.6 of the latest uploaded PDF. Moreover, if the model is not equivariant, it will overfit the variable/constraint orders in the training data, directly resulting in a decline in generalization performance and expressiveness. As shown in Table 1 of our previous paper, not every existing quantum graph neural network conforms to equivariance, and we make a comparison with QGCN in experiment section 4.8. The experiment shows that the performance of the model without equivariance will be worse. In Tabel 1, QGNN and EQGC are equivariant, but QGNN does not consider the feature encoding of nodes and edges, and EQGC does not consider the information encoding of edges, so they cannot be directly applied to the MILP task. In other words, the generalization and expressiveness brought by our model on the MILP task is unique.

On the second and third questions, Table 1 shows that our model (1) conforms to permutation equivariance, (2) can encode both node and edge features, and (3) can handle both node-level and graph-level tasks. More specifically, compared to these existing QGNNs, the main technical contributions of our approach are (1) the extension of multi-qubit encoding node features to address the high-dimensional feature coding problem, (2) the addition of an optional auxiliary layer to control model expressiveness, and (3) the association of a shared learnable parameter for each feature, and the use of repeated coding mechanisms. We have incorporated the above description into our revised paper. Thank you for your suggestions for making our article clearer.

On the last questions, the role of Table 1 is to compare existing quantum graph neural networks, but Schatzki et al. (2022) are not specialized for graph tasks. As we discuss in Appendix A, their equivariant QNNs can used to learn various problems with permutation symmetries abound, such as molecular systems, condensed matter systems, and distributed quantum sensors (Peruzzo et al., 2014; Guo et al., 2020). Moreover, they also haven't made any specific design for the high-dimensional features and edges of the graph, and they cannot be directly used for MILP tasks.

Question 4: How expressive is the proposed model?

A8: Thanks for your question. We use a tool presented in The Power of Quantum Neural Networks to quantify the expressiveness of our model. The authors introduce the effective dimension as a useful indicator of how well a particular model will be able to perform on the dataset. In particular, this algorithm follows 4 main steps:

1. Monte Carlo simulation: the forward and backward passes (gradients) of the neural network are computed for each pair of input and weight samples.
2. Fisher Matrix Computation: these outputs and gradients are used to compute the Fisher Information Matrix.
3. Fisher Matrix Normalization: averaging over all input samples and dividing by the matrix trace.
4. Effective Dimension Calculation: according to the formula in this picture (it is unable to display in this rebuttal box).

We test the proposed model with two blocks on the MILP dataset and use the random weight to evaluate the normalized effective dimension. As shown in Figure R2, with the increase in the number of data, the normalized effective dimension of the proposed model can converge to near 0.75, which shows that our model can achieve good expressiveness even with two blocks. We have added the above to the Appendix of our latest PDF.

We hope that the above answers can address your concern. If you have any further questions, please contact me again. I am looking forward to your reply.

In our model, we set the number of features encoded per qubit as $\omega$. If $\omega = 2$, each variable has four features, so it is represented by two qubits, as the node $v_0$ in Figure 5. Each constraint has two features, so it is represented by one qubit. Then, at the end of the circuit, we only measure the first qubit representing each node. The measurements corresponding to Figure 5 are shown at the end of the overall architecture diagram in Figure 3. As we can see, the measurement operation of the model acts on $q+p$ qubits, and we can obtain $q+p$ output values, where $q$ and $p$ are the numbers of decision variable nodes and constraint nodes, respectively. The proposed model can be described as a mapping $\Phi: \mathcal{G}_{q,p} \times \mathcal{H}^V_q \times \mathcal{H}^S_p \rightarrow \mathbb{R}^{q+p}$,

i.e., $\Phi(G,H) = \\{\langle 0|U_{\boldsymbol \theta}^{\dagger}(G,H)|O_i|U_{\boldsymbol \theta}(G,H)|0\rangle\\}_{i=0}^{q+p-1}$,

where ${\boldsymbol \theta}$ denotes the set of trainable parameters $(\boldsymbol \alpha, \boldsymbol \beta, \boldsymbol \gamma)$, and $U_{\boldsymbol \theta}(G,H)$ is the unitary matrix of the proposed model. $O_i$ represents $i$-th measurement, e.g., when $\omega$ is equal to 2, $O_0 = Z_0 \otimes I_1 \otimes ... \otimes I_{2q+p-1}$ indicates that Pauli-Z measurement is acted on the first qubit representing the first variable, and $O_1 = I_0 \otimes I_1 \otimes Z_2 \otimes... \otimes I_{2q+p-1}$ indicates that Pauli-Z measurement is acted on the first qubit representing the second variable, and $O_{q+p-1} = I_0 \otimes I_1 \otimes... \otimes Z_{2q+p-1}$ for the last constraint node. The output of EQGNN is defined as $\\{{\Phi_i(G,H)}\\}_{i=0}^{q+p-1}$.

For predicting feasibility, optimal value and optimal solution, we defined $\\phi_{sol}(G,H) = \\{{\\Phi_i(G,H)}\\}^{q-1}_{i=0}$,

and $\phi_\text{{fea}}(G,H) = \phi_\text{{obj}}(G,H) = \sum_{i=0}^{q+p-1} {\Phi_i(G,H)}$. As we can see, the three tasks use the same structure of EQGNN and the same measurements, but use different ways to utilize the information obtained by measurements.

For details on their loss function and evaluation metric, please see our highlights in Section 3.5 of our latest uploaded PDF, because some formulas are unable to show up in this rebuttal box.

Question 1: Is that possible for two MILP problems with different feasible/optimal solutions to be encoded into the same quantum state with the same VQC parameters?

A6: Thanks for your question. First of all, it is known that if the feasible/optimal solutions of two MILP problems are different, their MILP-induced graphs are also not the same. Then, we can answer your question by demonstrating that the proposed EQGNN can encode the two different graphs to different representations. Suppose there are two graphs $\mathcal G_1 = (G_1,H_1) =(V_1 \cup S_1, A_1, H_1)$ and $\mathcal G_2 = (G_2,H_2) = (V_2 \cup S_2, A_2, H_2)$, and the differences between the two graphs may appear in 3 places, i.e., the features of nodes, the weight and connectivity of edges. As we mentioned in Definition 4 in the latest submitted PDF, there are two types of layers in $U_{\boldsymbol \theta}(G, H)$, one is independent of the input graph, and the other is related to the input graph. In the proposed EQGNN, the feature encoding layer is $U_x^t(G,H,\alpha_t) = U_{x_1}(c,u,b) \cdot U_{x_2}(\alpha_t) \cdot U_{x_3}(l,\epsilon,\circ) \cdot U_{x_4}(\alpha_t)$, and the message interaction layer is $U^t_{g}(G,\beta) = U_{g_1}(G,\beta_{t,1}) \cdot U_{g_2}(G,\beta_{t,2}) \cdot U_{g_3} (\beta_{t,3}) \cdot U_{g_1}(G,\beta_{t,4}) \cdot U_{g_2}(G,\beta_{t,5}) \cdot U_{g_4} (\beta_{t,6})$. Among them, only $U_{x_1}(c,u,b)$ and $U_{x_3}(l,\epsilon,\circ)$ are related to the features of nodes, i.e., $H$. And $U_{g_1}(G,\beta_{t})$ and $U_{g_2}(G,\beta_{t})$ are related to the weight and connectivity of edges, i.e., $A$.

• For $A_1 = A_2, H_1 \neq H_2$, the only difference between two quantum circuits is $U_{x_1}(c,u,b)$ and $U_{x_2}(l,\epsilon,\circ)$. Taking $U_{x_1}(c,u,b)$ for example, $U_{x_1}(c,u,b) = \text{exp}(-**i**(\sum_{i=0}^{q-1}(c_i X_{2i}+u_i X_{2i+1}) + \sum_{j=0}^{p-1} b_j X_{2q+j}))$. If any feature changes, $U_{x_1}$ will change, causing $U_{\theta}(G_1,H_1)$ and $U_{\theta}(G_2,H_2)$ to be different.

• For $A_1 \neq A_2, H_1 = H_2$, there are two cases where the connections of the edges are different or the weights of the edges are different. We have already demonstrated the former in A2: Response to Weakness 2, and now we focus on the latter. $U_{g_1}(G,\beta) = \text{exp}(-**i**(\sum_{(i,j)\in \mathcal E} (A_{i,j}+\beta) Z_{2i}Z_{2q+j}))$. The change of $A_{i,j}$ will lead the change of $U_{g_1}(G,\beta)$, causing $U_{\theta}(G_1,H_1)$ and $U_{\theta}(G_2,H_2)$ to be different.

Weakness 3: In Section 4.1, it is unclear why and how the number of parameters are compared between neural networks in classic and quantum computers. The authors need to clarify how the criteria are chosen. In Section 4.2, it is also unclear why neural networks in classical computers with 2,096 and 7,904 parameters are chosen as the baseline. It will be better if the authors can compare the proposed model with a non-data-driven method for solving MILP.

A3: Thanks for your suggestion. The GNN of (Chen et al., 2023b) has one hyperparameter that controls the number of parameters, i.e., the embedding size $d$. In the tasks of predicting feasibility, the number of their parameters is $30d^2 +30d$. Our model also has a hyperparameter to control the number of parameters, i.e., the number of blocks $T$. Therefore, in our experiments, we vary these two hyperparameters separately to compare their performance results. In section 4.1, we choose $d = 4, 8, 16, 32$ as the embedding size of GNN, and $T= 2, 4, 6, 8$ as the number of blocks of EQGNN. In our experiments, we first gradually increase the embedding_size/blocks to test the performance of the models and then find the embedding_size/blocks $d^*$ or $T^*$ corresponding to the best performance. Then, we select the values near $d^*$ or $T^*$ and show their corresponding results in the experiments of the paper. In Section 4.2, for the sake of clarity in the curve graph, we only select two hyperparameters of each model for comparison. We trained EQGNN with the number of blocks at $6$, $10$, and the number of parameters is $88$ and $148$. The embedding size $d$ of GNN is chosen as $8$ and $16$ with the number of parameters at $2$,$096$ and $7$,$904$. The reasons for the setting are as follows. We observe that the train performance of GNN increases as the number of parameters increases, but the generalization performance decreases. The train performance of GNN with $d=8$ is worse than that of EQGNN, and the train performance of GNN with $d=16$ is better than that of EQGNN. Therefore, we choose the GNNs with these two hyperparameters for comparison. We have added the above description in the latest PDF. Furthermore, we have added more numerical experiments to show the performance change of the GNN with different parameters in the Appendix.

In addition, current non-data-driven methods for solving MILP, such as Gurobi or SCIP solver, can use brute-force search to find the exact solution of our datasets, because the size of our data set is relatively small. So we didn't compare them in the experiment.

Weakness 4: Throughout the paper, there are many typos/errors/inconsistencies. For example, at the beginning of the third paragraph of Section 2: "Foladable MILP instances"; in equation (4) the authors use $R_x$ gate for encoding $l$ and $\epsilon$ but in Figure 5 $R_z$ is used to encode $l$ and $\epsilon$ . It would be better to have consistent annotations for example in Section 2, the same $q$ decision variables and constraints $p$ can be used to index $V$ and $S$ in bipartite graph representation. There are many of these problems.

A4: Thanks for pointing out our mistakes. "Foladable MILP instances" should be corrected to "Foldable MILP instances", and Figure 5 should be corrected to this picture.

We have rewritten the equation (4) in our latest submitted PDF. In addition, we have modified the vertex set of bipartite graph representation as $V \cup S$, where $V=\\{v_0,\cdots,v_i,\cdots, v_{q-1}\\}$ with $v_i$ representing the $i$-th variable and $S=\\{s_0,\cdots, s_j, \cdots, s_{p-1}\\}$ with $s_j$ representing the $j$-th constraint, and $q$ and $p$ represent the number of decision variables and constraints, respectively. We apologize for not checking these errors carefully, and have checked the full text for such problems, and highlighted them in the latest submitted PDF.

Question 3: The author should clearly describe how the VQC learning problems are formulated for MILP feasibility/optimal solutions/optimal values based on the proposed EQGNN respectively either in the main text or appendix. The author should also clearly specify how the binary measurements can be used to approximate the optimal value which is often a non-discrete value. If multiple measurements are performed to recover the approximated quantum state of the VQC, the author should also clearly specify how many measurements are needed for the results provided in the experiment sections.

A5: Thank you for your suggestion. Next, we provide a detailed explanation and add it to our revised paper. First of all, it needs to be noted that our final measurement phase uses Pauli-Z measurement on the single qubit. That is, there is no need to recover the output quantum state of the entire quantum circuit by techniques like quantum state tomography, which may require an exponentially large number of measurements.

(to be continued)