Response to Reviewer HFxg

We sincerely thank the reviewer HFxg for their thorough insightful feedback and constructive suggestions. We deeply appreciate the effort you dedicated to assessing our work. We address the main concerns raised by the reviewer HFxg point by point.

1. Correct some typos (regarding Weakness 1 and Question 1)

We thank the reviewer for pointing out the typos. We have carefully re-checked all symbols, equations, and tables throughout the paper. Below, we list all the corrections we identified and will incorporate them in the camera-ready version:

We acknowledge that some of these typos may have caused confusion and may be directly responsible for the series of concerns you raised in Question 1. We sincerely apologize for this and now proceed to clarify your questions in detail. Regrettably, the rebuttal process does not allow us to submit a revised manuscript. Nevertheless, we will make sure to fix all the aforementioned issues and include a detailed discussion in response to your concerns in the final camera-ready version.

2. The necessity of the correction term in estimating the diagonal Hessian (regarding Question 1, (1))

In summary, this correction term can reduce the instability of training caused by the outliers from zeroth-order oracle with randomness. Specifically, due to the order of preconditioner is $-1/2$ ($\alpha=1/2$), too small value in the diagonal estimated Hessian may cause numerical instability (NaN) in training process. The correction term encourages numerical alignment between gradient magnitude and adaptive curvature scaling, similar in spirit to adaptive methods like Adam [1] and Adafactor [2]. To illustrate the issue more intuitively, we offer some insights through a non-strict analysis below.

In practice, we first obtain $\ell_+$, $\ell_-$ and $\ell$ through three forward pass. We obtain three numerical finite values and then $\mathcal{I}$, which, although dependent on $u$, are treated as a fixed scalar in our experiments. Thus if we do not have the correction term, we have $\mathbb{E}\left[\mathbf{\Sigma}\left(diag(uu^\top) - \mathbf{I}\right)\right] \approx 0,$

which means that in practice, we may obtain some extremely small values close to zero, causing the numerical outliers (NaN) considering the order $-1/2$. When adding the correction term, similarly ignoring the scaler, we have $\mathbb{E}\left[\mathbf{\Sigma}\left(\mathbf{I}+diag(uu^\top)\right)\left(diag(uu^\top) - \mathbf{I}\right)\right] \approx 2\mathbf{\Sigma},$

which means that $\Sigma_{t}$ is numerical remains numerically consistent with $\Sigma_{t-1}$. Thus we obtain a more stable update of $\Sigma$, which helps ensure the stability of training in our experiments. Moreover, we control the scale of correction term via $\beta_1$, which means that we can achieve numerical stability with only a small-scale correction term.

3. Why we compute $\mathbf{\tilde\Sigma} =\left( (1-\beta_1)\mathbf{\Sigma}_{t-1}^2 + \beta_1 \cdot diag^2(\tilde{g} \circ \tilde{g}) \right)^{1/2}$ (regarding Question 1, (2) and (3))

We sincerely apologize for the typo that we forgot the order of $diag$. As you mentioned, it is common to use linear-weighted average form $\mathbf{\tilde\Sigma} = (1-\beta_1)\mathbf{\Sigma}_{t-1} + \beta_1 \cdot diag(\tilde{g} \circ \tilde{g})$. However, in zeroth-order optimization, we observe that root-sum-square average $\sqrt{A^2+B^2}$ performs better compared with $A+B$. First, $\sqrt{A^2+B^2}$ is bounded by $A+B$ as $\frac{A+B}{\sqrt{2}} \leq \sqrt{A^2+B^2} \leq A+B$, which means that the two are close in order of magnitude. Second, the square-root form acts as a soft maximum. When $A$ is close to $B$, $\sqrt{A^2+B^2}$ are within the same order of magnitude, differing only by a constant factor. When $A \gg B$ or $B \gg A$, $\sqrt{A^2 + B^2}$ behaves like a softmax operator, approximating the larger value.

As you mentioned, applying such a correction causes bias in Equation (14). However, in practice, by the softmax operator $\sqrt{A^2+B^2}$ and the choice of $\beta_1$, we find that the deviation from unbiasedness is negligible empirically. On the contrary, such an almost negligible term enhances the robustness of the estimation against outliers, which is critical for the stability of zeroth-order optimization. We believe that for theory-inspired algorithms, making empirically motivated adjustments when applying them to practical tasks is both reasonable and necessary, so we named the algorithm "practical form" in Algorithm 1.

4. The choice of $\beta_1$ in the ablation study (regarding Question 1, (4))

At first we clarify that even though setting $\beta_1 =0$, PaZO is not similar to HiZOO due to the different order of estimated diagonal Hessian (different setting of $\alpha$). One of the most important contributions of this work is to demonstrate that the ideal choice of $\alpha$ is $-1/2$ instead of $0$ for MeZO and $1/2$ for HiZOO, no matter what $\beta_1$ we select. After correcting the typo in Table 6, we state the results of ablation study as below for OPT-1.3B on SST2.

The results show that when $\beta_1$ is set too large (e.g., 1, 1e-2), the optimization diverges and results in NaNs, because the proportion of correction term is too large to introduce non-negligible bias. When $\beta_1$ is too small (e.g., 0 or approaching 0), the correction term fails to take effect, which may lead to numerical instability during training and result in NaNs. Only in a proper range (around $1e\text{-}8$ to $1e\text{-}10$) does the algorithm yield reasonable performance. This behavior confirms our discussion earlier about the necessity of correction term, as we state in "2. The necessity of correction term". The use of a carefully chosen $\beta_1$ serves as a regularization mechanism, stabilizing the training process.

5. Ablation study about $T_0$ (regarding Weakness 2 and Question 2)

We briefly explain why we reset the $\Sigma$ frquently from line 228 to line 230 in the main body of our paper. This idea has been considered in periodic error feedback resetting used in optimization with communication compression [4].

As a supplement, we an ablation study on the effect of the reset period $T_0$. In this experiment, we fine-tune the OPT-1.3B model and evaluate its performance on the SST2 dataset, while varying the value $T_0 \in\{64,128,256, 512 \}$ and $\infty$ (without resetting). The results are summarized in the table below:

The reason behind this design and observation is as follows. As the number of iterations increases, the error in the Hessian estimation obtained from the zeroth-order oracle tends to accumulate over steps. This accumulated error may grow with the number of iterations. Therefore, when no resetting mechanism is used (as in the case of $\infty$ in the table above), the performance drops. On the other hand, considering that $\Sigma$ is initialized as the identity matrix $I$, there exists a gap between $I$ and the Hessian $H^*$. When the resetting frequency is too small (e.g., 128), the information accumulated in $\Sigma_t$ is insufficient to effectively bridge this gap. The hyperparameter $T_0$ plays the role of a trade-off between error accumulation and Hessian estimation accuracy. Our experimental results also indicate that the model achieves the best performance when $T_0$ falls within a certain range. We will discuss more about the ablation study on $T_0$ in the camera-ready version.

We sincerely thank the reviewer HFxg for the constructive and insightful feedback. We greatly appreciate the opportunity to improve our work through this review process, and we would be truly grateful if the reviewer HFxg could consider a positive adjustment to the scores.

Reference

[1] Kingma D P, Ba J. Adam: A method for stochastic optimization[J]. arXiv preprint arXiv:1412.6980, 2014.

[2] Shazeer N, Stern M. Adafactor: Adaptive learning rates with sublinear memory cost[C]//International Conference on Machine Learning. PMLR, 2018: 4596-4604.

[3] Chen X, Liu S, Xu K, et al. Zo-adamm: Zeroth-order adaptive momentum method for black-box optimization[J]. Advances in neural information processing systems, 2019, 32.

[4] Xie X, Lin Z, Toh K C, et al. Loco: Low-bit communication adaptor for large-scale model training[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2025.

Thank you for the detailed response, which has resolved my previous concerns.

Just two more minor points in your response:

1. Regarding the definition of $\mathcal{I}$, I think it should be $\mathbb{E}\left[ \frac{1}{2} \underbrace{u^\top \mathbf{\Sigma}^{-\frac{1}{2}} \mathbf{H} \mathbf{\Sigma}^{-\frac{1}{2}} u}_{\mathcal{I}} \cdot \mathbf{\Sigma}\left(diag(uu^\top) - \mathbf{I}\right)\right] = diag(\mathbf{H})$?

2. the statement "... the ideal choice of $\alpha$ is $-1/2$ instead of $0$ for MeZO and $1/2$ for HiZOO" requires clarification. In fact, the choice of $\alpha$ both in your algorithm and in HiZOO are $1/2$. The confusion likely arises from the way the moving average of the diagonal Hessian (inverse) is updated. In HiZOO, it is written as $\Sigma_t^{-1} = \cdots$, which is moving average of diagonal Hessian inverse, whereas in your algorithm it is expressed as $\Sigma_t = \cdots$, which is moving average of diagonal Hessian.

For the second point, I think multiple positions in the paper need to be fixed.

Thank you for your further clarification. The main point is that $\Sigma_t$ in PaZO serves as an estimate of the diagonal Hessian, whereas $\Sigma_t$ in HiZOO corresponds to an estimate of the inverse diagonal Hessian (as described in the begining of section 3.2 of HiZOO paper and the pseudocode of HiZOO). Given the definition of the preconditioner $H^{-2\alpha}$, where $H$ denotes the Hessian matrix, I still believe it is not appropriate to describe the choice of $\alpha$ in HiZOO as $-1/2$.

I have also checked the code implementation of HiZOO, which is a little bit different from the pseudocode description:

def efficient_Hessian_perturb_parameters(self, model: nn.Module, random_seed, Hessian_matrix=None, scaling_factor=1):
    torch.manual_seed(random_seed)
    for name, param in self.named_parameters_to_optim:
        z = torch.normal(mean=0, std=1, size=param.data.size(), device=param.data.device, dtype=param.data.dtype)
        param.data = param.data + scaling_factor / torch.sqrt(Hessian_matrix[name]) * z * self.args.zo_eps
    return model


def zo_Hessian_step_update(self, model, inputs, zo_learning_rate, Hessian_smooth):

    self.named_parameters_to_optim = []
    for name, param in model.named_parameters():
        if param.requires_grad:
            self.named_parameters_to_optim.append((name, param))

    random_seed = np.random.randint(1000000000)
    with torch.no_grad():
        loss_original = self.zo_forward(model, inputs)

        # first function evaluation
        model = self.efficient_Hessian_perturb_parameters(model, random_seed, self.Hessian_matrix, scaling_factor=1)
        loss1 = self.zo_forward(model, inputs)

        # second function evaluation
        model = self.efficient_Hessian_perturb_parameters(model, random_seed, self.Hessian_matrix, scaling_factor=-2)
        loss2 = self.zo_forward(model, inputs)
                 
        model = self.efficient_Hessian_perturb_parameters(model, random_seed, self.Hessian_matrix, scaling_factor=1)
        
        torch.manual_seed(random_seed)
        for name, param in self.named_parameters_to_optim:
            z = torch.normal(mean=0, std=1, size=param.data.size(), device=param.data.device, dtype=param.data.dtype)

            Hessian_temp = self.Hessian_matrix[name] * z * z
            Hessian_estimator = (torch.abs(loss1+loss2-2 * loss_original)* Hessian_temp * Hessian_smooth /(2 * self.args.zo_eps*self.args.zo_eps))
            
            self.Hessian_matrix[name] = ((1-Hessian_smooth) * self.Hessian_matrix[name] +  Hessian_estimator)

            grad = (loss1-loss2)/(2 * self.args.zo_eps) * z / torch.sqrt(self.Hessian_matrix[name])
            param.data = param.data - zo_learning_rate * (grad + self.args.weight_decay * param.data)

        loss_out = self.zo_forward(model, inputs)
    return loss_out

`

The implementation of HiZOO is with the similar way of PaZO with the choice of $\alpha=1/2$. Specifically, they compute the moving average of diagonal Hessian instead diagonal Hessian inverse, and the grad is computed through: $(l_{+} - l_{-}) / (2\mu) * \Sigma_t^{-1/2}$. The only difference lies in:

Hessian_estimator = (torch.abs(loss1+loss2-2 * loss_original)* Hessian_temp * Hessian_smooth /(2 * self.args.zo_eps*self.args.zo_eps))

In which they didn't minus the identity matrix $I$.

Response to Reviewer DhFc

We sincerely thank the reviewer DhFc for the thoughtful and constructive feedback. We appreciate the time and effort spent in evaluating our work. We provide detailed responses to each of the raised concerns in the following part, and we will incorporate the corresponding revisions and discussions into the camera-ready version.

1. The explanation on the relation between $D_\alpha$ and $d^2$ (regarding Question 1)

Thanks for your question on our clarification. In Theorem 3.5, the key parameter is $\alpha$. We define $\beta=2\alpha-1$. Let $\lambda\_1 \geq \lambda\_2 \geq \cdots \geq \lambda\_d > 0$ denote the eigenvalues of $\mathbf{H}\^*$. Consequently, $D\_\alpha$ can be written as:

Derivation of the Lower Bound. Applying the Cauchy-Schwarz inequality yields $D_{\alpha}\geq d^2$. Equality holds if and only if $\beta=0$ (which corresponds to $\alpha=1/2$), at which point $D_{\alpha}$ attains its theoretical minimum value of $d^2$.

Strict Inequality for $\alpha\neq 1/2$. When $\alpha\neq 1/2$ (i.e.,$\beta\neq 0$), the equality condition for the Cauchy-Schwarz inequality – namely, that all eigenvalues $\lambda_i$ must be equal – is generally not satisfied (unless the model is isotropic). Consequently, $D_{\alpha}>d^2$ strictly. Moreover, the magnitude of $D_{\alpha}$ relative to $d^2$ increases under two conditions: (1) as $|\beta|$ increases (i.e., as $\alpha$ deviates further from 1/2), and (2) as the dispersion of the Hessian eigenvalues $\{\lambda_i\}$ increases.

2. The scale of $\epsilon_0$ and the guarantee (regarding Question 2)

We thank the reviewer for highlighting this crucial assumption. The quantity $\epsilon_0$, introduced in the proof of Theorem 3.8, serves to bound the expected distance between $\mathbf{\theta}_t$ and $\mathbf{\theta}^*$.

Definition and Guarantee of $\epsilon_0$ The magnitude of $\epsilon_0$ reflects the degree to which the parameters $\mathbf{\theta}_t$ (for $t \geq T$) approach the optimal solution $\mathbf{\theta}^*$ after $T$ iterations. Specifically, we assume that at iteration step $t$, the expectation satisfies:

This inequality defines a neighborhood of radius $\epsilon_0$ around $\boldsymbol{\theta}^*$ within which the algorithm iterates. The achievability of such a neighborhood size $\epsilon_0$ is guaranteed by the algorithm's prior convergence properties.

Prior Convergence Properties and the Scale of $\epsilon_0$. Established analyses of Stochastic Gradient Descent (SGD) for non-convex optimization problems (e.g., [1-3]) demonstrate that, under appropriate conditions, the iterates $\mathbf{\theta}_t$ converge globally—either with high probability or in expectation—to a neighborhood containing the global optimum $\mathbf{\theta}^*$.

$\epsilon_0$ can typically be bounded by a polynomial function of the step size $\eta$. By selecting a sufficiently large total iteration number $T$, we can ensure that $\epsilon_0$ becomes sufficiently small. This, in turn, guarantees that the higher-order term $\bar{**Err**}$ in Equation (9) is asymptotically negligible relative to the dominant error term.

3. The scale of $\beta_1$ and $\beta_2$ in practice and the relation to low precision (regarding Question 3 and Limitations 1)

Thank you for pointing out the lack of clarification about the choice of $\beta_1$ and $\beta_2$. We explain the necessity of small choices of $\beta_1$ and $\beta_2$ in practice. For $\beta_1$, we first want to clarify that the correction term $\tilde g \circ \tilde g$ is designed to reduce the instability of training caused by the outliers from zeroth-order oracle with randomness. Specifically, due to the order of preconditioner is $-1/2$ ($\alpha=1/2$), too small value in the diagonal estimated Hessian may cause numerical instability (NaN) in training process. The correction term encourages numerical alignment between gradient magnitude and adaptive curvature scaling like Adam. However, applying such a correction causes bias in Equation (14). By the choice of small $\beta_1$ to control the bias, we find that such an almost negligible term enhances the robustness of the estimation against outliers and ensure the stability of training. As the results in ablation study, when $\beta_1$ is set too large (e.g., 1, 1e-2), the optimization diverges and results in NaNs, because the proportion of correction term is too large to introduce non-negligible bias. When $\beta_1$ is too small (e.g., 0 or approaching 0), the correction term fails to take effect, which may lead to numerical instability during training and result in NaNs. Only in a proper range (around $1e\text{-}8$ to $1e\text{-}10$) does the algorithm yield reasonable performance.

For $\beta_2$, it is originated from the large variance of estimated Hessian by zeroth-order oracle. We propose a non-strict analysis of variance to clarify it further as below. From an intuitive (though informal) perspective, consider the update rule: $\Sigma_t = (1 - \beta_2)\Sigma_{t-1} + \beta_2 |\widehat{\Sigma}_t|$, where $\widehat{\Sigma}_t$ denotes the current estimate of the (diagonal) Hessian, and the initialization $\Sigma_0 = I$. Since the iterative update forms the moving average, $\beta_2$ effectively controls the sensitivity to new information and hence the variance of the sequence ${\Sigma_t}$. Specifically, due to the identity initialization, the influence of $\Sigma_0 = I$ remains dominant unless $\beta_2$ is set large. In particular, when $\beta_2$ is small, the update rule behaves approximately as:

which can be interpreted as a weighted combination of the initial value $I$ and a series of estimated observations. We consider the variance of $\Sigma_t$ above by the analysis on the second term (though not a completely strict analysis, as some correlations are ignored, we hope to provide some intuition here). The expectation of each $\widehat{\Sigma}_{t-k}$ can be directly obtained as:

Thus, we can obtain the variance of each coordinate $i$ of $\widehat \Sigma$ as:

where $\hat h = \mathcal{I}{\widehat\Sigma}\left(\text{diag}(uu^\top) - I\right)$, $A = \widehat\Sigma^{-1/2}H\widehat\Sigma^{-1/2}$. Noticing that when $u \sim \mathcal{N}(0,1)$, the coupling term $\mathbb{E}[u_ju_ku_su_t(u_i^2-1)^2]$ is not 0 with special choice of $j,k,s,t$. For example, only considering the terms where $j=s$ and $k=t$, and denoting for any $k \leq t$, $\max_i (\widehat\Sigma_k)_{ii} \leq B$, we have:

It means that some values in estimated Hessian is numerically too large to introduce the risk of obtaining outliers in training. To solve this, considering:

a small $\beta_2$ yields the controlled variance of the second term to ensure a lower-variance, more stable update in the training process. Meanwhile, small $\beta_2$ still ensures that the Hessian estimates are valid from a global perspective, because when $t \to \infty$, the first term $(1 - \beta_2)^t \to 0$ and the second term $\beta_2 \sum_{k=0}^{t-1} (1 - \beta_2)^k \to 1$. Similar hyper-parameter setting is used in another zeroth-order optimizers for fine-tuning like [4]. Our experimental results also demonstrate that under the same parameter configurations, PaZO improves the performance of fine-tuned model across tasks. We will mention the sensitivity and choice of $\beta_1$ and $\beta_2$ in the main text of camera-ready version.

Moreover, your concerns about how small hyperparameters are combined with low precision training like bfloat16 are extremely important. The experiment setting in zeroth-order optimization is a little different while low-precision formats like bfloat16 are widely adopted in first-order optimization to reduce memory usage. In both the pioneering work of applying zeroth-order optimization to large language model finetuning [5], as well as subsequent follow-up studies such as [4], experiments are by default conducted in fp32 without bfloat16 to ensure numerical stability and reliable convergence. We follow the default experimental setting for fair comparison. Specifically, zeroth-order methods already lead to a significant reduction in memory usage due to its bp free—often far more than what bfloat16 achieves in standard first-order methods. In a sense, zeroth-order optimization can be viewed as an aggressive strategy to trade off between precision and memory efficiency at the algorithmic level. When applying bfloat16 in the context of zeroth-order optimization, the choice of hyperparameters may need to be adjusted accordingly to prevent loss of numerical precision. We believe this is an interesting direction for future research, but it is out of the scope of this paper.

Due to space limitations, I will immediately post the responses to the remaining questions starting from the next stage. Thank you for your insightful comments again.

6. The standard deviations in Table 6 and improvement with LoRA and prefix (regarding Question 6)

We thank the reviewer for this important question. In Table 6, our primary objective was to demonstrate the stable operating range for the hyperparameters $\beta_1$ and $\beta_2$, rather than to emphasize performance differences. Specifically, we observed that when $\beta_1$ and $\beta_2$ fall outside the proper range, the training often becomes unstable (e.g., resulting in NaNs). Therefore, we reported the mean performance within the stable region to guide practical usage. In fact, the standard deviations in these experiments are very small (below 0.2 in accuracy), and they do not affect the conclusions regarding proper $\beta_1$ and $\beta_2$ settings. For completeness, we provide the results including standard deviations below (here we correct the typo mentioned in 3).

This result clearly demonstrates that the algorithm only runs robustly when the hyperparameters are within a suitable range. When the parameters are set too large, the algorithm tends to become unstable and may even diverge.

Moreover, we would like to clarify that our primary goal is to design a zeroth-order optimizer with theoretically ideal convergence in zeroth-order optimization and practical memory efficiency compared with first-order optimization. As you mentioned, just for more difficult settings like PEFT (LoRA, prefix tuning), such improvements seem to be minor. However, these PEFT setups naturally constrain the model capacity, making it more difficult for any optimizer to yield large performance gains. Despite this, PaZO also achieves consistent improvements over HiZOO across types of models, dataset and tuning methods. Moreover, for RoBERTa-large(350M), the improvement of PaZO with PEFT is not marginal. Nonetheless, PaZO maintains stable and robust performance across various settings, and we believe this consistency, combined with its theoretical and practical advantages, makes it a compelling choice.

In summary, we truly appreciate the reviewer’s insightful questions and suggestions. In the camera-ready version, we will incorporate the clarifications above, including the discussion on theoretical analysis, hyperparameter choices and robustness (regarding Table 6), the explanation of memory overhead (regarding Figure 1), and figures regarding validation loss curves against steps which we cannot show here due to the rebuttal policy. We will also include the standard deviations. Your valuable comments have greatly helped us improve the clarity and completeness of our work. Hope our responses and discussions can address your concerns. We would be sincerely grateful if the reviewer DhFc would consider increasing the score in light of the above clarifications and revisions.

Reference

[1] HaoChen, J. Z., Wei, C., Lee, J., & Ma, T. (2021). Shape matters: Understanding the implicit bias of the noise covariance. In Conference on learning theory.

[2] Li, Y., Ma, T., & Zhang, H. (2018). Algorithmic regularization in over-parameterized matrix sensing and neural networks with quadratic activations. In Conference on learning theory.

[3] Xiong, N., Ding, L., & Du, S. S. (2023). How over-parameterization slows down gradient descent in matrix sensing: The curses of symmetry and initialization. In The twelfth international conference on learning representations.

[4] Yanjun Zhao et al. Second-order fine-tuning without pain for llms: A hessian informed zeroth-order optimizer. In: arXiv preprint arXiv:2402.15173 (2024).

[5] Sadhika Malladi et al.Fine-tuning language models with just forward passes.In:Advances in Neural Information Processing Systems 36(2023),pp.53038–53075.

Response to Reviewer bWDk

We sincerely thank reviewer bWDk for the positive evaluation and thoughtful comments on our work. We greatly appreciate your recognition of the theoretical contribution and the careful reading of both the algorithm design and the empirical results. We find the comments valuable in helping us improve the clarity of the paper. We provide detailed responses to each of the raised concerns in the following part, and we will incorporate the corresponding revisions and discussions into the camera-ready version.

1. Report on wall-clock time v.s. accuracy (regarding Cons 1 and Question)

We sincerely thank the reviewer for proposing this important point. We fully agree that comparing the wall-clock time to reach a target accuracy would more directly demonstrate the overall training efficiency of PaZO. However, due to rebuttal policy constraints, we regret that we are unable to include any PDF figures or external links in our response, and therefore cannot directly present the accuracy vs. wall-clock time plot at this stage. As a complementary reference you mentioned, we instead report the steps and wall-clock time (in seconds) required to reach 60% accuracy on a representative task RTE in the table below:

These results support our claim that PaZO achieves better convergence speed than existing zeroth-order methods by leveraging an ideal choice of the preconditioner order. Both the required number of steps and the total training time to reach the target accuracy are smaller for PaZO, validating our theoretical insights. Moreover, although the per-step cost of PaZO is slightly higher than MeZO—as we transparently report in the paper—this is more than offset by its improved convergence rate. In particular, PaZO reduces the total wall-clock time by approximately 10% compared to MeZO, demonstrating that it is not only theoretically sound but also efficient in practical settings. We will include these additional wall-clock time comparisons and corresponding wall-clock time v.s. steps plots in the camera-ready version.

2. The components of our algorithm and our main contributions (regarding Cons 2)

We sincerely thank the reviewer for the insightful comment. As you mentioned, some components of our practical implementation, such as diagonal Hessian estimation via zeroth-order oracles and the use of moving averages, are adapted from existing techniques. However, as emphasized in the Introduction, our main focus lies in addressing three questions including: (1) Is preconditioning theoretically beneficial in zeroth-order optimization? (2) If so, what is the optimal order to obtain ideal convergence? (3) How can we translate this theoretical insight into practical algorithm design? While previous works have proposed related techniques and components, they do not provide clear theoretical answers to the above questions. In contrast, our work presents a novel analysis framework that leads to a principled justification for preconditioning in zeroth-order optimization, particularly when the preconditioner order is set to $\alpha=1/2$, and demonstrates how this insight can be effectively implemented in practice. We believe that this new perspective offers substantial theoretical and practical value, and may serve as a useful work for better understanding and further advancing the field of zeroth-order optimization.

We sincerely appreciate the reviewer’s positive evaluation of our work. Your thoughtful comments and recognition of both the theoretical contributions and practical relevance of PaZO are highly encouraging to us. We will incorporate the additional experimental details on wall-clock time above into the camera-ready version to further strengthen the presentation. We would be deeply grateful if the reviewer bWDk could consider a higher score.

Response to Reviewer Gw49

We sincerely thank the reviewer Gw49 for the constructive and thoughtful feedback. We greatly appreciate your recognition of our work and your valuable suggestions. We provide point-by-point responses to each of your comments.

1. The hyperparameters in the practical version of PaZO (regarding Weaknesses 1)

Thank you very much for pointing out this important issue and helping us improve the clarity of our presentation. First, we would like to clarify that although our algorithm shares some structural similarities with HiZOO in terms of estimating the Hessian diagonal, there is a fundamental difference in the core design. Specifically, our method is guided by theoretical analysis, which leads to the choice of the preconditioner order $\alpha=1/2$ in contrast to HiZOO’s $\alpha=-1/2$. This results in essential differences in how the optimization trajectory evolves. While the Hessian estimation procedures may appear similar, the use of the preconditioner are theoretically justified and practically distinct in PaZO.

We also sincerely thank the reviewer for pointing out that our discussion of hyperparameters could be more comprehensive. We now provide additional clarification. First, we correct a typo in Table 6 of the appendix: the second column in the first row should be 1 instead of 0. We apologize for the confusion this may have caused. Then, regarding the role of the hyperparameters, it is sensitive to $\beta_1, \beta_2$ because $\beta_1$ controls the strength of the additive correction to the Hessian diagonal estimate. This correction ensures that the scaling effect in our algorithm behaves similarly to Adam-style updates, helping maintain numerical stability during training. $\beta_2$ governs the moving average applied to the Hessian diagonal estimates. This is a common technique to stabilize estimates over time and smooth the update directions from the zeroth-order oracle.

$T_0$ determines the frequency at which the preconditioner is reset. We agree that the original manuscript lacked a dedicated ablation for this parameter. We have now conducted additional experiments to investigate the effect of varying $T_0$, and we provide the results and analysis below.

In this experiment, we fine-tune the OPT-1.3B model and evaluate its performance on the SST2 dataset, while varying the value $T_0 \in\{64, 128,256, 512 \}$ and $\infty$ (without resetting). We aim to examine how the frequency of resetting the moving average of the preconditioning matrix influences training stability and final accuracy. The results are summarized in the table below:

The reason behind this design and observation is as follows. As the number of iterations increases, the error in the Hessian estimation obtained from the zeroth-order oracle tends to accumulate over steps. This accumulated error may grow with the number of iterations. Therefore, when no resetting mechanism is used (as in the case of $\infty$ in the table above), the performance drop when overly old historical accumulated error may mislead the current update direction, ultimately harming final performance. On the other hand, considering that $\Sigma$ is initialized as the identity matrix $I$, there exists a gap between $I$ and the Hessian $H^*$. When the resetting frequency is too small (e.g., 128), the information accumulated in $\Sigma_t$ is insufficient to effectively bridge this gap. The inaccurate estimation similarly leads to worse final model performance. The hyperparameter $T_0$ plays the role of a trade-off between error accumulation and Hessian estimation accuracy. Our experimental results also indicate that the model achieves the best performance when $T_0$ falls within a certain range. We will discuss more about the ablation study on $T_0$ in the camera-ready version. Thank you again for your feedback.

2. The memory overhead comparison between PaZO and HiZOO (regarding Weaknesses 2)

Thank you for pointing out the missing comparison in Figure 1. Our statement that PaZO maintains the memory efficiency of HiZOO is meant in the literal sense: both HiZOO and PaZO incur additional memory overhead primarily due to storing the diagonal estimate of the Hessian, which is of the same size in the same model. Therefore, the peak memory consumption of the two algorithms is nearly identical. For example, when fine-tuning OPT-1.3B on MultiRC, both PaZO and HiZOO have a peak memory usage of approximately 9 GB. We will update Figure 1 in the camera-ready version to include HiZOO for completeness and clarity. Thank you again for your careful review.

3. More comparison of PaZO with another methods (regarding Weakness 3 and Question 1)

We sincerely apologize for the limited scope of our experiments due to hardware constraints. However, we believe that the existing results are still indicative of the strengths of our method. While the improvement of PaZO over HiZOO on OPT is not dramatic, it is consistent across the majority of tasks. Moreover, PaZO shows improvement over HiZOO under different tuning paradigms, including LoRA and prefix tuning. This consistency suggests that the advantage of PaZO is not limited to a specific task or setup, but rather stems from a more principled design. This aligns well with the goal of our work to investigate how to better utilize preconditioning in zeroth-order optimization, and our theoretical analysis, which motivates a carefully chosen preconditioning order. We have indeed considered conducting experiments on LLaMA3. However, due to current hardware limitations, we were unable to properly apply this model. Now we are actively working to overcome these constraints and are making every effort to complete the experiments. If time permits, we may be able to report partial results during the review process, and we will include the complete results in the camera-ready version. We sincerely thank the reviewer for raising this point and appreciate your understanding.

4. The learning curve of experiments (regarding Question 2)

We sincerely thank the reviewer for raising this important question. However, due to the rebuttal policy, we are unfortunately not allowed to include any PDF files or external links in our response, and thus we are unable to directly provide the learning curves here. However, we would like to explain that all of our experiments are conducted in fine-tuning tasks. This means that, rather than focusing solely on minimizing the final loss value, our primary objective is to improve the performance of the fine-tuned model on downstream tasks—as shown in Tables 1 and 2. For fine-tuning, a lower loss does not necessarily imply better task performance.

As an alternative way, we can briefly describe an example of the validation loss v.s. steps: considering fine-tuning OPT-1.3B on RTE, the steps by which the validation loss achieves 0.62 are roughly 8000, 6000 and 5000 of MeZO, HiZOO and PaZO, respectively. These illustrate the efficiency of PaZO compared to others to some degree. We will include detailed learning curves for the experiments in Tables 1 and 2 in the camera-ready version to better illustrate the training dynamics of different methods.

Once again, we sincerely thank the reviewer Gw49 for the thoughtful comments and constructive suggestions. We have carefully addressed each of the concerns raised, and we truly appreciate the opportunity to clarify and improve our work. We are grateful for the reviewer’s recognition of our contributions. We would be deeply grateful if the reviewer Gw49 could consider a higher score.