SERENA: A Unified Stochastic Recursive Variance Reduced Gradient Framework for Riemannian Non-Convex Optimization

Thank you for your valuable comments. You are primarily concerned with the contributions of this paper compared to previous variance reduction algorithms in Euclidean space, as well as whether the assumptions made in this paper can be satisfied. Additionally, you are interested in the impact of the parameter $\beta$ on convergence and the simplification of symbols in the theorem. Below, we will address each of these points in detail.

1. Compared to previous work in Euclidean space, our contributions are as follows.

(1) Firstly, we propose the SRVRG estimator, which does not exist in Euclidean space. Secondly, we present a unified theoretical analysis for the first time. In non-convex settings, the theoretical analysis of various variance reduction algorithms shows significant differences even in Euclidean space, with some requiring the construction of Lyapunov functions.

(2) Our Corollary 5.4 indicates that SRVRG exhibits superiority in parameter selection. For instance, in finite-sum case where $n < \varepsilon^{-2}$, when the other parameter settings adhere to those specified in the Corollary, it is sufficient for $\beta$ and $b'$ to satisfy ${\beta}/{b'} = n^{-1}$ to achieve optimal complexity. Given that $\beta = m^{-1}$, it follows that $mb' = n$. For other algorithms, although some can also achieve optimal complexity, the parameters in their analyses are fixed; for example, in R-AbaSRG, $b' = m = n^{1/2}$. Our experimental results also validate the robustness of $\beta$. We will add a remark after Corollary 5.4 to clarify the following point, for example, "Corollary5.4 indicates that if $\beta$ and $b'$ satisfy $\sqrt {\frac{\beta }{{b'}}} = \max \\{ {\frac{1}{{\sqrt n }},\varepsilon } \\}$ (finite-sum case) or $\sqrt {\frac{\beta }{{b'}}} = {\varepsilon }$ (online case), our proposed algorithm can achieve optimal complexity. This highlights the advantages of our algorithm in terms of parameter selection.".

2. Due to space limitations, we omitted the explanation of the assumptions outlined below.

(1) Assumption 2.2 (1) is to avoid a bad case in Riemannian optimization. Namely, the sequence $\{x_k\}$ may converge to an optimum $x_*$, while the connecting retraction $\{R_{x_k}(\xi_k)\}$ does not converge where $x_{k+1}=R_{x_k}(\xi_k)$. This assumption is also employed in papers (Zhou et al., 2021, Han & Gao, 2021b, Kasai et al., 2018).

(2) In Assumption 2.2 (2), the boundedness of the Riemannian Hessian is unnecessary, and we will remove it. Assumptions 2.2 (1) to (3) clearly hold for compact manifolds, including sphere, (compact) Stiefel and Grassmann manifolds. For non-compact manifolds, like SPD matrices, the assumptions hold by choosing a sufficiently small neighbourhood $\mathcal{X}$.

(3) Assumption 2.2 (4) generalizes the notion of smoothness on the euclidean space to Riemannian manifold and is satisfied if $\frac{d^2 f(R_x(t \xi))}{d t^2} \leq L$ for all $x \in \mathcal{X}, \xi \in T_x \mathcal{M}$ with $\||\xi\||=1$. Assumption 2.2 (5) and (6) are necessary due to our reliance on vector transport to approximate parallel transport. These two assumptions can be derived by requiring the vector transport $\mathcal{T}$ to be isometric and satisfy $\|| \mathcal{T}\_{x}^y u- \mathrm{D}R_x(\xi)[u] \|| \leq c_0 \|| \xi \|| \|| u \||$, where $\mathrm{D}R_x(\xi)[u]=\frac{d}{d t} R_x(\xi+t u)|_{t=0}$ is the differentiated retraction. Assumption 2.2 (6) is ensured by Taylor approximation in a compact set $\mathcal{X}$.

(4) Assumption 2.3 (1) ensures that the Riemannian distance can be expressed in terms of the inverse exponential map. Assumption 2.3 (2) hence relates the exponential map with the retraction. Indeed, we have $\left\|R_x^{-1}(y)\right\| \leq \mu\left\|\operatorname{Exp}_x^{-1}(y)\right\|$ and $\left\|\operatorname{Exp}_x^{-1}(y)\right\| \leq \nu\left\|R_x^{-1}(y)\right\|$. This assumption remains valid when $\mathcal{X}$ is sufficiently small. These assumptions are also standard as in (Sato et al., 2019, Han & Gao, 2021b)

3. The experimental results in the appendix C.3 indicate that the algorithm is robust to $\beta$. Theoretically, Corollary 5.4 indicates that setting $\sqrt{\beta/b'} = \varepsilon$ can ensure the optimal complexity of the algorithm, which means that $\beta$ and $b'$ can balance each other, and when a larger $b'$ is chosen, the algorithm will be more robust to $\beta$.

4. We will simplify the symbols in the theorem, such as $\mathcal{L} = {{\tilde L}^2}{\mu ^2}{\nu ^2}$.

We sincerely appreciate your valuable comments.

1. We will increase the discussion on the motivation behind the proposed SRVRG algorithm. The motivation for proposing SRVRG is that when $u_k$ is an SVRG-type estimator, parameter $\beta$ can be larger, which allows the variance of $v_k$ to decrease more rapidly (see the formula in the right column on lines 168-170 of the paper). Our theory also demonstrates that the range of $\beta$ is relatively expansive. Corollary 5.4 indicates that setting $\sqrt{\beta/b'} = \varepsilon$ or $\sqrt{\beta/b'} = n^{-1}$ can ensure the optimal complexity of R-SRVRG, which means that $\beta$ and $b'$ can balance each other. In contrast, Corollary 5.10 indicates that the R-SRM algorithm requires $\beta = K^{-2/3} = \mathcal{O}(\varepsilon^{2})$ to achieve optimal complexity. This is the reason why our algorithms perform better. Experimental results also validate the robustness of $\beta$.

2. Although the updates in our proposed algorithm are somewhat more complex, the average time for a single IFO call is 0.634s, which is 3.65 times the average time for a Riemannian update and 4.7 times the average time for a vector transport. This indicates that a significant portion of the algorithm's runtime is consumed by IFO calls, especially when a larger batch size $b'$ is selected in the inner loop. Our algorithm allows for a trade-off between the inner loop size $m$ and the batch size $b'$ (because in Corollary 5.4, $m = \beta^{-1}$), meaning that for the same $m$, our algorithm can support smaller batch sizes. Therefore, even when the x-axis represents running time, our algorithm still performs well. The comparison of algorithm performance with time on the x-axis demonstrates that our algorithm is comparable to the best-performing algorithms. See https://anonymous.4open.science/r/SERENA-RiemanOptimization/.

3. For the parameters of the other algorithms, we adopted most of the settings from (Han&Gao,2021b). Additionally, we performed parameter optimization for certain algorithms, such as R-PAGE. The reason that algorithms R-SRVRG and R-SVRRM, as well as R-SRM and R-Hybrid-SGD, are selected with the same parameters is that algorithm R-SRVRG can cover R-SVRRM, and similarly, algorithm R-Hybrid-SGD can cover R-SRM.

4. Thank you very much for other comments or suggestions. The final part of the definition of the exponential mapping will be modified to $\dot{\gamma}(0)=\frac{d}{d t} \gamma(t)|_{t=0}=u$. The first condition of vector transport will be modified to ''1) $\mathcal{T}$ has an associated retraction $R$, i.e., for $x \in \mathcal{M}$ and $\iota, \xi \in \mathrm{T}\_x \mathcal{M}, \mathcal{T}\_\xi(\iota)$ is a tangent vector at $R_x(\iota)$'' with the symbol $y$ is changed to $\iota$.

We thank the reviewer for the encouraging feedback and valuable comments.

We introduce the SRVRG estimator, which is proposed for the first time in both Euclidean space and Riemannian space. Numerical experiments illustrate the superiority of the R-SRVRG algorithm. Furthermore, we present a unified framework for Riemannian stochastic variance reduction methods and provide a unified theoretical analysis by giving an upper bound on the variance of the Riemannian stochastic estimators. Our proposed algorithm can achieve optimal complexity while also offering a wide range of parameter choices. For example, when other parameter settings are as outlined in the Corollary 5.4, $\beta$ and $b'$ only need to satisfy $\sqrt {\frac{\beta }{{b'}}} = \max \{ {\frac{1}{{\sqrt n }},\varepsilon } \}$ (finite-sum case) or $\sqrt {\frac{\beta }{{b'}}} = {\varepsilon }$ (online case) for the R-SRVRG algorithm to achieve optimal complexity. Other methods often directly specify parameters in order to achieve optimal complexity, such as in the analysis of R-AbaSRG for online case, $b' = \varepsilon^{-1}$. Furthermore, experimental results demonstrate robustness of parameters.

We appreciate your thorough review and useful comments. Please find our responses below.

1. Methods And Evaluation Criteria:

Thank you for your valuable feedback regarding the applicability of our proposed methods. First, please allow me to explain why the RNGD method was not compared previously.

(1) The RNGD algorithm utilizes an approximation of the Riemannian Hessian, making it a second-order method, while this paper primarily focuses on Riemannian stochastic variance reduction methods, which are the first-order methods.

(2) This paper provide a unified framework and theoretical analysis, which cannot encompass the RNGD algorithm.

We compared our proposed R-SRVGR and R-SVRRM with RNGD on the LRMC and PCA problems, respectively. The results indicate that the performance of our algorithms is comparable to that of RNGD. In particular, when the x-axis represents running time, ours perform better in most experiments. See details in https://anonymous.4open.science/r/SERENA-RieOptimization

We will replace the last sentence in the conclusion and add references. "Furthermore, Riemannian second-order methods have garnered increasing attention, such as the RNGD[1], R-SQN-VR[2], and R-SVRC[3] algorithms. We will explore the integration of this framework with second-order information to enhance the convergence rate.".

2. Other Strengths And Weaknesses:

Although the optimal complexity of the Riemannian VR method has been established in previous work, R-SRVRG and R-SVRRM show better performance. Theoretically, we first provide a unified analytical framework, as the analyses of different algorithms previously varied significantly. Additionally, there are new theoretical insights; for instance, when the other parameter settings adhere to those specified in the Corollary 5.4 in finite-sum case, our algorithms can guarantee optimal complexity as long as $mb'= \min\\{n,\varepsilon^{-2}\\}$ is satisfied and R-AbaSRG needs $b' = m = n^{1/2}$. This indicates that our algorithms exhibits superiority in parameter selection, i.e., robustness and a broader range.

3. Other Comments Or Suggestions:

(1) We will correct the typo.

(2) Assumption 2.2(1) follows papers (Zhou et al., 2021, Han & Gao, 2021b, Kasai et al., 2018), and it is difficult to relax. In Assumption 2.2(2), the boundedness of the Riemannian Hessian is unnecessary, and we will remove it. More details for assumptions, see 2 in our rebuttal to Reviewer PqwF.

(3) We will change $\tilde x$ to $\tilde x_{s-1}$ and provide the necessary explanations as follows " $\tilde x_{s-1}$ represents the point at which the true gradient is calculated in the outer loop of SVRG algorithm". For Equation (3) and (7), we will add superscript $s$ and indicate that it refers to the $s$-th outer loop. In second line after Equation (7), it is "SG".

(4) We will add the caption for Table 1 as follows: Given $\beta$, the second and third columns represent the types of stochastic gradient estimators for $u$ and $w$, respectively, while the last column indicates the corresponding algorithms.

(5) In Algorithm 1, $|\mathcal{I}_k| = |\mathcal{J}_k| = b'$, we will add " of size $b'$ " at the end of line 7. After Algorithm 1, we will add the following note "For single-loop algorithms, such as R-Hybrid-SGD and R-SRM, as well as loopless algorithms like R-PAGE, $S=1$ in Algorithm 1. And we omit the superscript for these algorithms to simplify notation, as in Equations (8) and (10).".

(6) R-SRM, R-SVRG, and R-PAGE can be classified as special cases of Algorithm 1. These three algorithms are highlighted because they each exemplify distinct types: R-SRM has a single-loop structure with its stochastic estimate represented as a convex combination of SARAH-type and SG stochastic estimates. In contrast, R-SVRG employs a double-loop structure, while R-PAGE can be regarded as a loopless SARAH-type algorithm. We present these specific forms of the stochastic estimates to facilitate theoretical analysis.

(7) For our R-SRVRG and R-SVRRM, $S = \frac{K}{m} = \Theta (\frac{1}{\sqrt {mb'} {\varepsilon ^2}})$. The optimal complexity can be guaranteed as long as $mb'= \min\\{n,\varepsilon^{-2}\\}$ is satisfied in finite sum case, thus $S = \Theta (\max \\{ n^{-1/2} {\varepsilon ^{-2},\frac{1}{\varepsilon }} \\})$. Similarly, in the online case, $mb'= \varepsilon^{-2}$, $S = \Theta ({\varepsilon ^{ - 1}})$. For R-SVRG, $S = \Theta(\varepsilon^{-2}n^{-1/3})$ in finite-sum case and $S= \Theta(\varepsilon^{-4/3})$. We will include the choice of $S$ in corollary.

[1] Hu,J. et al. Riemannian natural gradient methods. SIAM J. Sci. Comput., 2024.

[2] Kasai, H. et al. Riemannian stochastic quasi-newton algorithm with variance reduction and its convergence analysis. In AISTATS, 2018.

[3] Zhang, D. and Davanloo Tajbakhsh, S. Riemannian stochastic variance-reduced cubic regularized newton method for submanifold optimization. J. Optim. Theory Appl., 2023.

We thank the reviewer for the encouraging feedback and valuable comments and suggestions. We will respond to each point below.

1. Methods And Evaluation Criteria:

Thank you for your insightful comments regarding the evaluation criteria and complexity analysis in our paper. As you mentioned, R-SRVRG makes $4b'$ calls to IFO computation and 1 call to Riemannian update in the inner iteration, while R-SVRRM makes $3b'$ calls to IFO computation and 1 call to Riemannian update in the inner iteration. In fact, for the R-SRVRG and R-SVRRM algorithms, the value of parameter $b'$ is closely related to $m$ (because Corollary 5.4 indicates that setting $\sqrt{\beta/b'} = \varepsilon$ or $\sqrt{\beta/b'} = n^{-1/2}$ can ensure the optimal complexity of R-SRVRG and R-SVRRM, and $\beta = m^{-1}$). When $m=n$ or $m=\varepsilon^{-2}$, $b'$ can be set to 1, the ratio $(b+ 4mb')/(m+1)= \Theta((n+\varepsilon^{-2})/m)$ between IFO calls and update calls can be a small constant $\mathcal{O}(1)$, our algorithm requires $S(1+m) = \Theta(\varepsilon^{-3})$ or $\Theta(n^{-1/2}\varepsilon^{-2})$ Riemannian updates. If we set $m = b'$, only $\Theta(\varepsilon^{-2})$ Riemannian updates are needed, but the ratio will be $\Theta(\varepsilon^{-1})$ or $\Theta(n^{1/2})$.

From an experimental perspective, we recorded the time taken for 1 IFO call, 1 Riemannian update, and 1 vector transport across different problems and datasets. The average time for a single IFO call is 0.634s, which is 3.65 times the average time for a Riemannian update and 4.7 times the average time for a vector transport. This indicates that a significant portion of the algorithm's runtime is consumed by IFO calls, especially when a larger batch size is selected in the inner loop. Thus we can choose a smaller batch size $b'$.

Furthermore, we also performed a comparison of the performance of various algorithms with the x-axis representing runtime. The results show that our proposed algorithm remains comparable to the best-performing algorithms. see https://anonymous.4open.science/r/SERENA-RiemanOptimization/.

2. Experimental Designs Or Analyses:

In the experimental design, we selected the RC problem on the SPD manifold. As you mentioned, this problem is indeed a convex optimization problem. We chose this problem because both papers (Han & Gao, 2021a, 2021, Han & Gao, 2021b, Kasai et al., 2018) tested their algorithms on this issue; thus, to facilitate comparisons between algorithms, we opted for the same problem.

3. Other Comments Or Suggestions:

(1) Thank you very much for your suggestion. We will reduce the discussion on manifolds and exponential mappings in subsequent versions.

(2) We apologize for any confusion caused. Due to space limitations, we provided only a simplified definition of IFO complexity in the submission. We will include a complete version and further explanations in subsequent papers. The definition "For problem (1), an Incremental first-order oracle (IFO) takes an index $i \in \\{1,\ldots,n\\}$ and a point $x$, and returns the pair $\left(f_i(x), \operatorname{grad} f_i(x) \in \mathrm{T}_x \mathcal{M}\right)$. " and the explanations "The IFO complexity effectively captures the overall computational cost of a first-order Riemannian algorithm, as the evaluations of the objective function and gradient typically dominate the per-iteration computations."

3. We will modify $[n]$ to $\\{1,\ldots,n\\}$.

4. The first line of the formula above Equation (2) denotes the stochastic estimator for the Stochastic Recursive Momentum (STORM) algorithm, while the second line corresponds to the stochastic estimator for the Hybrid-SGD algorithm. The Hybrid-SGD estimator can be interpreted as a convex combination of the SARAH-type estimator and the stochastic gradient estimator. When $i_k = j_k$, the Hybrid-SGD reduces to STORM. And we propose a Riemannian extension of Hybrid-SGD, called R-Hybrid-SGD, as outlined in Equation (2). We will modify the formula above Equation 2 as follows:

$v_k = \overbrace { (1 - \beta )( {v_{k - 1} - \nabla {f_{{i_k}}}({x_{k - 1}})} ) + \nabla {f_{{i_k}}}({x_k}) }^{STORM} \approx \overbrace{(1 - \beta )(\underbrace { {v_{k - 1} + \nabla {f_{{i_k}}}({x_{k }}) - \nabla {f_{{i_k}}}({x_{k - 1}})}}_{SARAH-type \\;estimator}) + \beta \nabla f\_{j_k} (x_k)}^{Hybrid-SGD}.$

5. Thank you for the reminder, we will modify it.