A Semi-smooth, Self-shifting, and Singular Newton Method for Sparse Optimal Transport

I want to know why the authors used only MNIST and Fashion-MNIST apart from Examples 1 and 2 to compare their algorithm experimentally. There are plenty of different datasets that the authors could explore and share the results of how their algorithm behaves compared to other algorithms. Moreover, why did they choose only image data to test their algorithm in addition to Examples 1 and 2?

We appreciate the various helpful comments. Firstly, OT is not sensitive to the type of data. Its input only contains two probability vectors and a cost matrix, and all types of data are eventually summarized into these three quantities. In this sense, we are more interested in designing the input probability vectors and the cost matrix, rather than choosing the original data set. For instance, our Example 1 studies two similar distributions, and Example 2 is for two distinct distributions.

The reason for choosing MNIST and Fashion-MNIST is that they have a simple and natural cost matrix, thus convenient for demonstration, enabling us to focus on algorithm development and validation. For higher-dimensional image data, it is commonly suggested that the cost matrix should be computed from some feature space rather than from the pixel space. Therefore, for more complicated image data, typically a feature extraction step is required before applying OT, which introduces another factor that needs to be controlled. Regarding this work, we focus more on the computing of sparse OT than on the application to various data types, but of course, we would be glad to explore more data sets later, once the technical questions have been properly handled.

Dear reviewer, per the suggestion we have added more experiments on OT applied to the higher-dimensional and larger-scale ImageNet data. In particular, we compute sparse OT between two categories of images from a subset of the ImageNet data set. The subset contains ten classes of ImageNet images, and approximately 1000 images per category are selected. We map each image to a 30-dimensional feature vector by first passing the image to a ResNet18 network, resulting in a 512-dimensional vector, then followed by a dimension reduction by principal component analysis. Let $x_i\in\mathbb{R}^{30}$ be the feature vector of an image in the first category, $i=1,\ldots,n$, and $y_j\in\mathbb{R}^{30}$ be the feature vector of an image in the second category, $j=1,\ldots,m$. Then $\nu=n^{-1}\mathbf{1}_n$, $\mu=m^{-1}\mathbf{1}\_m$, and the cost matrix is $C\_{ij}= \Vert x_i-y_j\Vert^2$.

The results are displayed in Figures 19 and 20 in the revised article.

In equation 2, how is the linear system solved? Do you solve it exactly (with usual decomposition techniques) or approximately, using iterative algorithms? If you use iterative algorithms, which algorithms do you use? To which tolerance / number of iterations do you solve the linear system?

Thanks for raising this question. Equation (2) is a generic statement for solving the linear system. Specific to sparse OT, the computation method is given in Appendix C, which is based on the conjugate gradient (CG) method. Regarding the setting of precision, in our implementation, the tolerance for CG is globally set to $10^{-6}$.

The step size. First, in the main theorem 3, no clear assumption is given on the step size. I think the authors should write clearly under which assumptions the proposed algorithm converges.

Our assumption on the step size is described in Section 2.2. Specifically, we require that the step sizes are bounded from both below and above: $\eta_k\in[\tilde{m},\tilde{M}]$, $0<\tilde{m}<\tilde{M}<+\infty$.

Then, I do not understand how the condition on the step size (loosely given in the main text) can be so loose: "In fact, the only requirement is that needs to be bounded", how is it possible?? If I take a sequence of exponentially decreasing step size, how can the proposed algorithm converge? This is usually why one has to resort to costly line-search techniques. Could authors comment on this?

We apologize for the confusion on the term "bounded". A more precise statement is that we require the step sizes to be bounded from both below and above: $\eta_k\in[\tilde{m},\tilde{M}]$, $0<\tilde{m}<\tilde{M}<+\infty$. We have clarified this point in the revision to avoid confusions.

Experiments. I think experiments should be reshaped: I do not know how significant/useful the experiments with run time below a second are: I guess that for these timings, implementation matters much more than the algorithms themselves: I would remove experiments with run time below 1 second. For clarity, I would select around 5 algorithms, and display all the benchmark algorithms on a single graph, and multiple values of the regularization parameter $\gamma$ , like in Figure 1 of [1]. The block coordinate descent algorithm does not seem to converge, could you comment on this?

Thanks for the suggestions. In fact, as can be inferred from Figures 1-3, for most first-order methods, they take much longer than one second to converge, and we truncate the iteration number for better illustration. Under the same setting, Newton-type methods are able to finish in shorter times.

We think our timing results are comparable, as we spent lots of efforts in implementing every algorithm in efficient C++ code, and tried our best to reduce the overhead caused by programming language.

For the illustration of results, since we have considered many existing algorithms, it would be visually undesirable if we put all curves in the same plot. Instead, we categorize the algorithms into three groups based on their characteristics, and then compare their performance with the S5N algorithm.

For the BCD algorithm, it is normal to see a very long 'flat' stage, as this is also observed in other articles such as Figures 1 and 2 in [2]. It will converge after a sufficiently large number of iterations.

[1] Lin, H., Mairal, J. and Harchaoui, Z., 2015. A universal catalyst for first-order optimization. Advances in neural information processing systems, 28.

[2] Pasechnyuk, D.A., Persiianov, M., Dvurechensky, P. and Gasnikov, A., 2023. Algorithms for euclidean regularised Optimal Transport. arXiv preprint arXiv:2307.00321.

Reviewer h4ty,

Thank you for the discussion with the authors. May I ask if you have any further comments? Would you keep your rating or update your evaluation?

Thanks,

AC

Weakness part

Authors use grid-search to find the best step size on each step, which can be time-consuming in practice.

Thanks for the comments. We shall point out that unlike classical line search methods, the choice of the step size in S5N is optional, since the theory only requires that the step sizes are bounded from both below and above. Therefore, users are allowed to set the number of candidates, and can keep it small. Also, Algorithm 2 is not a full grid search. It immediately stops when function value decreases, and it typically only takes one or two steps.

Proposed algorithm converges globally, but it is not clear, what is the speed of global convergence.

Thanks for raising this question. To the best of our knowledge, existing global convergence results for semi-smooth Newton methods are almost all asymptotic. On the other hand, although first-order methods have non-asymptotic rates, the empirical results shown in Figures 1-3 suggest that they are much slower than Newton-type methods. We agree that the global convergence rate of semi-smooth Newton is indeed a topic that deserves future research.

Lack of definitions of some important concepts. Maybe it is clear for those, who are very familiar with this topic, but I think, introduction of such definitions may improve the overall look of papers (see questions).

Thanks for the suggestions. We have added the relevant definitions in Appendix A in the revision.

Question part

Probably, in the abstract you have a typo and you meant "possibly non-smooth gradients" instead of "possibly non-differentiable gradients"

Thanks for pointing out this issue. In our context, "non-smooth" and "non-differentiable" have the same meaning. For consistency we have changed "non-differentiable" to "non-smooth" (in the PDF, as we are no longer able to edit the abstract shown in the page).

Please, increase the font size of text in the figures

Please, provide the definitions of local Lipschitz continuity and upper semi-continuity.

Thanks for all these helpful suggestions. We have updated the definitions related to local Lipschitz continuity and upper semi-continuity in Appendix A, and have enlarged the font size of the text in the figures.

First of all, we thank the reviewer for the comments, and would like to respectfully correct a factual misunderstanding in the review.

It seems that $f(x)$ in Eq. (4) is not a smooth function because $()_+$ operator is not smooth. Thus, Proposition 1 can not hold since $\nabla f(\alpha, \beta)$ does not exist.

The convergence analysis of this paper lies on Assumption 1 and Assumption 2 which requires that $f(x)$ is differentiable, so the convergence analysis can not be used in the problem that $f(x)$ is not smooth.

This paper claims that S5N is used to solve problems that have possibly non-differentiable gradients and non-isolated solutions. However, the convergence analysis requires Assumption 1 which supposes that $f(x)$ is differentiable.

We would like to point out that $f(x)$ in Eq. (4) in fact involves $h(x)=[(x)\_+]^2$, not $(x)\_+$. It can be verified that the former is differentiable, with the derivative $h'(x)=2\cdot(x)_+$. Therefore, $\nabla f(\alpha,\beta)$ in Proposition 1 indeed exists, and Assumption 1 holds for the problem we consider.

We humbly hope the reviewer could re-evaluate our main results and raise the score when appropriate, if this misunderstanding is resolved.

The results in Sec. 2.4 seems not new. They seem to be standard results of semi-smooth Newton. Thus, the author should cite some references.

Thanks for the comments. While the expressions of the main conclusions in Section 2.4 appear to be standard, we have proved them for our new algorithm, S5N, under completely new settings (non-smooth gradient and singular Hessian). Therefore, the proof is specific to S5N, which is not a trivial consequence of existing results.

We have added references to standard results, and have made our contributions clearer in the revision.

It is quite counterintuitive for me: why do we want to use second-order methods for almost quadratic problems (quadratic + linear)? Why not use simple CG or PCG to solve it?

We would like to point out that this problem is piecewise quadratic, which has properties quite distinct from quadratic functions. In particular, its KKT condition is not a linear system, so CG or PCG is not applicable here.

We focus on second-order methods here since (1) second-order methods in general exhibit fast convergence; (2) for sparse OT problems, the generalized Hessian is sparse, making the per-iteration cost very small.

Why in time plots the first-order methods are slower than S5N?

There are two major factors: (1) Newton-type methods have faster convergence speed in terms of iteration number; (2) according to our analysis in the response above, first-order methods might even have a higher per-iteration cost. Overall, S5N exhibits visible advantage over first-order methods.

Why do you not compare S5N with the classical baseline for OT, the Sinkhorn method?

Thanks for raising this question. Given that sparse OT is a highly structured problem, and Sinkhorn represents another approximation to OT that results in a dense optimal plan, our focus is on the computation of quadratically regularized OT. It is not comparable to Sinkhorn when sparsity is required.

Why does the gradient time computation differ for different methods? "the semi-dual L-BFGS has a fast convergence in iteration number, it suffers from a long run time since its gradient computation is more difficult and time-consuming than other methods."

Although L-BFGS can be applied to both the dual and semi-dual objective functions, their gradient computation is different. Semi-dual L-BFGS has a similar computation as BCD, and according to our analysis above, it takes $\mathcal{O}(n^2\log n)$ time in each gradient evaluation. In contrast, dual L-BFGS is $\mathcal{O}(n^2)$. Therefore, although the total number of iterations for semi-dual might be smaller than dual, the longer per-iteration time of semi-dual results in a higher overall run time.

Do you have any intuition: why “ASSN has a very slow convergence” in practice? It seems to have the same structure with gradient regularization.

We have two conjectures here: (1) although algorithmically ASSN can be applied to the QROT problem, the nature of the QROT objective function does not satisfy the convergence assumptions of the algorithm (fixed-point operator, singular Hessian, etc.); (2) from the perspective of ASSN's structure, it has no flexibility in choosing the step size. As a result, we have observed that the algorithm undergoes many rejection steps (i.e., Newton steps that do not behave well) during iterations, resulting in a slow convergence speed.

Relation to papers [1], [2]; comparison to first-order methods; improvement and theoretical challenge over the previous papers.

Thanks for the comments. Regarding the two articles mentioned above, we would like to make the following clarifications:

1. [1] introduced the ASSN algorithm, which only applies to fix-point problems of the form $T(z)=z$, where $T=(1-\alpha)I+\alpha R$, and $R$ is a nonexpansive operator. This imposes a restriction on the problem to solve, whereas in contrast, S5N has substantially broader use.
2. In addition, the local convergence result of [1] relies on the BD-regular condition, which essentially does not support singular problems. We have reflected this point in Table 1.
3. Also shown in Table 1, [2] assumes that the objective function is twice continuously differentiable and the Hessian is Lipschitz continuous. These are quite strong assumptions, and that is why it can achieve local cubic convergence.

In this sense, neither [1] nor [2] really applies to the general setting of non-smooth gradient and singular Hessian. In constract, we have obtained both global and local convergence properties of S5N under significantly weaker assumptions, and we think this is a major theoretical contribution that differs from existing literature.

As for first-order methods, although they have non-asymptotic global convergence rates, they are mostly sub-linear, thus resulting in slower convergence speed as evidenced in Figures 1-3.

More importantly, the per-iteration cost of first-order methods is not necessarily smaller than that of Newton-type methods for sparse OT problems. For example, block coordinate descent (BCD) needs to sort length-$n$ vectors $n$ times, with a total cost of $\mathcal{O}(n^2\log n)$. Newton-type methods constructs the $\sigma$ sparse matrix in Proposition 1 in $\mathcal{O}(n^2)$ time, and solves the linear equation using conjugate gradient (CG) in $\mathcal{O}(nnz\cdot\log(1/\varepsilon))$ time, where $nnz$ is the number of nonzero elements in $\sigma$, and $\varepsilon$ is the tolerance. Their total per-iteration cost can be substantially smaller than $\mathcal{O}(n^2\log n)$.

In the optimal transport (OT) section, the authors describe the OT problem and why it has the semi-smooth Hessian, as discussed in [3].

Again shown in Table 1, the algorithm used by [3], named GRSSN, does not provide proofs for the global and local convergence guarantees. Crucially, GRSSN introduces a very sensitive hyperparameter $\lambda$, which we have discussed in Figures 1-3, Section 5. S5N is shown to perform better than GRSSN with no hyperparameter tuning.

Numerical Experiments

To compare different algorithms, we not only consider the iteration number, but also the run time as a metric for the overall performance. We have implemented every algorithm in efficient C++ to minimize the overhead and to make the run times comparable. Also, the theoretical complexity analysis from our response above also indicates that second-order methods may even have a smaller per-iteration cost than some first-order methods in the context of sparse OT.

[1] Xiantao Xiao, Yongfeng Li, Zaiwen Wen, and Liwei Zhang. A regularized semi-smooth Newton method with projection steps for composite convex programs. Journal of Scientific Computing, 76:364–389, 2018.

[2] Weijun Zhou and Xinlong Chen. On the convergence of a modified regularized Newton method for convex optimization with singular solutions. Journal of Computational and Applied Mathematics, 239:179–188, 2013.

[3] Dirk A Lorenz, Paul Manns, and Christian Meyer. Quadratically regularized optimal transport. Applied Mathematics & Optimization, 83(3):1919–1949, 2021.

Reviewer tikC,

Thank you for the discussion with the authors. May I ask if you have any further comments? Would you keep your rating or update your evaluation?

Thanks,

AC