The Sparse-Plus-Low-Rank Quasi-Newton Method for Entropic-Regularized Optimal Transport

Thank you for your detailed review and thoughtful comments. Please see our point-by-point responpses below.

---

Not always, as ... on machine learning.

Biggest weakness is ... of this conference.

We sincerely appreciate your feedback on our manuscript. Your observation regarding the limited discussion on machine learning is invaluable and has encouraged us to further explore this perspective to enhance our work.

To provide a more comprehensive overview, we have reviewed recent advancements in computational optimal transport (OT) within the machine learning community, particularly those presented at ICML. Our findings indicate a substantial body of literature [1-4] that focuses on this topic, underscoring its relevance within the field.

Furthermore, the ICML 2025 Call for Papers webpage includes "Optimization (convex and non-convex optimization, matrix/tensor methods, stochastic, online, non-smooth, composite, etc.)" as a topic of interest, highlighting the strong connection between our research and ongoing discussions in the machine learning community.

Once again, we sincerely appreciate your insights, and we would add more discussions on the importance of OT optimization in the revised paper.

References

[1] Dvurechensky, Pavel, Alexander Gasnikov, and Alexey Kroshnin. "Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn’s algorithm." ICML, 2018.

[2] Lu, Haihao, Robert Freund, and Vahab Mirrokni. "Accelerating greedy coordinate descent methods." ICML, 2018.

[3] Lin, Tianyi, Nhat Ho, and Michael Jordan. "On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms." ICML, 2019.

[4] Guminov, Sergey, et al. "On a combination of alternating minimization and Nesterov’s momentum." ICML, 2021.

In Eq. (5), ... is it needed?

Thank you for your insightful question. We cannot guarantee that the low-rank approximation is always necessary. However, based on our numerical experiments, particularly in Figures 5 and 6, we have observed that adding the low-rank term generally leads to faster convergence of the algorithm. Specifically, when the curvature condition $(s^k)^T y^k > \varepsilon \|y^k\|^2$ holds, we consistently find it beneficial to include the low-rank term.

In Eq. (6) ... behind the method.

Thank you for pointing this out! The correct notation should be $u = y_k$, and we will fix it in the main text.

As for the variables $(a, b, u, v)$ in Eq. (6) and (7), they are mostly notational: $u$ and $v$ represent two low-rank components, derived from the BFGS quasi-Newton method, and $a$ and $b$ are their respective scaling coefficients, determined by the secant equation (Eq. 6).

These terms collectively help ensure a more stable and efficient update, improving the convergence behavior of the method. We will expand the explanation in the main text to better highlight the motivation behind these choices.

The experimental results ... in the first place.

Thank you for your insightful observation. Indeed, our results confirm that the SPLR method is as good as alternatives in virtually all scenarios, and significantly outperforms others when the Hessian is sparse. This demonstrates the robustness of SPLR in worst cases and its high efficiency in favorable cases.

This behavior is due to the sparse-plus-low-rank design, which ensures stable and efficient performance across different Hessian structures.

Figure 4 seems to ... shows similar patterns.

We appreciate your observation about the impact of the regularization weight $\eta$ on convergence speed. Indeed, as shown in Figures 4 and 11, the choice of $\eta$ significantly affects all methods, including ours. While our paper focuses on optimizing the OT problem given a fixed $\eta$, we acknowledge that selecting $\eta$ is critical in practice—though this lies outside our current scope. Investigating systematic strategies for choosing $\eta$ efficiently could be an interesting direction for future work.

Another concern is ... elaborate on this?

Thank you for your valuable observation. We appreciate your concern regarding the oscillation in the Log10 Gradient Norm observed in SPLR. However, we would like to clarify that such oscillations do not pose a risk of early or incorrect stopping. In fact, the gradient norm serves as a rigorous stopping criterion. Specifically, when the gradient norm is reduced below a predefined threshold, we can theoretically prove that the optimality gap is also bounded by a corresponding threshold. This ensures that the stopping condition is reliable and provides confidence in the convergence of the method.

On the other hand, our method guarantees that the dual objective function value is nonincreasing across iterations, so if the stopping criterion is on the function value, it will not have any oscillation.

Thank you for your detailed review and thoughtful comments. Please see our point-by-point responpses below.

---

Weakness 1/Comment 1

We thank you for raising this important point about our complexity analysis. Let us clarify the theoretical aspects of our algorithm's computational efficiency:

1. Theoretical per-iteration cost: As noted in Section 5 (Lines 267-270), each iteration's dominant cost is the inversion $(H_{\Omega}^{k+1} + \tau_{k+1} I)^{-1}$ in the computation of $B_{k+1}^{-1}$. According to Chapter 10 in (Shewchuk, J. R., 1994), when it is solved using conjugate gradient method, the time complexity is $O(\rho n m \sqrt{\kappa})$, where $\kappa$ is the condition number of $H_{\Omega}^{k+1}$.
2. Convergence Guarantees:
   • Theorem 5.1 establishes linear convergence
   • The convergence constant depends on $L$ and $U$

You are absolutely correct that the sparse pattern selection affects complexity, and we appreciate the opportunity to clarify this relationship. Our method provides explicit control over this through $\rho$ while maintaining theoretical convergence guarantees.

Weakness 2/Comment 2

We sincerely thank you for identifying this important gap between our theoretical and empirical convergence results. The observed superlinear-like behavior indeed requires deeper investigation, and we offer the following possible causes:

• Approximation quality improvement: As optimization progresses (especially near the solution), our SPLR approximation appears to capture increasingly accurate Hessian information. This is evidenced in Fig. 2 (columns 2-3) where the log gradient norm decreases sharply after certain iterations.

• Hybrid behavior: The method may initially follow linear convergence (as proved) but transition to superlinear once the approximation error becomes sufficiently small.

We greatly appreciate your suggestion to pursue this theoretically, as it could lead to significant new insights about Hessian approximation methods. The gap between our current theory and observations points to exciting future research directions.

Comment 3

We thank you for this valuable suggestion. Below is a proposed comparison table we will include in the revised manuscript to clearly demonstrate the advantages of our SPLR method:

Your suggestion significantly improves our paper's ability to communicate the method's advantages, and we appreciate this constructive feedback.

Comment 4

Thanks for raising this question. The basic idea is as follows: since our method guarantees that the objective function value is nonincreasing, the iterates $\{x_k\}$ are restricted to the level set $D = \{x: f(x) \leq f(x_0)\}$, which can be proved to be compact. Since $H(x)$ is continuous on $x$, it must be bounded on $D$. Of course, the bounds may depend on the initial value $x_0$.

Comment 5

We thank you for this insightful technical question. Our rank-2 approximation is motivated by the rank-2 modification scheme in the BFGS quasi-Newton method, and can indeed be naturally generalized to higher ranks with all theoretical gurantees. The question highlights an important implementation flexibility that deserves more explicit treatment and we will adjust our manualscript accordingly.

Comment 6

We thank you for pointing out this typo. It should be $auu^T + bvv^T$.

Thank you for your detailed review and thoughtful comments. We give our point-by-point responses below.

---

Weakness 1

Thanks for raising this question on motivation. Our approach is motivated by two key insights:

1. Problem-specific Sparsity: The entropic-regularized optimal transport (OT) problem naturally exhibits (approximatly) sparse Hessian structure, particularly when $\eta$ is small. This sparsity is unique to OT and related problems.
2. General Low-Rank Utility: Low-rank approximations are successfully applied in quasi-Newton methods (e.g., BFGS, L-BFGS) because they efficiently capture dominant curvature information while maintaining computational tractability.

Our intuition is that the sparse-plus-low-rank (SPLR) structure may better preserve the problem's intrinsic geometry compared to pure sparse or pure low-rank approaches, leading to a faster convergence speed, as is shown in our empirical results. For example, if the Hessian "residuals" (i.e., the elements removed during sparsification) have some low-rank structure, then SPLR naturally captures this information.

Weakness 2

We thank you for highlighting this important consideration. Our primary limitation is that although superlinear-like convergence performance is observed, we currently lack formal theoretical justification for this observation. We will clarify such limitation in the revised manuscript.

Note 1

We thank you for this helpful suggestion to improve clarity. In the introduction, we will add a brief explanation of the cost matrix in OT problems:

The cost matrix $C \in \mathbb{R}^{n \times m}$ encodes the pairwise transportation costs between source and target distributions, where $C_{ij}$ represents the cost of moving one unit of mass from source location $i$ to target location $j$. For example, when comparing two discrete probability distributions over spatial positions, $C_{ij}$ might be the Euclidean distance $\| x_i - y_j \|^2$ between points $x_i$ and $y_j$ in the source and target domains, respectively.

We agree that further clarification will help readers, particularly those less familiar with OT, and will incorporate it into the revised manuscript.

Note 2

We thank you for this observation, and appreciate the opportunity to clarify our notation system:

Indeed, we already use iteration number as subscripts:

• $H_k$: the true Hessian matrix at iteration $k$;
• $B_k$: the approximated Hessian matrix with SPLR at iteration $k$.

The point is that we use $H_{\Omega}^{k+1}$ to distinguish between the true Hessian and the sparsified ones. We acknowledge that better notations might be possible, but this is the clearest expression we have developed so far. When we need to express the power of $H$ (which only appears in Assumption 3.2), we have added a pair of parentheses around $H$. We will add a remark in the revised paper to avoid potential ambiguity.

Question 1

We thank you for this excellent question. Here are our comments:

1. Motivation: See our response to Weakness 1 above.
2. Connection to robust PCA: The reviewer astutely observes the conceptual parallel to robust PCA. Indeed, both approaches rely on a sparse + low-rank decomposition to provide a structured approximation of a given matrix. However, while robust PCA—formulated as a Principal Component Pursuit (PCP) problem [1]—aims to recover the low-rank and sparse components through optimization, SPLR adopts a more constructive approach, providing these components via closed-form solutions. Given that robust PCA has been shown to effectively capture the essential structure of a matrix, we argue that SPLR’s approximation should likewise preserve the key information in the Hessian, ensuring a well-justified representation. Importantly, this intuition is further supported by our experimental results, where SPLR shows a superlinear-like convergence speed.

Your comment has helped us recognize that this connection deserves more explicit treatment, and we appreciate the opportunity to elaborate on these important points.

Question 2

We thank you for raising this important point. As noted, the low-rank term (Eq. 7) is computed from the sparsified Hessian $H_{\Omega}^{k+1}$​, meaning the random selection in Algorithm 2 directly influences the low-rank approximation.

The rationale behind combining Algorithm 2 with the low-rank term is to strike a balance between computational efficiency and information retention. While Algorithm 2 deliberately sparsifies the Hessian to reduce computational cost—inevitably losing some structural information—the low-rank term serves as a compensatory mechanism. This hybrid approach ensures that we maintain a meaningful approximation to the Hessian while remaining computationally tractable.

[1] Candès, Emmanuel J., et al. "Robust principal component analysis?." Journal of the ACM (JACM) 58.3 (2011): 1-37.

Thank you for your detailed review and thoughtful comments. We give our point-by-point responses below.

---

Thm.3.3 and Cor.3.4

1. We respectfully clarify that positive definiteness cannot be directly guaranteed via the Gershgorin circle theorem in this case. Specifically, for the sparsified Hessian $H_{\Omega^*}$, the Gershgorin disc centered at $h_{n+1,n+1}$ has radius $\sum_{j \neq n+1} |h_{n+1,j}|$, admitting $\lambda = 0$ as a possibility. Thus, the theorem alone does not preclude singularity. Moreover, you are correct that explicitly highlighting the improvement in the condition number would strengthen the impact of our results and we will emphasize this point more clearly.

2. Thanks for this insightful suggestion. Indeed we have conducted numerical experiments on the condition number (this was actually done before we develop the theory), and we will add such experiments in the revised paper.

3. There are two key reasons:

   1. The assumption is analogous to the concept of irreducibility in Markov chains. Irreducibility is considered a mild condition because it only requires connectivity, not quantitative bounds on probabilities and it is necessary but not sufficient for stronger properties (e.g., ergodicity).

   2. Assum. 3.2 can be satisfied with extremely sparse matrices—e.g., retaining just one complete row and column (yielding a density lower bound of only $O(1/(nm))$).

   We will add a brief remark in the revision to make this intuition clearer.

Thm. 5.1, Cor. 5.2, Thm. 5.3

Indeed the current theoretical analysis still follows the classical framework of quasi-Newton methods, but we emphasize that one key ingredient of proving the linear convergence is to bound the condition number of the approximate Hessian matrix, and in our case it is highly related to Thm.3.3 and Cor.3.4, which are unique to OT problems. Your invaluable comments have encouraged us to actively explore novel technical tools to further optimize the linear rate constant.

We will revise the paper to acknowledge the significance of the theory in Carlier (2022).

Weakness 1

See our responses above.

Weakness 2

You are absolutely correct that the Sinkhorn algorithm has a linear convergence. Our phrasing was indeed misleading, and we will correct it in the revision. Our intended point was that, while the convergence is theoretically linear, the rate can practically resemble sub-linear behavior when the regularization parameter $\eta$ is small ($1 - e^{-24\\|C\\|_{\infty}/\eta}$ approaches 1).

Also, Sinkhorn is indeed an alternating maximization procedure on the dual. We will revise for precision.

Weakness 3

We will update the manuscript accordingly.

Weakness 4

$\tilde{T}$ represents the first $(m-1)$ columns of matrix $T$, as defined in Notation (Line 68).

Weakness 5

While thresholding the top x% of elements would indeed allow direct density control, we emphasize that existing Hessian sparsification methods do not employ this approach, nor do they provide theoretical guarantees for arbitrary density selection. Specifically:

1. SNS: Uses an elementwise threshold $\rho$. Though we can set the Hessian to a desired density, there are no theoretical gurantees on invertibility or convergence rate.

2. SSNS: Sets a Hessian error threshold $\delta_k = \nu_0 \\|g_k\\|^{\gamma}$ at each iteration, which depends on the gradient norm for theoretical guarantees. This scheme cannot be set to an arbitrary density a priori, as $\delta_k$ varies with optimization progress.

In contrast, our proposed scheme explicitly enables any desired density level.

Weakness 6

The sparsification pattern $\Omega^*$ consists of all coordinates in the first row and the first column. For example, if $n=4$ and $m=4$, then $\Omega^*=\\{(1, 1), (1, 2), (1, 3), (2, 1), (3, 1), (4, 1)\\}$.

Could you please clarify in what specific way $\Omega^*$ appears ill-defined? We would be happy to provide additional explanations.

Weakness 7

We sincerely apologize for this unintentional redundancy and we will carefully revise this section to eliminate repetition.

Weakness 8

The discrepancy in our current plots likely stems from that each Sinkhorn "iteration" in our implementation includes both the $\alpha$ and $\beta$ optimization steps (two full matrix scaling operations).

Weakness 9

We thank you for this important technical observation. We will make the proof more rigorous in the revised manuscript. Below are the main ideas: we can first show that the optimal function value is finite, i.e., $f^*>-\infty$. Then the set $D^*=\\{x:f(x)\le f^*\\}=\\{x^*\\}$ is non-empty and bounded. Corollary 8.7.1 of [1] shows that $D_c=\\{x:f(x)\le c\\}$ is bounded for every $c$, which implies that $D = \\{x: f(x) \leq f(x_0)\\}$ is bounded.

[1] R Tyrrell Rockafellar. Convex Analysis. Princeton University Press, 1970.