Thanks for your response! We provide our responses as follows:

**Q1. In general, you won't be able to eliminate a term by using a different derivation method. Your method and OM essentially derivate the same variational lower bound, so the term needs to be somewhere in your result.**

We want to clarify our VLB differs from OM's more than just eliminating one term. There are three differences between the two VLBs: (1) Ours does not need to compute $|\mathit{\Pi}(G_{1:n})|$; (2) Ours does not include the entropy of node ordering; (3) Ours has an additional KL divergence term between the diffusion ordering and the generation ordering.

The VLBs of the two methods are different because their modeling processes are different. The key is how $p_\theta\left(G_{1: n}, \sigma_{1: n}\right)$ is treated. In OM, $p_\theta\left(G_{1: n}, \sigma_{1: n}\right)$ is factorized as: $$p_\theta\left(G_{1: n}, \sigma_{1: n}\right) = p_\theta\left(G_{1: n}\right) p_\theta\left(\sigma_{1: n} \mid G_{1: n} \right).$$

As $p_\theta\left(\sigma_{1: n} \mid G_{1: n} \right)$ is a uniform distribution, OM inevitably needs to involve $|\mathit{\Pi}(G_{1:n})|$ in their VLB. In contrast, in our method (please refer to line 6 in the derivation of VLB in our previous response), we factorize $p_\theta\left(G_{1: n}\right)$ as:

$$ p_\theta\left(G_{1: n}, \sigma_{1: n}\right) = p_\theta\left(G_{1: n-1} \mid \sigma_{1: n}\right) p\left(\sigma_{1: n} \mid G_n\right) p\left(G_n\right).$$

This factorization leads to the following consequences / benefits:

- The challenge of modeling the node ordering variable is now pushed into the term $p\left(\sigma_{1: n} \mid G_n\right)$. This is a key reason that we do not have the $\left|\Pi\left(G_{1: n}\right)\right|$ term in our VLB, instead we have a KL term. Now how do we learn this conditional distribution? It turns out we do not need to! During training, the generation ordering is set to the exact reverse of the diffusion ordering and thus the KL-divergence term is equal to 0. At test time, the generator network will implicitly result in an ordering matching the diffusion ordering, and the KL-divergence term tends to approach 0 (elaborated in Q5). In other words, the node ordering network $q_\phi\left(\sigma_{1:n} \mid G_0\right)$ itself is a proxy for $p\left(\sigma_{1: n} \mid G_n\right)$ and the KL term will approach zero. This is why we do not have the entropy term and the term $\left|\Pi\left(G_{1: n}\right)\right|$ in our VLB, and meanwhile we can also ignore the KL term but only optimizing the first term in Eq 3.

<!--
Due to the symmetry between the forward and backward diffusion processes, the generation ordering $p\left(\sigma_{1:n} \mid G_n\right)$ will be the same as the diffusion ordering $q_\phi\left(\sigma_{1:n} \mid G_0\right)$. In other words, the node ordering network $q_\phi\left(\sigma_{1:n} \mid G_0\right)$ itself is a proxy for $p\left(\sigma_{1: n} \mid G_n\right)$ and the KL term will approach zero. This is why we do not have the entropy term and the term $\left|\Pi\left(G_{1: n}\right)\right|$ in our VLB, and meanwhile we can also ignore the KL term but only optimizing the first term in Eq 3.
-->

- Still because of the factorization, our generator network $p_{\theta}(G_{0:n-1}|G_n,\sigma_{1:n})$ models the graph generation **conditioned on** a fixed node ordering. In other words, given an oracle ordering of how the nodes should be generated (which is given by the diffusion ordering network $q$), we use the generator network to model how the graph sequences are generated. In contrast, OM models $p(G_{1:n})$ with no exact node ordering information (elaborated in Q3). Our conditional generator $p_{\theta}(G_{0:n-1}|G_n,\sigma_{1:n})$ is further factorized by autoregressive generation: $\log p_{\theta}(G_{0:n-1}|G_n,\sigma_{1:n}) = \sum_{t=0}^{n-1} \log p_\theta\left(G_t\mid G_{t+1}, \sigma_{1: n}\right).$

Deriving the lower bounds of the same objective from different perspectives can indeed lead to different results. This also happens in the vallina VAE and diffusion model. VAE and diffusion model are both deriving a lower bound of $p(\mathbf{x})$. However, their VLBs are different because their latent variables are different. Similarly, this happens between our method and OM. The modeling processes of the two are different, and their VLBs are different.

**Q2. What does $G_t(\sigma_{1:n})$ exactly mean?**

We used the notation $p_\theta\left(G_t\left(\sigma_{1: n}\right) \mid G_{t+1}\left(\sigma_{1: n}\right)\right)$ to mean the autoregressive generation is conditioned on a node ordering. Namely,
$p_\theta\left(G_t\left(\sigma_{1: n}\right) \mid G_{t+1}\left(\sigma_{1: n}\right)\right) \equiv p_\theta\left(G_t \mid G_{t+1}, \sigma_{1: n}\right)$.
We can change this notation in the final version to avoid confusion.

**Q3. Compared with the OM paper, which works on adjacency matrices, your conditional probability $p_{\theta}(G_t(\sigma_{1:n})|G_{t+1}(\sigma_{1:n}))$ seems to correspond the probability of $p(A_{1:(n-t)}|A_{1:n-t-1})$ but differ by the probability of a particular node leading the growth of the adjacency matrix from $A_{1:(n-t),1:(n-t)}$ to $A_{1:(n-t-1),1:(n-t-1)}$. Because of automorphism, the probability is one over the number of nodes that lead $A_{1:(n-t-1),1:(n-t-1)}$ to $A_{1:(n-t),1:(n-t)}$ . Because of this term, you don't have $\text{log}|\mathit{\Pi}(G_{1:n})|$, while OM's result has.**

$p_{\theta}(G_t|G_{t+1},\sigma_{1:n})= p_{\theta}(\{e_{\sigma_t,\sigma_j}\}_{j=t+1}^{n}|\sigma_{1:n})$,

The physical meaning of $p_\theta\left(G_t\left(\sigma_{1: n}\right) \mid G_{t+1}\left(\sigma_{1: n}\right)\right) \equiv p_\theta\left(G_t \mid G_{t+1}, \sigma_{1: n}\right)$ is different from $p(A_{1:(n-t),1:(n-t)}|A_{1:(n-t-1),1:(n-t-1)})$ and the probability of a particular node ordering leading the growth of the adjacency matrix. It corresponds to the probability of the edges between the node $v_{\sigma_t}$ with the set of unmasked nodes that have been generated so far $\{v_{\sigma_j}\}_{j=t+1}^{n}$, namely $p_{\theta}(G_t|G_{t+1},\sigma_{1:n})= p_{\theta}(\{e_{\sigma_t,\sigma_j}\}_{j=t+1}^{n}|G_{t+1},\sigma_{1:n})$, where we use $e_{\sigma_t,\sigma_j}$ to denote the edge between node $v_{\sigma_{t}}$ and $v_{\sigma_{j}}$. Note that this is different from the growth of the adjacency matrix, as we illustrate below.

In OM, it considers the adjacency matrix $A_{1:n,1:n}$ by autoregressively modeling it as $p(A_{1:n,1:n})=\prod_tp(A_{1:(n-t),1:(n-t)})|A_{1:(n-t-1),1:(n-t-1)})$ using a neural network. However, the node ordering information is lost in $A_{1:n,1:n}$ since multiple ordering can lead to this matrix. That is why OM needs to compute automorphism to obtain $p(A_{1:n.1:n},\sigma_{1:n})=p(\sigma_{1:n}|A_{1:n,1:n})p(A_{1:n,1:n})=\frac{1}{|\mathit{\Pi}(A_{1:n,1:n})|}p(A_{1:n,1:n})$.

In constrast, we factorize $p_{\theta}(G_{1:n},\sigma_{1:n})$ as $p_\theta\left(G_{1: n}, \sigma_{1: n}\right) = p_\theta\left(G_{1: n-1} \mid \sigma_{1: n}\right) p\left(\sigma_{1: n} \mid G_n\right) p\left(G_n\right)$, where our generation procudure $p(G_{1:n}|\sigma_{1:n})$ is **conditioned on** a fixed node ordering $\sigma_{1:n}$. In our framework, the node ordering $\sigma_{1:n}$ is sampled in the forward diffusion process which essentially assigns a unique ID to each node of $G_0$. In the generation procedure, our starting point is the masked graph $G_n$ with $n$ nodes. Our generator network sequentially denoises the node $v_{\sigma_t}$ in $G_n$ in the reverse order of the diffusion ordering. That is, at each generation step, it computes the connectivity of node $v_{\sigma_t}$ with all unmasked nodes that have been generated so far:
$p_{\theta}(G_t|G_{t+1},\sigma_{1:n})=p_{\theta}(\{e_{\sigma_t,\sigma_j}\}_{j=t+1}^{n}|G_{t+1},\sigma_{1:n})$. As for the generation ordering $p(\sigma_{1:n}|G_n)$, we do not need to model it since the diffusion ordering network is a proxy (illustrate in Q1 and Q5). Therefore, by our new factorization and modeling process, we eliminate the need of computing the graph automorphism.

**Q4. In your implementation, your network seems to also compute $p(A_{1:(n-t-1),1:(n-t-1)}|A_{1:(n-t-1),1:(n-t-1)})$ as other autoregressive models do, so your derivation of the lower bound does not seem to match your implementation.**

Our derivation of the lower bound matches our implementation. At each time step $t$, our network predicts the connectivity of the node $v_{\sigma_t}$ with all the unmasked nodes that have been generated so far, i.e., $p_{\theta}(G_t|G_{t+1},\sigma_{1:n})= p_{\theta}(\{e_{\sigma_t,\sigma_j}\}_{j=t+1}^{n}|G_{t+1},\sigma_{1:n})$.

The physical meaning is different from $p(A_{1:(n-t-1),1:(n-t-1)}|A_{1:(n-t-1),1:(n-t-1)})$ as we have illustrated in Q3. $A_{1:(n-t-1),1:(n-t-1)}$ is the partial adjacency matrix with no exact node ordering information. While our network denoises the node $v_{\sigma_t}$ in the masked graph $G_t$ with $n$ nodes
in the reverse order of the diffusion ordering.

**Q5. There is a separate issue in your derivation: $p(\sigma_{1:n}|G_n)$ is unknown. According to your model, this distribution cannot be defined: it needs to be derived by marginalizing out $G_{0:n-1}$ from $p(G_{0:n-1},\sigma_{1:n}|G_n)$. By an argument of symmetry, it seems that the probability $p(\sigma_{1：n}|G_n)$ is $\frac{1}{n!}$, but it needs formal proof. The submission just says "it is the distribution of the generation ordering" without a concrete calculation.**

$p(\sigma_{1:n}|G_n)$ is the denoising/generation ordering of the generator network. During training, the generator network denoises/generates the nodes in the graph $G_n$ using the exact reverse order of the diffusion ordering $\sigma_{1:n}$. Therefore, in training, $p(\sigma_{1:n}|G_n)$ is equal to $q(\sigma_{1:n}|G_0)$; the KL divergence term is equal to 0 and we can ignore it.

Now the question is, at test time, how do we get $p(\sigma_{1:n}|G_n)$ as we do not have a diffusion ordering? The answer is we do not need to, because the autoregressive procedure of the generator network will implicitly result in such an ordering. During training, the training objective $E_{q(\sigma_{1:n}|G_0)}\sum_{t}\log p(G_t|G_{t+1},\sigma_{1:n})$ encourages the generator network to denoise/generate the nodes in the reverse ordering of the diffusion process. Therefore, at test time, the generation ordering will be close to the diffusion ordering of the target graph, and the KL-divergence term tends to approach 0 (we have shown this in Appendix A.9 of the updated manuscript). Since the diffusion ordering learns an data-dependent node ordering, the matching generator network can also result in generation orderings that reflect graph topology/regularity.