On the ergodic convergence properties of the Peaceman-Rachford method and their applications in solving linear programming

We sincerely appreciate your positive feedback and thoughtful suggestions! To address your concerns, we provide detailed responses below.

Response to the weaknesses

1. Weakness 1: Could the authors provide some convergence theorems or results for the proposed EPR-LP (Algorithm 3)?
   Response:
   Thank you for your question. Our current theoretical results are in the settings of a fixed penalty parameter without restart. We include these heuristics in EPR-LP because they can improve the practical performance of the algorithms. In the implementation of the algorithm, we usually set a maximum number of restarts. Therefore, the convergence guarantee in Corollary 2.1 still holds under practical settings with parameter tuning and restart as a direct consequence. However, the iteration complexity now only holds for iterations between two restarts. Here, we would like to further highlight the key contributions of this paper:

• This paper presents significant progress in understanding the convergence properties of the PR method by proving the convergence of its ergodic sequence. More importantly, numerical results demonstrate the superior performance of the PR method’s ergodic sequence compared to that of the DR method. This finding resolves the long-standing puzzle of why the non-ergodic sequence of the PR method outperforms that of the DR method when the PR method does converge. In particular, we study the convergence properties for a more general semi-proximal PR algorithm, which includes the PR method as a special case. Here, we highlight the motivation and importance of the extra semi-proximal terms in the considered semi-proximal PR algorithm: for solving large-scale convex optimization problems like LP problems, appropriate selection of the proximal operators $\mathcal{T}_1$ and $\mathcal{T}_2$ are necessary for solving the subproblems efficiently.

• We have established the $O(1/k)$ ergodic iteration complexity for the semi-proximal PR method, with respect to feasibility violation, objective error, and KKT-residual based on $\varepsilon$-subdifferential. A key challenge arises from the fact that the concerned function may not be sub-differentiable at the ergodic point, even if it is sub-differentiable at all the points in the non-ergodic sequence. We leverage the $\varepsilon$-subgradient to overcome this challenge.

We also want to highlight that the preliminary numerical results shown in this paper have demonstrated the superior numerical performance of EPR-LP for solving general large-scale LPs, compared to the award-winning solver cuPDLP [Applegate et al. (2021), Applegate et al. (2023), Lu and Yang (2023)].

References

1. Applegate, David, Mateo Díaz, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan O'Donoghue, and Warren Schudy. "Practical large-scale linear programming using primal-dual hybrid gradient." Advances in Neural Information Processing Systems 34 (2021): 20243-20257.
2. Applegate, David, Oliver Hinder, Haihao Lu, and Miles Lubin. "Faster first-order primal-dual methods for linear programming using restarts and sharpness." Mathematical Programming 201, no. 1 (2023): 133-184.
3. Lu, Haihao, and Jinwen Yang. "cuPDLP.jl: A GPU implementation of restarted primal-dual hybrid gradient for linear programming in Julia." arXiv preprint arXiv:2311.12180 (2023).

We hope our response addresses your concerns. We sincerely appreciate your re-evaluation of our paper.

Response to the weaknesses

1. Weakness 1: Comparison with the works of Chen et al. (2024) and Sun et al. (2024).
   Response: Thank you for your comments. We want to emphasize that the main results of this paper are different from those of Chen et al. (2024) and Sun et al. (2024). Below, we compare the differences in detail.

   This paper presents significant progress in understanding the convergence properties of the PR method by proving the convergence of its ergodic sequence. More importantly, numerical results demonstrate the superior performance of the PR method’s ergodic sequence compared to that of the DR method. This finding resolves the long-standing puzzle of why the non-ergodic sequence of the PR method outperforms that of the DR method when the PR method does converge.

   In particular, we study the convergence properties for a more general semi-proximal PR algorithm, which includes the PR method as a special case. Here, we highlight the motivation and importance of the extra semi-proximal terms in the considered semi-proximal PR algorithm: for solving large-scale convex optimization problems like LP problems, appropriate selection of the proximal operators $\mathcal{T}_1$ and $\mathcal{T}_2$ are necessary for solving the subproblems efficiently.

   Moreover, we have established the iteration complexity guarantees of the ergodic sequence generated by the semi-proximal PR algorithm for solving the convex optimization problems. More specifically, we have established the $O(1/k)$ ergodic iteration complexity with respect to feasibility violation, objective error, and KKT-residual based on $\varepsilon$-subdifferential (See Theorem 2.2 in our manuscript).

   In contrast, the studies by Sun et al. (2024) focused on studying the accelerating pADMM algorithm with Halpern iteration. The theoretical results contained in their paper cannot provide a convergence guarantee for semi-proximal PR or PR algorithms. The work of Chen et al. (2024) focuses on the implementation of the Halpern PR method proposed by Sun et al. (2024), which incorporates the restart strategy and dynamic penalty parameter updates.

   We hope the above discussion helps to clarify the key differences between this paper and the two main references mentioned in your comments.

2. Weakness 2: Achieving non-ergodic convergence is generally more challenging than ergodic convergence.
   Response: Thank you for your comments. First, we want to mention that it is known that the non-ergodic sequence of the PR method does not converge for general convex optimization problems. This is one of the key motivations for us to study the ergodic convergence. As discussed above, we established the global ergodic convergence of the more general semi-proximal PR method in this paper.

   Here, we want to highlight the challenges in the analysis of the ergodic complexity. Compared to the non-ergodic case, the analysis is more challenging. One key challenge is that the concerned function may not be sub-differentiable at the ergodic point, even if it is sub-differentiable at all the points in the non-ergodic sequence. To overcome this challenge, we leveraged the concept of an $\varepsilon$-subgradient. More specifically, we established the ergodic $O(1/k)$ complexity of the semi-proximal PR algorithm for solving convex optimization problems in terms of feasibility violation, objective error, and KKT residual using the $\varepsilon$-subdifferential framework.
   Additionally, as detailed in Appendix C in our manuscript, we establish the $o(1/\sqrt{k})$ non-ergodic iteration complexity of the pADMM method for $\rho \in (0, 2)$, in terms of feasibility violation, objective error, and KKT residual. It is worthwhile mentioning that pADMM with $\rho = 2$ is the semi-proximal PR method.

Response to the questions

1. Question 1: Clarify the difference between the DR method and the modified PDHG method.
   Response: Thanks for your comments. The modified PDHG method is equivalent to the DR method (which is pADMM with $\rho = 1$ as described in Algorithm 1 of our manuscript) with $\mathcal{S}_1= \lambda_1(AA^\*)I - AA^{*}$. For further details, the reviewer may refer to the works of Esser et al. (2010) and Chambolle & Pock (2011).

2. Question 2: In Assumption 2.2, it seems that the coefficient $\sigma$ is missing before $B_1^\* B_1$ and $B_2^\* B_2$.
   Response: Thanks for your comments. We omitted the $\sigma$ for simplicity in Assumption 2.2. For self-adjoint positive semidefinite matrices $P$ and $Q$, and a constant $\sigma > 0$, it holds that $P + \sigma Q \succ 0$ if and only if $P + Q \succ 0$. Below, we include proof for the sake of completeness:

   • For $P + Q \succ 0 \Rightarrow P + \sigma Q \succ 0$:

     • If $\sigma \geq 1$, then $P + \sigma Q \succeq P + Q \succ 0$.
     • If $\sigma < 1$, then $P + \sigma Q \succeq \sigma(P + Q) \succ 0$.

   • For $P + Q \succ 0 \Leftarrow P + \sigma Q \succ 0$:

     • If $\sigma \leq 1$, then $P + Q \succeq P + \sigma Q \succ 0$.
     • If $\sigma > 1$, then $P + Q \succeq \frac{1}{\sigma}(P + \sigma Q) \succ 0$.

3. Question 3: Proposition 2.1 establishes that an admissible preconditioner ensures Lipschitz continuity, but Line 260 presents these properties as assumptions.
   Response: Thank you for your comments. In Proposition 2.1, we established the Lipschitz continuity of $(\mathcal{T} + \mathcal{M})^{-1}$ for $\mathcal{T}$ and $\mathcal{M}$ as defined in (6) and (7). However, the convergence results for the degenerate PPM in Theorem 2.1 hold for general $\mathcal{T}$ and $\mathcal{M}$. To improve clarity, we have updated the corresponding parts in Section 2 of the revised manuscript.

References

1. Chen, Kaihuang, Defeng Sun, Yancheng Yuan, Guojun Zhang, and Xinyuan Zhao. "HPR-LP: An implementation of an HPR method for solving linear programming." arXiv preprint arXiv:2408.12179 (2024).
2. Sun, Defeng, Yancheng Yuan, Guojun Zhang, and Xinyuan Zhao. "Accelerating preconditioned ADMM via degenerate proximal point mappings." arXiv preprint arXiv:2403.18618 (2024).
3. Monteiro, Renato DC, and Chee-Khian Sim. "Complexity of the relaxed Peaceman–Rachford splitting method for the sum of two maximal strongly monotone operators." Computational Optimization and Applications 70 (2018): 763-790.
4. Esser, Ernie, Xiaoqun Zhang, and Tony F. Chan. "A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science." SIAM Journal on Imaging Sciences 3, no. 4 (2010): 1015-1046.
5. Chambolle, Antonin, and Thomas Pock. "A first-order primal-dual algorithm for convex problems with applications to imaging." Journal of Mathematical Imaging and Vision 40 (2011): 120-145.

We hope our response addresses your concerns. We sincerely appreciate your re-evaluation of our paper.

We sincerely appreciate your positive feedback and thoughtful suggestions! To address your concerns, we provide detailed responses below.

Response to the weaknesses

1. Weakness 1: High similarity between this submission and the work by Chen et al. (2024). The unique contribution of this paper remains unclear.
   Response: Thank you for your comments. While Chen et al. (2024) focus on implementing the HPR method introduced by Sun et al. (2024) for solving linear programming (LP) problems, with enhancements like adaptive restart strategies and dynamic penalty parameter adjustments, our contributions are distinct and outlined as follows:

• This paper presents significant progress in understanding the convergence properties of the PR method by proving the convergence of its ergodic sequence. More importantly, numerical results demonstrate the superior performance of the PR method’s ergodic sequence compared to that of the DR method. This finding resolves the long-standing puzzle of why the non-ergodic sequence of the PR method outperforms that of the DR method when the PR method does converge.
  In particular, we study the convergence properties for a more general semi-proximal PR algorithm, which includes the PR method as a special case. Here, we highlight the motivation and importance of the extra semi-proximal terms in the considered semi-proximal PR algorithm: for solving large-scale convex optimization problems like LP problems, appropriate selection of the proximal operators $\mathcal{T}_1$ and $\mathcal{T}_2$ are necessary for solving the subproblems efficiently.

• We have established the $O(1/k)$ ergodic iteration complexity for the semi-proximal PR method with respect to feasibility violation, objective error, and KKT-residual based on $\varepsilon$-subdifferential. A key challenge arises from the fact that the concerned function may not be sub-differentiable at the ergodic point, even if it is sub-differentiable at all the points in the non-ergodic sequence. We leverage the $\varepsilon$-subgradient to overcome this challenge.

2. Weakness 2: Including HPR-LP in the numerical experiments as a baseline for comparison.
   Response: Thank you for your comments. Compared to the numerical results presented by Chen et al. (2024), we observe that EPR-LP is slightly slower than HPR-LP, though the difference is not significant. This performance gap is understandable, given that the HPR method is a more advanced algorithm developed only recently.

3. Weakness 3: More detailed comparison with Monteiro and Sim (2018).

   Response: Thank you for bringing this paper to our attention. Monteiro and Sim (2018) employed a hybrid proximal extra-gradient (HPE) framework to analyze the ergodic iteration complexity of the PR method, which does not prove the ergodic convergence of the PR method for solving general convex optimization problems. In contrast, we established the ergodic convergence of the semi-proximal PR method (and the PR method) based on the degenerate proximal point method (dPPM). In terms of the iteration complexity, we outline the detailed results for comparison below:

• Proposition 2.2 in our paper, the ergodic iteration complexity of dPPM:

  $$\mathcal{M}(w_a^k - \bar w_a^k) \in \mathcal{T}^{\bar \varepsilon_a^k}(\bar w_a^k), \quad 0 \leq \bar \varepsilon_a^k \leq \frac{1}{4(k+1)} \\|w^0 - w^*\\|_{\mathcal{M}}^2, \quad \forall w^* \in \mathcal{T}^{-1}(0),$$

  and

  $$\\|\bar w_a^k - w_a^k\\|_{\mathcal{M}} \leq \frac{1}{(k+1)} \\| w^0 - w^* \\|\_{\mathcal{M}};$$

• Theorem 3.4 in Monteiro and Sim (2018), the ergodic iteration complexity of HPE (Adopting the notation used in our paper, we omit $\mathcal{M}$ as $\mathcal{M} \succ 0$ in this case):

  $$(w^{k-1}_a - w^k_a) \in T^{\varepsilon_k^a}(\tilde{w}_k^a), \quad 0 \leq \varepsilon_k^a \leq \frac{3}{k} \left( \\|w^0 - w^*\\|^2 + \frac{\sigma^{\prime} \\|w^0 - w^*\\|^2 }{2(1-\sigma^{\prime})} \right),$$

  for any $w^* \in \mathcal{T}^{-1}(0)$ and

  $$\\|w^{k-1}_a - w^k_a\\| \leq \frac{2 \\|w^0 - w^*\\|}{k},$$

  where $\sigma^{\prime} < 1$ is the parameter controlling the inexactness of solving the subproblems in the HPE framework.

  In our paper, we assume the operator $\mathcal{M}$ to be positive semidefinite, which allows for a broader and more general analysis. This generalization is particularly important because it enables the analysis of the ergodic iteration complexity of the semi-proximal PR method. For certain convex optimization problems, such as general LPs, appropriately chosen semi-proximal terms in the semi-proximal PR method are essential to solve the subproblem efficiently. We have included a Remark 2.1 in the revised manuscript.

We sincerely appreciate your positive feedback and thoughtful suggestions! To address your concerns, we provide detailed responses below.

Response to the weaknesses

1. Weakness 1: The proposed EPR-LP uses restart and dynamic penalty tuning. Are these heuristics covered in the convergence analysis?
   Response: Thank you for your question. Our current theoretical results are in the settings of a fixed penalty parameter without restart. We include these heuristics in EPR-LP because they can improve the practical performance of the algorithms. In the implementation of the algorithm, we usually set a maximum number of restarts. Therefore, the convergence guarantee in Corollary 2.1 still holds under the practical settings with parameter tuning and restart as a direct consequence.

2. Weakness 2: In Line 1 of Algorithm 1, compared with Algorithm 2, the variable $z^{0}$ may not be in the effective domain of the augmented Lagrangian.
   Response: Thank you for your comments. We require $(z^0, y^0, x^0) \in {\rm dom}(f_1) \times {\rm dom}(f_2) \times \mathbb{X}$; otherwise, the augmented Lagrangian function will take the value of $+\infty$. We have revised Algorithm 1 and Algorithm 3 accordingly. Here, we want to mention that finding an initial point $(z^0, y^0, x^0) \in {\rm dom}(f_1) \times {\rm dom}(f_2) \times \mathbb{X}$ is not difficult for most optimization problems.

3. Weakness 3: Is this GADMM related to the GADMM of Davis-Yin in Table 2? Elaborate on the subtle differences between the algorithms referenced in Table 2.
   Response: Thank you for your questions and suggestions. The pADMM in Algorithm 2, with $\mathcal{T}_1 = 0$ and $\mathcal{T}_2 = 0$, corresponds to the GADMM derived from the PR method (Eckstein & Bertsekas, 1992). Below, we highlight key differences between our results and those in the existing works:

   • The iteration complexity of GADMM ($\mathcal{T}_1 = 0$ and $\mathcal{T}_2 = 0$) established by Davis & Yin (2016) focuses on feasibility violation and objective error. In contrast, our work additionally establishes the iteration complexity for the KKT residual.
   • We also establish the ergodic convergence rate for pADMM with semi-proximal terms ($\mathcal{T}_1 \succeq 0, \mathcal{T}_2 \succeq 0$), encompassing feasibility violation, objective error, and KKT residual based on $\varepsilon$-subdifferential. For solving large-scale convex optimization problems, such as general LPs, appropriately chosen semi-proximal operators $\mathcal{T}_1$ and $\mathcal{T}_2$ are necessary for solving subproblems efficiently.

   Please refer to the updated Table 2 and Remark E.1 in the revised manuscript for more details.

4. Weakness 4: Compared with other non-ergodic results in Table 2, why is the non-ergodic analysis in Appendix B only of Big-O type rather than small-o?
   Response: Thank you for your question. We have refined the analysis of the non-ergodic case and obtained an $o(1/\sqrt{k})$ iteration complexity result. Please see Appendix C in the revised manuscript for details.

5. Weakness 5: For the figures in Section 4, what is the definition of the horizontal axis?
   Response:
   Thank you for your question. In Section 4, we use performance profiles as proposed by Dolan and Moré (2002), which are cited in our paper, to compare the performance of cuPDLP and EPR-LP. Performance profiles are a classical and widely accepted tool for benchmarking optimization software, including LP solvers.

   Disregarding the $\log$ scaling on the horizontal axis, $x$ represents the performance ratio, measuring how much worse a solver performs compared to the best solver for a given problem. Specifically, for a problem $p$:

   • $t_{p,s}$: The computational cost (e.g., time, iterations) of solver $s$ on problem $p$.
   • $t_{p,\text{best}} = \min_s t_{p,s}$: The computational cost of the best-performing solver on $p$.

   The performance ratio for solver $s$ on problem $p$ is defined as:

   $$r_{p,s} = \frac{t_{p,s}}{t_{p,\text{best}}}.$$
   • $r_{p,s} = 1$: Solver $s$ is the best solver for problem $p$.
   • $r_{p,s} > 1$: Solver $s$ is slower than the best solver by a factor of $r_{p,s}$.

   The horizontal axis $x$ represents the upper bound for the performance ratio $r_{p,s}$:

   • At $x = 1$, only the best-performing solvers ($r_{p,s} = 1$) are considered.
   • At $x = 2$, solvers that are at most twice as slow as the best solver ($r_{p,s} \leq 2$) are included.
   • As $x$ increases, the comparison becomes more inclusive, eventually considering all solvers for all problems.

   For further details, please refer to the paper by Dolan and Moré (2002).