We greatly appreciate the reviewers insightful comments and suggestions. To begin, we are thankful for the suggested references in W3, and will add these to the manuscript.

Regarding Q1, if the score matrix is close to a permutation $P^\*$, then the rounding scheme can indeed produce a solution better than $P^\*$. All that is required for this to happen is that some permutation in the decomposition is better than $P^\*$. While Property 4-2 gives that at least one permutation will be $P^\*$, there is nothing preventing the decomposition from containing better permutations.

Theorem 2.11 shows the existence of an $S$ that allows for escaping a local minimum. However, the proof is not constructive, as, indeed, the score $S$ used in the proof is a permutation minimizing the permutation function $f$. Thus, the hardness of finding this $S$ is the same as optimizing $f$. However, as suggested in Q2, an interesting avenue for future research would be to investigate whether there are score matrices with desirable theoretical properties that are not computationally hard to find. For now, Theorem 2.11 only shows that the approach of changing $S$ always has the potential to escape local minima but does not provide any procedure for finding such an $S$.

We thank the reviewer for their constructive comments. We have performed a variety of ablation studies in response.

Including evaluation of the dependence on truncation and score update frequency is a great suggestion. We have performed these experiments and will update the manuscript to include curves showing these dependencies. Our results show that, for $k \in \{3,5,10,15,20\}$ the greatest difference in performance is about 0.88% in terms of avg. gap from optimal values. Furthermore, the difference from $k = 5$ to $k = 10$ (which achieves the best performance) is about 0.24%. Note that we do not include the plots here due to restrictions from Neurips on the response format.

We also performed ablation experiments on score update frequency that we will amend the paper to include. In these experiments update_score is set to True every $m\in \{1,5,20,10,50\}$ epochs. These experiments show that the performance improves with higher frequency updates, that is, smaller $m$. Therefore, the best choice is always set update_score = True.

Despite the $O(n^3)$ complexity (with truncation), our method does scale to large instances. For example, for QAP, we test on QAPLIB which has an instances of size $n =256$. For this problem we get within 0.8% of the best known solution using a runtime of 417 seconds.

We agree that additional analysis on the truncated case would be beneficial. To this end, we will introduce a section formally introducing the truncated extension $F_S^{(k)} := \frac{1}{\sum_{i=1}^{k} \alpha_i} \sum_{i=1}^{k}\alpha_iP_i$. proving that all rounding properties (Properties 3, 4-1, and 4-2) hold in the presence of truncation and renormalization of coefficients.
This result follows easily since these properties only utilize characteristics of the first term in the decomposition and that extension is formed from a convex combination of discrete objective values. We will also add a proof of Theorem 2.11 in the truncating case. Again this is straightforward as the proof of this theorem only requires evaluating properties of the first term in the decomposition (ee duplicate comment to reviewer Qomh).

The reviewer brings up a valuable point that there are many ways to optimize over the Birkhoff polytope: Frank--Wolfe, gradient descent + Sinkhorn, or, potentially, projected gradient descent. We performed experiments with both Frank--Wolfe and gradient descent + Sinkhorn, and found that Frank--Wolfe offered superior optimization results. In particular, using Frank-Wolfe yields an avg. gap from optimal of 4.14%, whereas gradient descent + Sinkhorn gives a gap of 4.55%. We will edit the manuscript to include these in the appendix.

For each of these Birkhoff-polytope optimization methods, an objective over the polytope is required. Although there might be many options for this objective, we propose Birkhoff extension for its several advantageous properties, including differentiability and guaranteed rounding. For a comparison with an alternative objective, see the TSP experiments, where we evaluate optimization over a quadratic program relaxation. Here we demonstrate that superior results are achieved by initially optimizing the QP relaxation and then using this solution as the score matrix initialization. We invite further comments if there is a more direct objective we should consider.

We appreciate the thoughtful review and now respond to points raised.

We agree that our results could be strengthened by closing the gap between theory and practice created by truncation. To this end, we will introduce a section formally introducing the truncated extension $F_S^{(k)} := \frac{1}{\sum_{i=1}^{k} \alpha_i} \sum_{i=1}^{k}\alpha_iP_i$, proving that all rounding properties (Properties 3, 4-1, and 4-2) hold in the presence of truncation. This result follows easily since these properties only utilize characteristics of the first term in the decomposition and that extension is formed from a convex combination of discrete objective values. We will also add a proof of Theorem 2.11 in the truncating case. Again this is straightforward as the proof of this theorem only requires evaluating properties of the first term in the decomposition (see duplicate comment to reviewer S7f4).

Thank you for suggesting the need for additional clarity with regards to $\mathcal P$. We define (Definition 2.2) $\mathcal P$ only as an operator on matrices $X$ with non-negative entries. For such matrices $\mathcal P(X)$ is the set of all permutations $P$ such that $X(i,j) >0$ if $P(i,j) = 1.$ This set is different than the set of permutations appearing in the Birkhoff decomposition; it contains all permutations that could potentially appear in a Birkhoff decomposition of $X.$ Line 219 discusses finding the permutation in $\mathcal P(X)$ that comes first in the total order. This operation does not require computing the decomposition, so there is not a remainder to address.

Regarding the necessity of random initialization, since the extension is a.e. differentiable, random initialization ensures we initialize at a point that is differentiable (almost surely). This approach also has the advantage that optimizing from multiple initializations may yield superior results. However, it is very interesting to consider starting from a fixed doubly stochastic matrix. In this setting we can use multiple score matrix initializations to get multiple optimization trajectories. We performed experiments showing the effectiveness of initializing at the barycenter. We picked 14 problem instances from the QAPLIB with sizes from $n = 12$ to $n = 35$. We use identical hyperparameter (truncation depth $k = 5$, learning rate $\eta = .001$, time budget $t = 2n$ seconds) and, for random and barycenter initialization, we evaluate the final average assignment costs across all 14 problem instances. These experiments show that random initialization version has avg. gap to optimal of 4.14%, whereas the barycenter initialization has an avg. gap of 4.42%. We will include these experiments in the manuscript.

We do, indeed, use automatic gradients in practice and will include a comment stating this in the manuscript. As the expression for the $k$th term in the decomposition includes $k$ nested min functions, the gradient is quite challenging to express.

We are thankful for the especially thorough review, and glad the reviewer has a strong understanding of our work. We now address the questions raised.

The reviewer's suggestion of warm starting the decomposition for $X$ using the previous GD iterate $X'$ is very interesting. In an extreme case, e.g. if $X'$ is a permutation matrix, the previous iterate's decomposition may contain very little information about the current decomposition (as the previous decomposition only has a single element). Therefore it seems challenging to reduce the worst-case complexity from $O(n^3)$ using this approach. On the other hand, in most cases the decompositions of $X$ and $X'$ should be very similar, so using a warm-start could offer a substantial practical speedup. This is an encouraging direction for future improvements. We have explored this observation to some degree in our implementation, where we have made the engineering choice to use a hash table that saves the results of common maximum weight matching problems. Using this technique, any maximum weight matching problems solved for computing the decomposition in the previous iterate will not be recomputed if these problems arise in computing the new decomposition.

It is compelling to consider, as the reviewer suggests, how the properties of the extension $F_S$ are dependent on $S$. We note a bound on the Lipschitz constant of the extension need not depend on $S$. The relevant quantity is instead $\Delta_f = \max_P f(P) - \min_Pf(P)$. Roughly, this arises because if two doubly stochastic matrices $A$ and $A'$ are within an $L_1$ distance $\epsilon$, then the $L_1$ difference in the weights $\alpha_i$ of their decompositions is $O(\epsilon)$ (note that each coefficient is a 1-Lipschitz function of the doubly stochastic matrix). Thus, in the worst case, $|F(A) - F(A')|$ is $O(\epsilon \Delta_f)$, where we have suppressed the dependence on $n$. We will add a rigorous proof to the manuscript showing a Lipschitz bound independent of $S$. There may be other interesting ways in which the extension is dependent on $S$. These are compelling avenues for future work.