We appreciate the recognition of the contribution of our paper. We will reply to your questions as follows.

1. Thank you for the suggestion of studying the lower bound for second-order stationary point. We can reach some conclusions based existing works. Since the hard instance constructed by Zhang et al (2021) automatically satisfies Assumption 3 and second-order stationary point is a special case of first-order stationary, the lower bound of second-order stationary point should be no smaller than first-order stationary point. In addition, the complexity of perturbed gradient method matches the result of searching first-order stationary point if ignoring logarithmic term, which indicates our method is near-optimal to achieve second-order stationary point.
2. Yes, we think the flat region is probably the stage of escaping from a bad critical point. When $x$ is far away from a first-order stationary point, the algorithm will stay in the descent phase and the convergence is fast in this stage. When $x$ is approaching a first-order stationary point, perturbations and escaping phases will occur and sometimes they could be slower than the original descent phase, which can be a possible reason for the slow convergence after $4 \times 10^4$ gradient oracles.

Thank you for the effort to review our work and we will address your concerns as follows.

1. Thank you for the suggestion of discussing the relation between saddle point of $\Phi(x)$ and local minimax point. Here we will prove that under Assumption 1 and Assumption 2, a saddle point of $\Phi(x)$ is not a local minimax point. Suppose $x$ is a saddle point of $\Phi(\cdot)$, then for $\forall \delta > 0$, there exists $x’$ such that $|| x’ – x || \le \frac{\delta}{\kappa}$ and $\Phi(x’) < \Phi(x)$. Let $y = y^*(x)$ and $y’ = y^*(x’)$. By Proposition 1 in Appendix F we have $|| y’ – y || \le \delta$. As $y$ and $y’$ are maximizers for $x$ and $x’$, we know $g_{\delta} (x) = \Phi(x)$ and $g_{\delta} (x’) = \Phi(x’)$ since $|| y’ – y || \le \delta$. Therefore, we have $g_{\delta} (x’) \le g_{\delta} (x)$ according to the definition of $x’$, which implies that $(x, y)$ is not a local minimax point.
2. We summarize some fundamental challenges as follows. Firstly, the estimation of $v_t - \nabla \phi(x_t)$ is different to the original minimax optimization as we introduce the perturbation and adopt a larger stepsize in the escaping phase. It should be estimated in these new situations. Secondly, since in minimax problem $y$ also produces noise rather than just one variable $x$, the noise of gradient estimator is larger than conventional single-level optimization. Meanwhile, the noise is also influenced by the parameters of the escaping phase (such as $r$ and $\bar{D}$) as we have mentioned, which indicates that the selection of parameters is harder than the work of the LENA algorithm. Thirdly, the dependence of $\kappa$ is not included in related works of Perturbed GD and its variant. But we need to analyze the relation between the parameters of escaping phase and $\kappa$ and provide proper options.
3. In StocBiO + iNEON, a subroutine iNEON is run to compute the negative curvature. In each inner loop iteration of this subroutine, the full batch loss function value is estimated approximately via the stochastic minibatch. In theoretical analysis the batchsize is $O(\epsilon^{-2})$, which makes the estimation precise. However, it is too expensive in practice. The precision of this estimation is crucial to StocBiO + iNEON, and hence it suffers poor performance if the batchsize used for this estimation is not large enough. In our experiment, all regular batchsize is fixed as 40. In PRGDA, a large batch of $S_1 = S_2 \cdot q = 1000$ is loaded periodically for q = 25. The averaged batchsize is small, while StocBiO + iNEON needs the large batch in every iteration in the subroutine. In our experiment, the batchsize used for estimating the loss function value in iNEON is set to 40, in which condition StocBiO + iNEON fails to escape saddle point. We conduct some additional experiments to increase the batchsize to 1/10 full batch and StocBiO + iNEON still fails to escape saddle point. When we increase the batchsize to 1/4 full batch and other hyperparameters keep the same, StocBiO + iNEON escapes saddle point and achieves the ratio of distance of 0.42 after $8 \times 10^4$ gradient oracles and 0.29 after $1.2 \times 10^5$ gradient oracles when d = 100. The result is still slower than PRGDA in Figure 2 (c).

Thank you for pointing out the problems of our presentation. We will improve our presentation and address your concerns as follows.

1. In our final version, we will refer to the SREDA algorithm in the Appendix and provide more details to ensure that readers can compare our work with previous works more easily. The most significant difference between StoBiO + iNEON and our work is that StoBiO + iNEON applies a subroutine (iNEON) to compute the negative curvature when reaching a first-order stationary point, while our method adds a perturbation. A weakness of StoBiO + iNEON is that large batch of $O(\epsilon^{-2})$ is required to estimate the loss function value in every iteration in the iNEON subroutine. In contrast, our method only requires an averaged batchsize of $O(\epsilon^{-1})$. Our experimental results show that StoBiO + iNEON cannot escape saddle point effectively with small batch, under which setting our method succeeds to find second-order stationary point.
2. In our final version, we will specify the assumptions that is used for each method in Table 2 and 3. For example, in Table 3 Acc-MDA and SREDA require Assumption 1 and 2. SGDA requires Assumption 2 (more complexity) or Assumption 1 and 2 (less complexity). Our PRGDA (minimax) requires Assumption 1, 2 and 3. The assumption of Lipschitz continuous second-order derivative is common in previous works that analyze the second-order optimality, including PGD, SSRGD, Cubic-GDA, MCN and StocBiO + iNEON. Besides, without the escaping phase, algorithms such as SREDA still cannot reach second-order stationary point even Assumption 3 is hold (considering the example that GD is trapped at saddle points).
3. We summarize some fundamental challenges as follows. Firstly, the estimation of $v_t - \nabla \phi(x_t)$ is different to the original minimax optimization as we introduce the perturbation and adopt a larger stepsize in the escaping phase. It should be estimated in these new situations. Secondly, since in minimax problem $y$ also produces noise rather than just one variable $x$, the noise of gradient estimator is larger than conventional single-level optimization. Meanwhile, the noise is also influenced by the parameters of the escaping phase (such as $r$ and $\bar{D}$) as we have mentioned, which indicates that the selection of parameters is harder than the work of Pullback algorithm. Thirdly, the dependence of $\kappa$ is not included in related works of Perturbed GD and its variant. But we need to analyze the relation between the parameters of escaping phase and $\kappa$ and provide proper options.