Weaknesses

Thanks for the comments. First of all, we shall clarify that SSNS is directly inspired by the regularized Newton method [R6] and the Levenberg–Marquardt algorithm [R7,R8], which have close relations to the trust-region methods but with some subtle differences. Specifically, standard trust-region methods first set a trust-region radius $\Delta_k$ in each iteration, and then (approximately) solve a subproblem to determine the search direction. In SSNS, we update the shift parameter $\lambda_k$ instead, and the search direction has a closed form. The formula of search directions in trust-region methods has a similar form, but its shift $\lambda_k$ is implicitly determined by $\Delta_k$, and typically requires some root-finding techniques. Also, SSNS allows for a step size selection procedure as in Algorithm 3, thus enhancing the flexibility.

Also per the suggestions of reviewers, we have improved the numerical experiments as explained in the global rebuttal.

Questions

1. Recall that in Theorem 2 we allow $\delta=0$, which corresponds to the case that we use the genuine Hessian matrix. In this setting, Algorithm 2 reduces to a conventional regularized Newton method, but its per-iteration cost would be huge. Theorem 2 is an improvement to existing sparsifed Newton methods such as [R3], since previous works do not guarantee that the sparsified Hessian is invertible (but note that the true Hessian matrix is always invertible by Theorem 2).
2. In our original design of the algorithm, "safe" mainly refers to the positive definiteness of the sparsified Hessian matrix. But from the experiments, we also find that the shift parameter (which has a close relation to trust-region techniques) is very useful, as the classical Newton method fails in some examples due to numerical instability.
3. Thanks for the suggestion. We have added such information in our local revision. Basically it is $n=m=784$ for (Fashion-)MNIST and $n\approx m\approx 1000$ for ImageNet data.
4. Yes. This can be inferred from the fact that the number of iterations of the Newton method is very small, but its runtime is large. The huge per-iteration cost is the consequence of a dense linear system.
5. Yes. We have defined this notation in the beginning of Section 2, meaning removing the last column of $A$.

[R3] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., & Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse Newton iterations.

[R6] Li, D. H., Fukushima, M., Qi, L., & Yamashita, N. (2004). Regularized Newton methods for convex minimization problems with singular solutions.

[R7] Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares.

[R8] Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters.

Thank you for the comments. Below are our point-by-point responses for the questions.

Weaknesses

1. To enhance the reproducibility, we do not pick the pairs of images. Instead, the randomly selected image IDs are taken from the prior literature [R1] that studies quadratically regularized OT, and we follow their experiment setting to study entropic-regularized OT. In fact, we have already included other image pairs in Figures 4 and 5 in Appendix A.3, and the IDs of image pairs can be seen from the titles of images. In the rebuttal PDF we have added more test cases in Figure R1.

2-4. Thanks for the suggestions. We have fixed these issues in our local revision.

5. Thanks for the suggestion. We call our proposed algorithm a safe and sparse Newton method, in the sense that the linear systems for computing the search directions are always positive definite. This property addresses the invertibility issues of existing sparsified Newton methods, and is crucial for practical implementation. In our local manuscript, we have added explanations in the introduction, above the contribution section.

Questions

We are trying to demonstrate two typical uses of OT. In the first experiment, images are vectorized as density vectors, and OT is a used as a tool for image morphing and interpolation [R2]. The experiment setting is already used in previous literature such as [R1] and [R3]. The second experiment uses OT as a statistical distance to measure the difference between two distributions. In this setting, each image is one observation of a distribution, and we use OT to compute the (approximate) Wasserstein distance between two classes of images. In our local revision, we have added a section in the appendix to explain the motivation of these experiments.

[R1] Pasechnyuk, D. A., Persiianov, M., Dvurechensky, P., & Gasnikov, A. (2023). Algorithms for Euclidean-regularised optimal transport.

[R2] Papadakis, N. (2015). Optimal transport for image processing.

[R3] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., & Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse newton iterations.

Thanks for the comments. Below are our point-by-point responses.

Weaknesses

1. The dimension of the features does not impact the scale the problem. The extracted features are only used to compute the cost matrix, which has a size of $n\times m$. We do not make the feature dimension $d$ too large mainly because it is commonly accepted that the Euclidean distance suffers from the curse of dimensionality. In the global rebuttal, we have added experiments to study the impact of feature dimension and the scalability of algorithms on very large problems.
2. As we have mentioned in the global rebuttal and other threads, to enhance the reproducibility, we do not pick the pairs of images. Instead, the randomly selected image IDs are taken from the prior literature [R1] that studies quadratically regularized OT, and we follow their experiment setting to study entropic-regularized OT. We have included different examples in Figures 4 and 5 in Appendix A.3. Also, to investigate the scalability of algorithms on very large problems, we have added experiments to test the performance for $n=m=1000,5000,10000$.
3. We agree that the overall framework of SNS [R3] is easy to understand, but the major challenges for directly comparing with SNS are: (a) there is no publicly available code of SNS to reproduce the experiments; (b) in [R3], the algorithm has many hyperparameters, and we have no clear scheme to set those hyperparameters. For example, how many Sinkhorn iterations to run before switching to the Newton step, how large the sparsification parameter should be, what the fallback strategy should be when the linear system is not invertible, etc. In fact, one of the major motivations of this article is to make the SNS framework more adaptive and practical, and we position our SSNS algorithm as a concrete and practical variant of the original SNS method.

Questions

1. A trivial upper bound of the computational complexity of Algorithm 1 is $O(n^2\log(n))$, assuming $n=m$. This is done by sorting the values of each column or row of $T$. However, if we have an estimated upper bound $\varrho$ of the density of the sparsified $T$, meaning that each column or row of $T$ has at most $\varrho n$ nonzero elements after sparsification, then the cost can be reduced to $O(n^2+\varrho n^2\log(\varrho n))$. This is done by selecting the largest $\varrho n$ elements in each column or row of $T$, and only sort values within these $\varrho n$ elements. In a typical setting, $\varrho$ is very small due to the approximate sparsity of $T$.
2. Each iteration of Algorithm 2 has a computational cost of $O(n^2)$ to compute the gradient and objective function value, plus the cost in Algorithm 1, plus the cost in solving the sparse linear system. When $T$ is very sparse, the number of nonzero elements of $T$ may be as small as $O(n)$, so the overall computation is dominated by the first and second parts, i.e., $O(n^2)$.
3. As we have explained above, there are certain practical difficulties in directly comparing with SNS. But we have tried our best to provide more experiment results similar to those in the SNS paper, as introduced in the points below.
4. A typical definition of the optimality gap is the value $f(x_k)-f(x^*)$, where $f(\cdot)$ is the objective function, and $x^*$ is the optimal point. The marginal error in entropic-regularized OT coincides with the gradient norm $\Vert g(x_k) \Vert$. We do not use the optimality gap because the ground-truth optimal point $x^*$ is in general unknown, and needs to be approximated using existing solvers. But this process also introduces rounding errors, and may be very sensitive when $x_k$ is close to convergence. Instead, $\Vert g(x_k) \Vert$ is easy to compute, and has an absolute lower bound of zero. Overall, we think that the marginal error is a more objective and easy-to-compute criterion to evaluate convergence.
5. Yes, we have added such experiments in the global rebuttal.
6. Per the suggestion, we have added such experiments in the global rebuttal (Figure R2).

[R1] Pasechnyuk, D. A., Persiianov, M., Dvurechensky, P., & Gasnikov, A. (2023). Algorithms for Euclidean-regularised optimal transport.

[R3] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., & Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse newton iterations.

Weaknesses

Thanks for the comments. As we have explained in the global rebuttal and other threads, the experiments are intentionally designed to reflect two typical uses of OT, one for image morphing and interpolation, and the other for computing statistical distances. In our local revision, we have added a section in the appendix to explain the motivation of these experiments.

As for the scalability, we have added a new set of experiments at the scales of $n=m=1000,5000,10000$, as explained in the global rebuttal. The results show that SSNS is capable of solving entropic-regularized OT on very large problems, and is more efficient than BCD in terms of run time and number of iterations.

Questions

1. We suppose the "runtime" here refers to the theoretical computational complexity based on the global convergence speed. First of all, to the best of our knowledge, most convergence rate results for Newton-type methods are local, and there is few on the global convergence rate. However, there is a simple "ensemble" method to safeguard the global convergence speed. Note that our Algorithm 2 never increases the objective function value in any iteration, so we can append a Sinkhorn iteration after every $N$ Newton-type steps, where $N$ is a fixed integer. In this way, we obtain at least a Sinkhorn-like convergence speed, but are likely to make larger progresses during the Newton steps. This idea is also mentioned in Section 3.2, page 41 of [R4].
2. Following the formulation in [R5], entropic-regularized unbalanced OT also has a smooth dual objective function (equation (4) of [R5]). Its difference with the balanced version is that the linear term in equation (4) of our article, $-\alpha^T a-\beta^T b=-\sum_{i=1}^n \alpha_i a_i-\sum_{j=1}^m \beta_j b_j,$ is replaced by sums of exponentials, $\tau \sum_{i=1}^n e^{-\alpha_i/\tau} a_i+\tau \sum_{j=1}^m e^{-\beta_j/\tau} b_j.$ Clearly, this means that only the diagonal elements of the Hessian matrix would be different from the setting in our article, and we think that the proposed sparsification scheme and the Newton-type algorithm still apply.

[R4] Nocedal, J., & Wright, S. J. (2006). Numerical optimization.

[R5] Pham, K., Le, K., Ho, N., Pham, T., & Bui, H. (2020). On unbalanced optimal transport: An analysis of Sinkhorn algorithm.