Adaptive Riemannian ADMM for Nonsmooth Optimization: Optimal Complexity without Smoothing

Thank you for your positive feedback on the structure and clarity of the paper.

Reply to Weakness

• In the revised version, we will remove such redundant labels to improve readability and streamline the presentation.
• The mention of “Euclidean” in line 60 was a typo — we intended to refer only to smoothing-based Riemannian ADMM methods. We have corrected this in the revised manuscript and added the reference [38].
• We will include a brief explanation of the $\tilde{\mathcal{O}}$ notation in the revised manuscript. Specifically, we will clarify that $\tilde{\mathcal{O}}(\cdot)$ hides logarithmic factors in the standard $\mathcal{O}(\cdot)$ notation.
• We have revised the formula in line 172 to correctly reflect the expression $\lambda_k - \lambda\_{k+1}$, and changed the final inequality to an equality as appropriate.
• We have revised the equation to explicitly define the symbol used in the paper:

where $\lambda\_{\max}: = \frac{\beta_0 \pi^2 }{6} \\|\mathcal{A}x_0 - y_0\\|$.

• In the revised manuscript, we have clarified the definitions of all relevant parameters: The parameter $M$ is defined in Theorem 3.1. The parameters $\alpha$ and $\beta$ are defined in Equation (5). The constant $G$ is defined in Lemma A.2. The symbol $\mathcal{D}$ in line 431 now explicitly refers back to its definition in line 142, where we define $\mathcal{D} = \max\_{x, y \in \mathcal{M}} \\|x - y\\|.$ We have revised the manuscript to ensure that all such symbols are clearly defined and consistently referenced.
• We have corrected all of them in the revised manuscript. We have also thoroughly proofread the entire manuscript to eliminate other potential typographical errors.

Reply to Questions

• Regrading the parameter $\beta_k$. We thank the reviewer for the thoughtful question about the design of the dual step size $\beta_k$. The choice of $\beta_k$ is composed of two parts:

  • The first part follows a design in [10] (refer to manuscript). Its original purpose is to ensure that the dual variable in the augmented Lagrangian method remains bounded, i.e., $\\|\lambda_k\\|$ is uniformly controlled. We incorporate this term, along with the penalty parameter and the dual stepsize, to derive a sharp bound on $\\|\lambda_{k+1} - \lambda_k\\|$, as shown in Lemma 3.1. The constant $\log^2 2$ included in this part does not play a critical role in the analysis—it is included solely to maintain consistency with the $\log(k+1)$ term in the denominator.

  • The second part is introduced to achieve optimal iteration complexity. This term arises from the analysis in the proof of Theorem 3.1, where it is chosen to balance the convergence terms optimally. In particular, the denominator $\log^2(k)$ is introduced to avoid the appearance of a $\log K$ factor in the final convergence rate — as reflected in Equation (45) of our paper. We will provide a more detailed explanation of the motivation and design of $\beta_k$ in the revised manuscript to clarify both its theoretical and practical implications.

• Regarding the parameter c_\{\tau}

  • We acknowledge that the choice of c_\{\tau} may depend on certain problem-specific constants, such as the Lipschitz constant. However, we would like to clarify that this dependence is not essential. Since $\tau_k = c\_\tau k^{-1/3}$, the condition on $c_\tau$ is essentially imposed to ensure that $\tau_k$ satisfies the required bound, i.e., $\tau_k < \min\\{ \frac{1}{C}, \frac{1}{M}, \frac{1}{16 c\_{\beta}\alpha^2 \sigma_A^2}\\}$ (Theorem 3.1). This condition can always be satisfied after a sufficiently large number of iterations without requiring $c_\tau$. That is, for any $c_\tau$, there exists an integer $k_0 > 0$ such that for all $k > k_0$, $\tau_k$ satisfies the condition. Importantly, this does not affect the final convergence guarantees—the same convergence rate still holds.
  • In our experiments, we tested several different values of c_\{\tau} and selected the one that yielded the best performance. While the results do exhibit some variation across different choices of c_\{\tau}, we observed that the algorithm is not particularly sensitive to this parameter. We will clarify this point after the main theorem and provide additional details in the experimental section.

• Regarding the experiments.
  While preparing additional experiments as suggested, we identified an implementation error in one of the experiments reported in Table 2 (comparison with ManPG in our manuscript): due to an inadvertent use of an output from another method as input to our algorithm, the reported runtime was significantly underestimated. We emphasize that this issue occurred only in this particular experiment; all other experiments were correctly implemented and remain unaffected. We have corrected the implementation and re-run this specific experiment (see Table 1 below). As suggested by the reviewer, we have replaced ManPG with a stronger baseline, AManPG, in this experiment. The updated results still confirm the competitive performance of our method, and our main findings remain unchanged. All other experiments and theoretical results are unaffected by this issue. We sincerely apologize for the oversight and appreciate the reviewer’s comments, which helped us catch and resolve this error promptly.

  • Comparison with AmanPG. We compare ARADMM with the Riemannian subgradient method (RSG) and the accelerated manifold proximal gradient method (AManPG) [1]. ARADMM uses the same settings as in our original experiment. We set the step size $\eta_k \equiv 10^{-2}$ for RSG, while the code of AManPG is provided by [1]. Table 1 reports the averaged results across 10 repeated experiments with random initializations. The termination rules are the KKT conditions with an accuracy tolerance of $10^{-8}$.

  • Comparison with ALM-based methods. We compare two ALM-based algorithms: ALMSSN [2] and ALMSRTR [3], which use a semismooth Newton method and a Riemannian trust region method, respectively, to solve the augmented Lagrangian subproblem on manifolds. The codes for these algorithms are obtained from related work. Table 2 reports the function value (denoted by “Obj”), CPU time (in seconds), and the percentage of zero entries in the output (denoted by “Spa”), where the sparsity level is defined as the proportion of entries with a magnitude less than $10^{-4}$, and the results are averaged across 10 repeated experiments with random initializations. ARADMM uses the same settings as in our SPCA experiment. The termination rules are the KKT conditions with an accuracy tolerance of $10^{-8}$.

• Regrading the choice of $c\_{\rho}$

  • We have re-run the experiments with $c\_{\rho} = 1$. The result is shown in Table 1 (See reply to 3rd Reviewer xPLN due to word limit). It is shown that the results remain consistent with those reported using $c\_{\rho} = 0.6$.
  • In our response to 3nd Reviewer xPLN, we realized that the condition on $c_\tau$ was missing its dependence on $c_\rho$. We have now corrected it to:

• We have carefully addressed and corrected all the typographical errors. We sincerely thank the reviewer for the detailed and thoughtful comments, which have significantly improved the clarity and accuracy of our manuscript.

Table 1: Comparison of RSG, AManPG, and ARADMM for solving the SPCA problem with $\mu = 0.01$

Table 2: Comparison of the ALM-type methods for solving the SPCA problem with different settings

Reference （We have cited all those papers in our revised manuscript）

[1] Huang, W., & Wei, K. (2023). An inexact Riemannian proximal gradient method. Computational Optimization and Applications, 85(1), 1-32.

[2] Zhou, Y., Bao, C., Ding, C., & Zhu, J. (2023). A semismooth Newton based augmented Lagrangian method for nonsmooth optimization on matrix manifolds. Mathematical Programming, 201(1), 1-61.

[3] Zhang, C., Xiao, R., Huang, W., & Jiang, R. (2024). Riemannian Trust Region Methods for SC 1 Minimization. Journal of Scientific Computing, 101(2), 32.

Thank you for the positive feedback and for recognizing the core innovation of our work. We have addressed each concern carefully and made corresponding revisions where appropriate, as detailed below.

Reply to Questions

• Thank you for the comment. To avoid confusion, we will revise the sentence to: “Our approach involves only one Riemannian gradient evaluation and one proximal update per iteration.”
• We have removed the sentence about the indicator function.
• We have added its definition in Section 2.1 in the revised version.

Vector Transport. Given a Riemannian manifold $\mathcal{M}$, the vector transport $\mathcal{T}_x^y$ is an operator that transports a tangent vector $v \in T_x \mathcal{M}$ to the tangent space $T_y \mathcal{M}$, i.e., $\mathcal{T}_x^y(v) \in T_y \mathcal{M}$. In this paper, we assume $\mathcal{T}_x^y$ is isometric.

• Thank you for pointing this out. The inequality was missing a norm, which may have led to confusion. It should be written as:

We will correct this in the revised version of the paper. Thank you again for carefully catching this issue.

• The second equation in (30) should use "$\leq$" instead of equality. The correct derivation is as follows:

We will correct this in the revised version. Thank you again for carefully catching this oversight.

• We thank the reviewer for the careful reading. We have added the appropriate reference for Lemma A.1, which cites Lemma 1 from [10]. Regarding the vector transport: it is not required to be an isometry for the result in Lemma A.1 to hold. In fact, Line 419 can be relaxed to the more general condition

In our analysis, we assumed an isometric vector transport to simplify the notation and presentation. We will clarify this assumption and the rationale for it in the revised manuscript.

• Consistency of Assumptions on $c_\rho$.
  • Regrading the choice of $c\_{\rho}$ in robust subspace recovery, we have re-run the experiments with $c\_{\rho} = 1$. For ARADMM, we reset $c_\rho = 1$ and $\beta_0 = 700$, and retain the original values for $\rho_0$, $c_\tau$, and $c_\beta$. The implementations of other methods follow the code provided by [28] in the manuscript, where we set the ADMM stepsize $\eta = 10^{-5}$. All methods terminate when $|F(X_{k+1}) - F(X_k)| \leq 10^{-6}$ or the number of iterations exceeds 5000. The results in Table 1 are averaged over 10 runs with random initialization. It is shown that the results remain consistent with those reported using $c\_{\rho} = 0.6$.
  • We want to emphasize that the condition $c_\rho \geq 1$ is not essential. The condition appears initially in Lemma A.2, where it was used to simplify the expression in the following inequality:

Assuming $c\_\rho \geq 1$ implies $\rho_k \geq 1$, which simplifies the bound to:

However, this assumption is not strictly necessary. In fact, we can relax it by identifying a threshold index $\tilde{k} > 0$, such that for all $k > \tilde{k}$, we have $\rho_k \geq 1$. This does not affect the final convergence guarantees—the same convergence rate still holds. This type of strategy has been employed in prior work (e.g., [1]) to mitigate the reliance on precise problem constants. This is why $c\_{\rho} = 0.6$ also works well empirically. In fact, the condition for other parameters can also be eliminated by this strategy, such as $\tau_k = c_\tau k^{-1/3}$.

• Regarding lines 408 and 420, our intention was to derive the desired bounds using conditions on $c_\tau$, but we overlooked the existence of $c_\rho$ in $\rho_k$. We will revise Equation (20) to explicitly include this dependence, i.e.,

With this correction, it is sufficient to assume $c_\rho \geq 1$ throughout the analysis. In addition, for line 418, inequality should read $\rho_{k-1} = c_\rho (k-1)^{1/3} \geq 1$, not $> 1$. We will correct this in the revision.

We will revise the theoretical discussion accordingly and clarify how the convergence still holds under more general parameter settings.

Table 1: Comparison of SOC, MADMM, RADMM and ARADMM for solving the robust subspace recovery problem

• We thank the reviewer for pointing out the typos and missing terms in the proof of Theorem B.1. We have corrected the subscript inconsistencies (e.g., $\rho_k \to \rho\_{k+1}$, $\rho\_{k-1} \to \rho_k$) and added the missing constant term $(\ell_h + \lambda\_{\max})^2$, which was inadvertently omitted in the derivation. Since this term is independent of the iteration index $k$, its omission does not affect the final complexity result. We have carefully reviewed the entire proof and corrected all related issues to improve clarity and correctness in the revised manuscript.

• Thank you for the helpful suggestions. We have addressed the symbol-related issues as follows:

  • We have now added the definition of $T\mathcal{M}$ (the tangent bundle of the manifold $\mathcal{M}$) before its first use in Definition 2.1.

  • In inequality (44), the symbols $f_*$ and $h_*$ were previously undefined. We have clarified them explicitly in Assumption 3.1 as: $f_* = \inf_x f(x) > -\infty, \quad h_* = \inf_x h(x) > -\infty.$

  • Regarding the symbol $\mathcal{D}$, it is intended to denote the same quantity throughout the paper, specifically: $\mathcal{D} = \max_{x, y \in \mathcal{M}} \\|x - y\\|,$ as introduced in line 142. We will explicitly reference this definition when $\mathcal{D}$ appears in inequality (44) to ensure consistency and avoid confusion.

• We have carefully proofread the entire manuscript and corrected all identified typos and formatting issues to improve the overall clarity and presentation. We appreciate the reviewer’s suggestion and will ensure that the revised version meets a high standard of writing quality.

Reply to Limitations

Thank you for pointing this out. We acknowledge this limitation, as also noted in the paper. Our current focus is on establishing convergence without smoothing and analyzing the iteration complexity in the nonconvex, nonsmooth setting. Incorporating structural properties such as the KL inequality is indeed a promising direction and could lead to sharper rates, which we leave for future work.

Reply to Paper Formatting Concerns

We will move the assumptions to appear before the convergence analysis section. Additionally, we will add an organization paragraph to improve the clarity and guide the reader through the structure of the paper.

[1] Lu, Z., Mei, S., & Xiao, Y. (2024). Variance-reduced first-order methods for deterministically constrained stochastic nonconvex optimization with strong convergence guarantees. arXiv preprint arXiv:2409.09906.

Thank you for the positive feedback and for recognizing the core innovation of our work. We have addressed each concern carefully and made corresponding revisions where appropriate, as detailed below.

1.Regarding the Comparison to the Euclidean Counterpart:

We thank the reviewer for this insightful question. We discuss that our method differs from Euclidean ADMM in the following two aspects:

• On one hand, when the manifold reduces to $\mathbb{R}^n$, several ADMM variants with convergence guarantees have been developed (e.g., [1, 2, 3]). These methods typically rely on using the optimality condition of the $x$-subproblem and the smoothness of $f$ to derive a bound on $\\|\lambda_{k+1} - \lambda_k\\|$. However, this strategy fails to generalize to the manifold setting due to the presence of projection operators, see Section 3.1 for details. This also explains why the original Riemannian ADMM (e.g., RADMM) lacks convergence guarantees.

• On the other hand, for general nonlinear constraints (including manifold constraints), some prior works [4,5,6] apply ADMM by penalizing the constraint. However, this approach does not exploit the underlying manifold geometry and essentially reduces to penalty-based methods, which are generally less efficient than manifold-aware algorithms.

• In contrast, our method adopts a fundamentally different strategy from Euclidean ADMM to control the difference $\\|\lambda_{k+1} - \lambda_k\\|$. We achieve this through careful design of the penalty parameter and dual step size updates, which is of independent interest. Notably, our analysis does not rely on the smoothness of $f$ to control $\\|\lambda_{k+1} - \lambda_k\\|$, which opens up a potential future direction: developing ADMM-type methods for problems where both terms in the objective are nonsmooth.

We will clarify this comparison more explicitly in the revised manuscript.

2.Regarding the statement of Theorem 3.1

We thank the reviewer for pointing out this important issue. We agree that, as originally stated, the theorem could imply that the minimum of each of the three terms occurs at a different index $k$, which would not guarantee a single iterate where all terms are simultaneously small.

It follows from (46) that:

Let us denote

Now we can use Lemma B.2 to bound $R_k$. Then we multiply $\tau_k$ into both sides and sum it up from $k=1$ to $K$.
Following the proof of Theorem B.1, we can obtain the following result:

with some constant $\Gamma$. Note that $\tau_k = c\_{\tau}k^{-1/3}$ and

Then we have that

This guarantees the existence of a single iterate $k$ at which the sum of all three quantities is small, and thus the desired convergence rate holds at that point. We will update the theorem and corresponding discussion in the revised manuscript.

3.Regarding the proof of Theorem B.1

We thank the reviewer for carefully identifying the indexing inconsistency in the proof. At the beginning of the proof, the augmented Lagrangian should be written as:

with $\rho$ taking indices consistent with the corresponding variables. This correction does not affect the validity of the proof or the final result.

4.Regarding the claim that our method is adaptive

Thank you for your comment. We agree that although our algorithm is described as adaptive, it still relies on certain problem-dependent constants.

However, we would like to clarify that this dependence is not essential. Specifically, the algorithm involves three hyperparameters—$c\_\tau$, $c\_\rho$, and $c\_\beta$—which determine the sequences $\tau_k = c\_\tau k^{-1/3}$, $\rho_k = c\_\rho k^{1/3}$, and $\frac{c\_{\beta}}{k^{1/3} \log^2(k+1)}$, respectively. For example, since $\tau_k = c\_\tau k^{-1/3}$, the condition on $c_\tau$ is essentially imposed to ensure that $\tau_k$ satisfies the required bound, i.e., $\tau_k < \min\\{ \frac{1}{C}, \frac{1}{M}, \frac{1}{16 c\_{\beta}\alpha^2 \sigma_A^2}\\}$ (Theorem 3.1). This condition can always be satisfied after a sufficiently large number of iterations without requiring $c_\tau$. That is, for any $c_\tau$, there exists an integer $k_0 > 0$ such that for all $k > k_0$, $\tau_k$ satisfies the condition. A similar argument applies to the other parameters. Importantly, this does not affect the final convergence guarantees—the same convergence rate still holds. This type of strategy has been employed in prior work (e.g., [7]) to mitigate the reliance on precise problem constants.

Therefore, our use of the term adaptive was meant to emphasize that both the penalty parameter and the dual step size evolve with the iteration count $k$, rather than being fixed. A more accurate description might be that these parameters are dynamic, rather than fully adaptive in the traditional sense. We will clarify this point in the revised manuscript.

Reference

[1] Boţ, R. I., & Nguyen, D. K. (2020). The proximal alternating direction method of multipliers in the nonconvex setting: convergence analysis and rates. Mathematics of Operations Research, 45(2), 682-712.

[2] Guo, K., Han, D. R., & Wu, T. T. (2017). Convergence of alternating direction method for minimizing sum of two nonconvex functions with linear constraints. International Journal of Computer Mathematics, 94(8), 1653-1669.

[3] Hien, L. T. K., Phan, D. N., & Gillis, N. (2022). Inertial alternating direction method of multipliers for non-convex non-smooth optimization. Computational Optimization and Applications, 83(1), 247-285.

[4] Hien, L. T. K., & Papadimitriou, D. (2024). An inertial ADMM for a class of nonconvex composite optimization with nonlinear coupling constraints. Journal of Global Optimization, 89(4), 927-948.

[5] El Bourkhissi, L., Necoara, I., & Patrinos, P. (2023, December). Linearized ADMM for nonsmooth nonconvex optimization with nonlinear equality constraints. In 2023 62nd IEEE Conference on Decision and Control (CDC) (pp. 7312-7317). IEEE.

[6] Hien, L. T. K., & Papadimitriou, D. (2024). Multiblock ADMM for nonsmooth nonconvex optimization with nonlinear coupling constraints. Optimization, 1-26.

[7] Lu, Z., Mei, S., & Xiao, Y. (2024). Variance-reduced first-order methods for deterministically constrained stochastic nonconvex optimization with strong convergence guarantees. arXiv preprint arXiv:2409.09906.

We sincerely thank the reviewer for the positive feedback on the clarity of the writing and the contributions of our work. We are glad that the theoretical and practical aspects of the proposed ADMM-based algorithm were found valuable. We have addressed each concern carefully and made corresponding revisions where appropriate, as detailed below.

• Our current version mainly focuses on summarizing existing Riemannian ADMM approaches. In the revised version, we will include a more comprehensive discussion of Euclidean-space ADMM methods for nonsmooth nonconvex problems. We will discuss how our method differs from Euclidean ADMM in the following two aspects:

  • On one hand, when the manifold reduces to $\mathbb{R}^n$, several ADMM variants with convergence guarantees have been developed (e.g., [1, 2, 3] and the mentioned reference). These methods typically rely on using the optimality condition of the $x$-subproblem and the smoothness of $f$ to derive a bound on $\\|\lambda_{k+1} - \lambda_k\\|$. However, this strategy fails to generalize to the manifold setting due to the presence of projection operators, see Section 3.1 for details. This also explains why the original Riemannian ADMM (e.g., RADMM) lacks convergence guarantees.

  • On the other hand, for general nonlinear constraints (including manifold constraints), some prior works [4,5,6] apply ADMM by penalizing the constraint. However, this approach does not exploit the underlying manifold geometry and essentially reduces to penalty-based methods, which are generally less efficient than manifold-aware algorithms (e.g., Riemannian gradient and retraction).

• We appreciate the suggestion of adding an experiment on nonnegative PCA. However, we would like to clarify that the nonsmooth term in nonnegative PCA is an indicator function of the nonnegative orthant, which is not Lipschitz continuous. This violates Assumption 3.1 in our paper, which requires the nonsmooth component $h$ to be Lipschitz. For this reason, our algorithm and theoretical analysis do not directly apply to this setting, and we are unable to include it as a valid benchmark in our experiments.

• To further demonstrate the applicability of our algorithm to broader settings, we have added a new experiment: Regularized Linear Classification over the Sphere [7,8]. Consider a classification problem with $N$ training examples $\{a_i,b_i\}_{i=1}^N$ where $a_i \in \mathbb{R}^{m \times 1}$ and $b_i \in \\{-1,1\\}$ for all $i \in [N]$. The goal is to estimate a linear classifier parameter on the sphere manifold, $x \in \mathbb{S}^{m-1} := \{x \in \mathbb{R}^m : x^\top x = 1\}$, that minimizes a smooth nonconvex loss function with an $\ell_1$-regularization:

For data generation, the true parameter $x$ is sampled from $\mathcal{N}(0, I_m)$ and projected onto $\mathbb{S}^{m-1}$. The features $\\{a_i\\}\_{i=1}^N$ are sampled independently, and the labels $b_i$ are set to 1 if $x^\top a_i + \epsilon_i > 0$, where noise $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$, and -1 otherwise. All algorithms use the identical random initialization and terminate when $|F(X_{k+1}) - F(X_k)| \leq 10^{-8}$ or after 500 iterations, with $\mu = 0.2$ fixed in (1).

For the three fixed methods, we set the penalty parameter $\rho = 150$ and the step size $\eta = 0.01$ for SOC, MADMM, and RADMM. For the two adaptive methods, we set $\beta_0 = \rho_0 = c_\beta = 100$, $c_\rho = 1$, and $c_\tau = 0.05$ for ARADMM, and use the same settings as in our SPCA experiment for OADMM. Here, the choices of parameters for all the algorithms in the experiment are consistent with the corresponding theories. Before presenting the numerical comparisons, we remind the reader that there is no convergence guarantee for SOC and MADMM.

Table 1 reports the function value (denoted by “Obj”) and CPU time (in seconds), where the results are averaged over 10 repeated experiments with random initializations. From Table 1, we can see that ARADMM and MADMM quickly decrease the objective value, whereas both ARADMM and SOC achieve a lower objective value for the outputs. Moreover, ARADMM is more advantageous than existing methods in more challenging scenarios. In short, ARADMM is more efficient in terms of CPU time and objective value for test instances.

Table 1: Comparison of the ADMM-type methods for solving problem (1) with different settings. The best and second-best results are marked in bold and italics, respectively.

[1] Boţ, R. I., & Nguyen, D. K. (2020). The proximal alternating direction method of multipliers in the nonconvex setting: convergence analysis and rates. Mathematics of Operations Research, 45(2), 682-712.

[2] Guo, K., Han, D. R., & Wu, T. T. (2017). Convergence of alternating direction method for minimizing the sum of two nonconvex functions with linear constraints. International Journal of Computer Mathematics, 94(8), 1653-1669.

[3] Hien, L. T. K., Phan, D. N., & Gillis, N. (2022). Inertial alternating direction method of multipliers for non-convex non-smooth optimization. Computational Optimization and Applications, 83(1), 247-285.

[4] Hien, L. T. K., & Papadimitriou, D. (2024). An inertial ADMM for a class of nonconvex composite optimization with nonlinear coupling constraints. Journal of Global Optimization, 89(4), 927-948.

[5] El Bourkhissi, L., Necoara, I., & Patrinos, P. (2023, December). Linearized ADMM for nonsmooth nonconvex optimization with nonlinear equality constraints. In 2023 62nd IEEE Conference on Decision and Control (CDC) (pp. 7312-7317). IEEE.

[6] Hien, L. T. K., & Papadimitriou, D. (2024). Multiblock ADMM for nonsmooth nonconvex optimization with nonlinear coupling constraints. Optimization, 1-26.

[7] Zhao, L.; Mammadov, M.; and Yearwood, J. 2010. From convex to nonconvex: a loss function analysis for binary classification. In 2010 IEEE international conference on data mining workshops, 1281–1288. IEEE

[8] Zhang, D.; and Davanloo Tajbakhsh, S. 2023. Riemannian stochastic variance-reduced cubic regularized Newton method for submanifold optimization. Journal of Optimization Theory and Applications, 196(1): 324–361.

I appreciate the reviewer's detailed and thoughtful reply to my review. Particularly the clarifications regarding nonsmooth / nonconvex ADMM algorithms. It makes me more confident about the contribution and I maintain my support for this paper's acceptance by raising my confidence.