qNBO: quasi-Newton Meets Bilevel Optimization

The initial matrix $H_0$ plays a significant role in quasi-Newton methods. Here, a scalar multiple of the identity matrix is used. Could the authors provide specific results or theoretical analysis on how different choices of $H_0$ might affect QNBO’s convergence rate and computational efficiency?

$**Response:**$ Thank you for your thoughtful questions. As discussed in Chapter 6 of [1], the initial matrix $H_0$ is often set as a scalar multiple of the identity matrix, but there is no good general strategy for selecting this scalar. We have conducted experiments to evaluate the impact of the scalar multiple of $H_0$ on the performance of qNBO algorithms in toy and data hyper-cleaning experiments. The results indicate the following: (1) different choices of the scalar multiple for $H_0$ do not affect the convergence rate or computational efficiency of qNBO (SR1) in the toy experiment; and (2) for qNBO (BFGS), the scalar multiple of $H_0$ should not be too small in the data hyper-cleaning experiment. These findings are presented in Figures 1 and 2 of a single-page PDF (containing only figures and available at the anonymous link: https://drive.google.com/file/d/1rKXLVTyE-_iSna_8xBpSknMCSYmqMmlQ/view?usp=sharing).

In theory, for the convergence result in the quadratic case, the condition $H_0 = LI$ should be explicitly stated in the theorem’s description (see Theorem 3.4 in the revised paper). For the general case, we have added a new theorem for qNBO (BFGS) in the Appendix (Theorem D.25 of the revised paper) to justify the selection of $H_0 = LI$. The proof of this theorem builds on the local convergence result of Rodomanov and Nesterov (2021b) (see Lemma D.3 of the revised paper).

[1]Nocedal, Jorge, and Stephen J. Wright, Numerical optimization. New York, NY: Springer New York, 1999.

Given the success of stochastic adaptations in other optimization frameworks, could the authors discuss potential challenges or advantages of incorporating stochastic quasi-Newton methods into QNBO? Specifically, how might this adaptation impact convergence rates or computational efficiency in practical implementations?

$**Response:**$ Thank you for your thoughtful questions. qNBO incorporates a generic decoupling structure, wherein all search directions are linear with respect to the objective functions. This structure builds upon existing methods such as HOAG (Pedregosa, 2016), AID-BIO (Ji et al., 2021), AMIGO (Arbel & Mairal, 2022), and SOBA/SABA (Dagréou et al., 2022), enabling qNBO to extend effectively from deterministic to stochastic settings.

To elaborate, qNBO consists of three components, with the first two employing quasi-Newton recursion schemes. A straightforward stochastic adaptation involves replacing deterministic quasi-Newton recursion schemes with stochastic variants, such as K-BFGS [1] or Stochastic Block BFGS [2]. Additionally, in Part 3 of qNBO, deterministic gradients are replaced with stochastic gradients throughout the iterative process, consistent with the transition from deterministic to stochastic methods discussed by Dagréou et al. (2022).

Stochastic quasi-Newton methods could be better equipped than stochastic gradient methods to handle ill-conditioning. Incorporating these methods into qNBO offers the advantage of mitigating the adverse effects of high non-linearity and ill-conditioning in the lower-level objective function by leveraging second-order information. Another potential benefit is the improvement of the constants involved in the sublinear convergence rates of stochastic methods. However, the full potential of the stochastic quasi-Newton schemes discussed (and possibly others) is not yet known.

Furthermore, this straightforward extension may compromise convergence rates and computational efficiency under the same smoothness assumptions as in the deterministic setting. For an intuitive comparison (without using quasi-Newton methods), refer to the analysis of SOBA (Dagréou et al., 2022) and its improvement by MA-SOBA (Chen et al., 2023, Optimal Algorithms for Stochastic Bilevel Optimization under Relaxed Smoothness Conditions).

In our opinion, the main challenges in maintaining comparable convergence rates and computational efficiency in a stochastic setting include:

$**Constructing effective estimators:**$ How can we develop unbiased or biased estimators in Part 3 of qNBO by integrating techniques such as variance reduction or momentum (e.g., moving averages)? Potential candidates for these estimators include SGD, SAGA, SARAH, STORM, PAGE, and even Adam. However, analyzing this extension, particularly its convergence properties, requires significant effort. For recent progress (without using quasi-Newton methods), see SRBA (Dagréou et al., 2024) and SPABA (Chu et al., 2024).

$**Analyzing convergence rates and complexity:**$ How can we evaluate the proposed stochastic algorithms in a bilevel setting while incorporating noisy second-order (curvature) information? Addressing these difficulties may require theoretical breakthroughs that go beyond existing first-order techniques, leaving this an unresolved challenge and an interesting future extension.

[1] Goldfarb, D., Ren, Y., & Bahamou, A. (2020). Practical quasi-newton methods for training deep neural networks.

[2] Gower, R., Goldfarb, D., & Richtárik, P. (2016). Stochastic block BFGS: Squeezing more curvature out of data.

The SR1 method, despite its potential in quasi-Newton methods, is noted for its numerical instability in general functions without correction strategies.

$**Response:**$ Thank you for your concern. The main drawback of SR1 updating is that the denominator in its update formula $H_{t+1}=H_t+\frac{(s_t-H_t g_t)( s_t-H_t g_t)^T}{( s_t-H_t g_t)^T g_t}$ can vanish. In fact, even for convex quadratic functions, there may be steps where no symmetric rank-one updating formula satisfies the secant equation. Such situations can lead to numerical instabilities and even a breakdown of the SR1 algorithm (see Section 6.2 of [1]).

To address these issues, a correction strategy can be applied to improve the numerical stability of SR1. The goal of such strategies is to keep the Hessian approximations under control. For detailed discussions on this topic, see Section 4.1 of Ye et al. (2023), which employs a correction strategy introduced in [2]. However, the correction strategy relies on a Hessian-based constant, which varies during iterations of $(x, y)$ in the bilevel setting. This introduces the need for an adaptive correction strategy in qNBO (SR1) to handle the dynamic nature of the Hessian in the bilevel setting. Developing such an adaptive strategy presents a distinct challenge in the bilevel optimization and could be an interesting future extension.

[1]Nocedal, Jorge, and Stephen J. Wright, Numerical optimization. New York, NY: Springer New York, 1999.

[2]Rodomanov, Anton, and Yurii Nesterov, Greedy quasi-Newton methods with explicit superlinear convergence. SIAM Journal on Optimization (2021): 785-811.

QNBO offers a flexible choice between subroutine A and B based on performance and computational cost. However, the conditions for selecting one over the other are not rigorously outlined.

$**Response:**$ Thanks for the concern. Unfortunately, at this point, we do not have a rigorous selection strategy for choosing between the two options: $Q_k=1$ or $Q_k>1$ in the iteration of $u_{k+1}$. Theoretically, if the ULI assumption described in (Ramzi et al., 2022) holds, then $Q_k=1$ would suffice, and we can share $\\{s_t, g_t\\}_{t=0}^{T-1}$ with subroutine $\mathcal{A}$. However, this assumption is quite strong (cf. Ramzi et al., 2022).

To avoid potential incorrect inversion, we adopt $Q_k-1>0$ quasi-Newton steps in subroutine $\mathcal{B}$ for the direction $\nabla_y F(x_k, y_{k+1})$. Thanks to Reviewer 5z6s’s suggestion on employing a warm-start strategy for $u_{k+1}$ , we have conducted additional experiments incorporating this strategy. See Appendix C.5.1 (Figure 6) of the revised paper for details. The empirical results show that qNBO performs better when applying the warm-start strategy for $u$. A rigorous analysis would require a more detailed characterization of the relationships between quasi-Newton iterations and their sensitivity to variations in $(x, y)$, which would be an interesting future extension.

In Remark 3.4, the authors stated that the qNBO (SR1) algorithm lacks convergence guarantees without specific corrections used to achieve numerical stability in the general case. It was also mentioned that the performance of qNBO (SR1) on the 20News dataset was omitted due to its ineffectiveness in the experiment in Section 4.2. What are the underlying reasons for these issues? A bit more discussion on these points would be helpful.

$**Response:**$ Thank you for your concern. The SR1 update formula is given by: $H_{t+1}=H_t+\frac{(s_t-H_t g_t)( s_t-H_t g_t)^T}{( s_t-H_t g_t)^T g_t}.$ However, even for convex quadratic functions, there are instances where no symmetric rank-one updating formula satisfies the secant equation. This issue arises when $s_t-H_t g_t \neq 0$ but $( s_t-H_t g_t)^T g_t =0$ (see Section 6.2 of [1]). Such scenarios can lead to numerical instabilities and even a breakdown of the SR1 algorithm. To improve the numerical stability of the SR1 method, a correction strategy should be applied. For discussions on this topic, refer to Section 4.1 of Ye et al. (2023).

[1]Nocedal, Jorge, and Stephen J. Wright, Numerical optimization[M]. New York, NY: Springer New York, 1999.

From the first terms on the right hand sides of both (6) and (7), it seems that Theorems 3.3 and 3.6 require a lower bound assumption for $\Phi(x)$. Is this correct? Additionally, in the proof sketch in Appendix D.2, $\tilde{\nabla}\Phi$ is mentioned but not defined. It would be helpful to specify its definition.

$**Response:**$ Yes, your understanding is correct. As noted in the conclusions of the proofs of Theorems 3.3 and 3.6 (Lines 1669 and 1942 of the original version), $\Phi(x^*)$ can be replaced with any lower bound of $\Phi$. This correction has been incorporated into the revised version. To clarify, we have added the standard lower bound assumption for $\Phi(x)$ in the revised version.

Thank you for pointing out this notation issue. $\tilde{\nabla} \Phi$ represents the approximate hypergradient, defined as $\tilde{\nabla} \Phi := \nabla_{x} F(x_k, y_{k+1}) - [\nabla^2_{xy} f(x_k, y_{k+1})]^T u_{k+1}$. This definition has been included in Line 165 of the revised paper.

In the experiments, how to select the step sizes and the number of inner iterations for the proposed algorithms?

$**Response:**$ The implementations and hyperparameter settings for the proposed algorithms in all numerical experiments are detailed in Section C of the Appendix.

Based on our experience, the warm-up step size $\beta$ for the inner loop $\mathcal{A}$ is typically set to 0.1. However, for logistic regression, a step size of $0.0001/(j+1)$ has been found to be more effective. The step size $\gamma$ was tuned within ${1, 0.5, 0.1}$. Specifically, we set $\gamma = 1$ in the toy example, while $\gamma = 0.1$ was more effective in other experiments. The number of inner iterations, $T$, significantly influences the parameters $\omega$ and $\tau$, which are critical for satisfying the convergence conditions. Adjustments to $T$ are made based on the condition number of the objective function in the bilevel problem. Experiments, including ablation studies on $T$ from the data hyper-cleaning experiment (Figure 9 in the Appendix of the original version), demonstrate that smaller values of $T$ lead to greater efficiency.

In Table 1, why does qNBO (BFGS) save so much time compared to SHINE for data hyper-cleaning?

$**Response:**$ The reason is that qBNO employs a constant step size strategy instead of the time-consuming line search used in SHINE for the inner solver. To clarify, consider the data hyper-cleaning experiment on the MNIST dataset. The SHINE-OPA algorithm requires 23.19 seconds to achieve a test accuracy of 91.51%. (The 20.25 seconds reported in Table 1 of the paper is the average time across ten runs.) Notably, the line search method in the lower-level solver accounts for 19.61 seconds of the total time.

Although quasi-Newton methods are used to estimate the inverse Hessian-vector product, the proposed algorithms still require computing Jacobian-vector products, which can be computationally expensive in large-scale cases. Could this computation pose a potential problem?

$**Response:**$ Thank you for the interesting question. Apart from the computation of Jacobian-vector products, all other operations in qNBO involve inner products and vector summations. Consequently, when Jacobian-vector products are not computationally expensive, qNBO can significantly reduce computational costs.

Modern automatic differentiation frameworks, such as PyTorch and JAX, make the computational complexity of Jacobian-vector products comparable to that of gradient computation. Alternatively, recent methods in the bilevel optimization literature offer different strategies for estimating Jacobian-vector products. For instance, PZOBO (Sow et al., 2022b) computes these products using differences between two optimization-based trajectories, while FdeHBO (Yang et al., 2023) and HJFBiO (Huang, 2024) employ finite-difference approximations. Additionally, ZDSBA (Aghasi & Ghadimi, 2024) uses Gaussian smoothing as an effective alternative.

SHINE incorporates two different OPA (Outer Problem Awareness) corrections as described in Ramzi et al. (2022). $**The first correction involves the cross-derivatives of the inner-level function$f_{xy**(x, y)$when the upper-level variable$x is one dimentianal} (i.e., $f_{xy}(x, y)$ is a vector). This correction is outlined in equation (5) of Ramzi et al. (2022) and explained above it. It corresponds to the term $f_{xy}[f_{yy}]^{-1}$ in the computation of the hypergradient and requires $x$ to be one-dimensional.

$**The second OPA correction, described in equation (8) of Ramzi et al. (2022), involves$F_y(x,y)$**$, which is a vector regardless of the dimensionality of $x$. This correction corresponds to $[f_{yy}]^{-1}F_y$ in the computation of the hypergradient. $**This is the version we used in both the toy experiment and other experiments.**$

The specific file paths for verifying the implementation of SHINE are as follows:

$**SHINE:**$ From ICLR2025\qNBO\qNBO\logisticregression\hoag\hoag.py (lines 308 and 383) and ICLR2025\qNBO\qNBO \logisticregression\hoag\lbfgs.py (line 84), it can be confirmed that the code for SHINE aligns with equation (8) of Ramzi et al. (2022). Additionally, it matches the original SHINE open-source code (lines 147–148 in https://github.com/zaccharieramzi/hoag/blob/shine/hoag/hoag.py and lines 80–83 in https://github.com/zaccharieramzi/hoag/blob/shine/hoag/lbfgs.py).

(2)Since the implementations of AID-TN, AID-BIO, and AMIGO-CG in the above experiments used the closed-form expression $f_{yy} = A$ for the Hessian, $**it would not be**$ $** fair to compare them with other algorithms that do not use the explicit Hessian information.**$ Therefore, $**we also conducted the toy experiment**$ $**using the more practical torch.autograd.grad for computing Hessian-vector products.**$ Note that this approach is typical in the machine learning literature, as the closed form of the Hessian is generally difficult to compute in practical applications. After performing a grid search to fine-tune the hyperparameters for step sizes tested across [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.2, 0.3, 0.4, 0.5], TN/CG iterations evaluated at [1, 10, 20, 30, 40, 50], and inner gradient iterations between 1 and 10, the new results—using torch.autograd.grad for computing the Hessian-vector product instead of the closed Hessian form—are presented in Figure 2 of the single-page PDF at the anonymous link and in the following tables.

These results show that: (i) for AID-BIO and AMIGO-CG, 50 CG iterations yielded the best results, while a single CG iteration caused divergence, resulting in much worse performance compared to the scenarios using the closed-form Hessian; (ii) compared to other algorithms, both qNBO (BFGS) and qNBO (SR1) exhibit smaller hypergradient errors and produce iterates that are closer to the optimal solutions. $**Therefore, after implementing the debugging suggestions from Reviewer yCCW and using the more practical torch.autograd.grad for**$

$**Hessian-vector products, the results for the toy experiment align with the conclusions of this work.**$

To facilitate the reproduction of the results, the specific hyperparameter settings are as follows:

$**AID-BIO/AMIGO-CG:**$ The maximum number of outer iterations is $K = 5000$, CG iterations is $N = 50$, the inner gradient iteration is $T = 1$, the inner step size is $\beta = 0.01$, and the outer step size is $\alpha = 0.01$.

$**AID-TN:**$ The maximum number of outer iterations is $K = 5000$, TN iterations is $N = 50$, the inner gradient iteration is $T = 1$, the inner step size is $\beta = 0.01$, and the outer step size is $\alpha = 0.2$.

The hyperparameters of qNBO (SR1) are the same as those used in the above experiment, and the settings for the other algorithms remain unchanged from those in the paper.

To address these issues, we have re-conducted the toy experiment, as detailed below:

(1)$**By retaining the closed-form expression$f_{yy** = A of the Hessian in the implementations of AID-TN, AID-BIO, and AMIGO-CG}, we performed a grid search to fine-tune the hyperparameters of all algorithms, including inner loop iterations and step sizes. As noted by Reviewer yCCW, AID-BIO and AMIGO-CG achieve a relative error of $10^{-5}$ on the upper-level variable (compared to the values reported in the paper: $10^{-2}$ for AID-BIO and $10^0$ for AMIGO-CG) after setting the inner loop iterations for CG steps and gradient steps to 1 while maintaining the outer step size at 0.01.

By further conducting a grid search over TN/CG iterations in the range [1, 5, 10, 15, 20], inner gradient iterations between 1 and 10, and step sizes within [0.1, 0.2, 0.3, 0.4, 0.5], $**the relative error of AID-BIO and AMIGO-CG can be reduced to below$10^{-6**$within 1.5 seconds when the outer step size is adjusted to$\alpha = 0.2}. At the same time, after fine-tuning the hyperparameters of qNBO (SR1), $**the new implementation of qNBO (SR1) can also reduce the error to below$10^{-6** within 1.5 seconds}.

The updated results are presented in Figure 1 of a single-page PDF (containing only figures and available at the anonymous link: https://drive.google.com/file/d/15x5Fp1XYlo1DKrKs4RcquXU6khqv9Y1l/view?usp=drive_link) and in the following tables.

Note that the reason the performance of the AID-BIO and AMIGO-CG algorithms differs slightly in the updated results is due to incomplete memory release. To facilitate reproduction of the results, the specific hyperparameter settings are as follows:

$**AID-BIO/AMIGO-CG:**$ The maximum number of outer iterations is K = 5000, the CG iterations is $N = 1$, the inner gradient iteration is $T = 1$, the inner step size is $\beta= 0.01$ and $**the outer step size is$\alpha = 0.2$**$.

$**AID-TN:**$ The maximum number of outer iterations is $K = 5000$, the TN iteration is $N = 1$, the inner gradient iteration is $T = 1$, the inner step size is $\beta = 0.01$ and $**the outer step size is$\alpha = 0.2$**$.

$**qNBO (SR1):**$ The maximum number of outer iterations is $K = 5000$, the inner iterations are $T = 6$, $P = 9$, $Q_k = 25$, the inner step sizes are $\beta = 0.01$, $\gamma = 1$, the initial matrix is $H_0 = I$, and the outer step size is $\alpha = 0.5$.

We sincerely thank Reviewer yCCW for the thorough and valuable feedback on our code. Before addressing the issues identified with AID-TN, AID-BIO, and AMIGO-CG in the toy experiment, we first assessed their potential impact on other experiments.

Our analysis indicates that these issues do not affect the outcomes of other experiments. This conclusion is based on the fact that $**the toy experiment and other experiments use separate algorithm files.**$ Specifically, all the code for the algorithms in the toy experiment is contained in the file ICLR2025\qNBO\qNBO\toy, which was rewritten exclusively for the toy experiment. Notably, the algorithms AID-TN, AID-BIO, and AMIGO-CG in the toy experiment utilize a closed-form Hessian (discussed later). In contrast, their implementations in other experiments rely on code from the file ICLR2025\other algorithms\solvers, which originates from the open-source implementation cited in the paper.

The following provides the specific file paths for thoroughly checking the bugs mentioned by Reviewer yCCW for AID-TN, AID-BIO, and AMIGO-CG during the toy experiment. $**Notably, these bugs are absent from the codes in**$ ICLR2025\other algorithms\solvers $**used for other experiments.**$

$**AID-TN:**$ From ICLR2025\other algorithms\solvers\stocbio.py (lines 84–95) and ICLR2025\other algorithms\benchmark_utils\hessian_approximation.py (lines 116–121), it can be verified that the code of AID-TN used in other experiments aligns with Reviewer yCCW’s debug on AID-TN for the toy experiment.

$**AID-BIO:**$ From ICLR2025\other algorithms\solvers\stocbio.py (lines 84–95), it can be verified that the code of AID-BIO used in other experiments aligns with Reviewer yCCW’s debug on AID-BIO for the toy experiment.

$**AMIGO-CG:**$ From ICLR2025\other algorithms\solvers\amigo.py (lines 87–98 and lines 99–102), it can be verified that the code of AMIGO-CG used in other experiments aligns with Reviewer yCCW’s debug on AMIGO-CG for the toy experiment.

We greatly appreciate Reviewer yCCW’s insightful debugging of AID-TN, AID-BIO, and AMIGO-CG in our toy experiment. The suggestions are both accurate and helpful, aligning with the code used in other experiments. Taking this opportunity, we further reviewed the toy experiment and identified two additional directions for improvement in the toy experiment.

First, we need to fine-tune the hyperparameters of all algorithms (e.g., inner loop iterations and step sizes) to ensure a fair comparison. Previously, the inner loop setting (10 steps) for AID-TN, AID-BIO, and AMIGO-CG was adopted directly from other experiments without optimization.

Second, in the toy experiment, we used the closed-form expression $f_{yy} = A$ for the Hessian in the implementation of AID-TN, AID-BIO, and AMIGO-CG. Specifically, the code ICLR2025\qNBO\qNBO\toy\toy3 (lines 750, 814, and 876) calls the function f_yy1 (lines 99–100 of toy3) to directly return the Hessian matrix $A$. However, a more practical approach would involve computing Hessian-vector products using torch.autograd.grad. It is important to note that only the implementations of AID-TN, AID-BIO, and AMIGO-CG require the Hessian.

We compared the performance of AID-BIO with different numbers of CG iterations ($T=1, 10$). The results are presented in Figure 3 of the single-page PDF at the anonymous link and in the following tables. Using the number of outer-loop iterations of the algorithm on the x-axis, we observed that both CG iteration strategies exhibited nearly identical performance. This suggests that, due to the closed form of the Hessian provided in the code, a single CG iteration is sufficient to reach the optimal solution. Additional CG iterations do not improve computational accuracy, making $T=1$ more time-efficient than $T=10$.

For the theory, the results in Theorems 3.3 and 3.6 are limited to the setting where $Q_k=k+1$. As a result, qNBO is a double-loop algorithm. Is it possible to design a single-loop version? Recent progress has been made toward single-loop bilevel optimization algorithms, especially in the nonconvex-strongly-convex setting, by using a warm-start strategy.

$**Response:**$ Thank you for the suggestion. A warm-start strategy is reasonable based on the Lipschitz continuity of $u^*(x, y)$ under the given assumptions. We have incorporated experiments using the warm-start strategy you suggested for $u$. As shown in Appendix C.1 (Figure 6 on Page 17) of the revised paper, the empirical results demonstrate that qNBO(BFGS) performs better when applying the warm-start strategy for $u_k$ (denoted as $Q_k = \text{ws}$ in the Figure 6). A rigorous analysis would require a more explicit characterization of the relationships between quasi-Newton iterations and their sensitivity with respect to $(x, y)$, which would be an interesting future extension. In our opinion, this could serve as a first step toward developing a single-loop version of qNBO. Thank you for your insightful question; we will continue to explore this idea.

For the experiments, since the value of $Q_k$ affects the running time, it would be beneficial to empirically demonstrate how increasing $Q_k$ influences performance gains. Note that the numerical results in Fig. 1(d) illustrate only the impact of $Q_k$ as the outer iteration $k$ varies.

$**Response:**$ Thank you for the suggestion. In the original submitted version (Figure 10 on Page 22 in Appendix C.5), an ablation study on the iteration number $Q_k$ of qNBO (BFGS) in the data hyper-cleaning experiment was performed using running times. We have added additional ablation study experiments based on running times for the toy experiment; see Appendix C.5.1 (Page 17, Figure 6) in the revised version. As shown in Figure 6, when $Q_k$ is constant, an increase in $Q_k$ initially slows the decline in various indicators over time but ultimately leads to a smaller minimal error. However, $Q_k = k+1$ demonstrates superior performance in terms of both the rate of decline and the minimal error achieved, even when its initial number of steps is fewer than that of constant $Q_k$. Furthermore, employing a warm-start strategy for $u_k$ improves the performance of qNBO (BFGS), highlighting the potential advantages of this approach, as noted by the reviewer.

Page 6: Do Theorems 3.3 and 3.6 require the assumption that an optimal solution $x^∗$ exists? Note that $\Phi(x^*)$ appears in both (6) and (7).

$**Response:**$ Thank you for your questions. Theorems 3.3 and 3.6 do not assume the existence of optimal solutions. As noted at the conclusions of their proofs (Lines 1669 and 1942 of the original version), $\Phi(x^*)$ can be replaced with any lower bound of $\Phi$. Thus, we add the lower bound assumption on $\Phi(x)$ (Assumption 3.3 in the revised paper) and replace $\Phi(x^*)$ with $\inf_x \Phi(x)$ in Theorems 3.4 and 3.7 and the corresponding proofs in the appendix of the revised version.

Page 7: Why is it stated that qNBO is more efficient than AID-BIO? Note that the notation $\tilde{\mathcal{O}}$ omits the $\log(1/\epsilon)$ term.

$**Response:**$ Thank you for your concern. Your understanding is correct—it is inaccurate to claim that qNBO is more efficient than AID-BIO in terms of matrix-vector complexity. According to Table 1 of AID-BIO (Ji et al., 2021), the number of Jacobian-vector products for AID-BIO is of the order $\mathcal{O}(\kappa^{3}\epsilon^{-1})$, and the number of Hessian-vector products is of the order $\mathcal{O}(\kappa^{3.5}\epsilon^{-1})$. When the dimension of the lower-level variable exceeds that of the upper-level variable (e.g., the model parameter’s dimension is larger than that of the hyperparameter), Hessian-vector products become more computationally intensive. Consequently, the matrix-vector complexity for AID-BIO is $\mathcal{O}(\kappa^{3.5}\epsilon^{-1})$, compared to $\tilde{\mathcal{O}}(\kappa^{3}\epsilon^{-1})$ for qNBO. Thus, qNBO is not more efficient than AID-BIO in terms of $\epsilon$-order; it is only more efficient in the $\kappa$ term within matrix-vector complexity. This correction has been included in the revised version.

In Appendix C.5, the authors report various ablation studies on the parameters $T$ and $Q$, but they do not specify which experiment these ablation studies are based on.

$**Response:**$ The ablation studies described in Appendix C.5 are based on the data hyper-cleaning experiment. This information has been added to Appendix C.5.2 of the revised version of the paper. Thank you for bringing this to our attention.

Page 6: In Eq. (6) of Theorem 3.3, the constant Mfxy can be computed, as f takes the quadratic form given in (5).

$**Response:**$ Thank you for your suggestion. Since $\nabla^2_{xy} f(x, y) = -I$ for the lower-level objective function given in Eq. (5), the constant $M_{f_{xy}} = 1$. We have made the corresponding revisions in the paper.

Page 8: In Eq. (9), ${a_i}’$, ${b_i}’$” should be ai′, bi′”.

$**Response:**$
Thank you for pointing this out. The corrections to Eq. (9) have been made in the revised version.

Questions: Can this method be effectively applied in a stochastic setting while maintaining comparable convergence rates and computational efficiency? The required adjustments in iteration rounds and step size suggest potential difficulties in achieving the same convergence speed as existing methods.

$**Response:**$ Thank you for your thoughtful questions. This work introduces a flexible algorithmic framework, qNBO, designed to improve hypergradient approximation. qNBO incorporates a generic decoupling structure, wherein all search directions are linear with respect to the objective functions. This structure builds upon existing methods such as HOAG (Pedregosa, 2016), AID-BIO (Ji et al., 2021), AMIGO (Arbel & Mairal, 2022), and SOBA/SABA (Dagréou et al., 2022), enabling qNBO to extend effectively from deterministic to stochastic settings.

To elaborate, qNBO consists of three components, with the first two employing quasi-Newton recursion schemes. A straightforward stochastic adaptation involves replacing deterministic quasi-Newton recursion schemes with stochastic variants, such as K-BFGS [1] or Stochastic Block BFGS [2]. Additionally, in Part 3 of qNBO, deterministic gradients are replaced with stochastic gradients throughout the iterative process, consistent with the transition from deterministic to stochastic methods discussed by Dagréou et al. (2022).

However, this straightforward extension may compromise convergence rates and computational efficiency under the same smoothness assumptions as in the deterministic setting. For an intuitive comparison (without using quasi-Newton methods), refer to the analysis of SOBA (Dagréou et al., 2022) and its improvement by MA-SOBA (Chen et al., 2023, Optimal Algorithms for Stochastic Bilevel Optimization under Relaxed Smoothness Conditions).

Addressing your question about maintaining comparable convergence rates and computational efficiency in a stochastic setting, we consider this a promising research direction. The main challenges in achieving this include:

$**Constructing effective estimators:**$ How can we develop unbiased or biased estimators in Part 3 of qNBO by integrating techniques such as variance reduction or momentum (e.g., moving averages)? Potential candidates for these estimators include SGD, SAGA, SARAH, STORM, PAGE, and even Adam. However, analyzing this extension, particularly its convergence properties, requires significant effort. For recent progress (without using quasi-Newton methods), see SRBA (Dagréou et al., 2024) and SPABA (Chu et al., 2024).

$**Analyzing convergence rates and complexity:**$ How can we evaluate the proposed stochastic algorithms in a bilevel setting while incorporating noisy second-order (curvature) information?Addressing these difficulties may require theoretical breakthroughs that go beyond existing first-order techniques, leaving this an unresolved challenge and an interesting future extension.

[1] Goldfarb, D., Ren, Y., & Bahamou, A. (2020). Practical quasi-newton methods for training deep neural networks.

[2] Gower, R., Goldfarb, D., & Richtárik, P. (2016). Stochastic block BFGS: Squeezing more curvature out of data.

Weaknesses: The paper is challenging to follow, with a vague explanation of the quasi-Newton recursion scheme. The design of the proposed method appears overly complex compared to existing methods, and there is limited explanation regarding the validity of $u_{k+1}$ as an accurate approximation for the Hessian-vector product.

$**Response:**$ Thank you for your comment. Appendix B provides a comprehensive explanation of the quasi-Newton recursion scheme, including the algorithmic framework and key details omitted from the main text for brevity. Below, we elaborate on the necessity of subroutine $\mathcal{B}$ and and present explicit lemmas to justify why $u_{k+1}$ serves as an accurate approximation of the Hessian-vector product.

In quasi-Newton methods, the convergence of the quasi-Newton matrix $H^{-1}$ to the true Hessian matrix $\nabla_{yy} f$ depends on strong assumptions [1], which are often not satisfied in practice [2]. Nevertheless, even without these assumptions, the quasi-Newton matrix $H^{-1}$ can still converge to the true Hessian $\nabla_{yy} f$ along specific directions [1,3]. Accordingly, in Step 2 of Algorithm 1, we apply the quasi-Newton recursion along the direction $\nabla_y F$ to compute $u_{k+1}$ , which serves as an approximation of the Hessian-vector product $[\nabla_{yy} f]^{-1} \nabla_y F$ . Recent studies [4,5,6] have explicitly characterized the accuracy of the quasi-Newton matrix in approximating the Hessian along specific directions. Building on these findings, Lemmas D.16 (Lemma D.18 in the revised paper) and D.21 (Lemma D.23 in the revised paper）in Appendix D demonstrate how $u = H^{-1} \nabla_y F$ effectively approximates $[\nabla_{yy} f]^{-1} \nabla_y F$ .

[1] Nocedal, Jorge, and Stephen J. Wright, Numerical optimization[M]. New York, NY: Springer New York, 1999.

[2] Mannel, F. On the convergence of the Broyden-like matrices[J]. 2020.

[3] Dennis J E, Moré J J. A characterization of superlinear convergence and its application to quasi-Newton methods[J]. Mathematics of computation, 1974, 28(126): 549-560.

[4] Rodomanov A, Nesterov Y. Rates of superlinear convergence for classical quasi-Newton methods[J]. Mathematical Programming, 2022: 1-32.

[5] Rodomanov A, Nesterov Y. New results on superlinear convergence of classical quasi-Newton methods[J]. Journal of optimization theory and applications, 2021, 188: 744-769.

[6] Jin Q, Mokhtari A. Non-asymptotic superlinear convergence of standard quasi-Newton methods[J]. Mathematical Programming, 2023, 200(1): 425-473.