A Lazy Hessian Evaluation Framework for Accelerating Stochastic Bilevel Optimization

Comment 3: Normally, strong convexity will have a better convergence rate or better sample complexity. In the current version, both non-convex and strongly convex settings have the same complexity order in the partial gradients in your corollaries. Why is that?

Our Response: Thank you for your question. After addressing the typo mentioned earlier in the results for the strongly convex setting, we find that the strongly convex case requires $\mathcal{O}(\epsilon^{-1}\log\epsilon^{-1})$ partial gradient evaluations, which is fewer than the $\mathcal{O}(\epsilon^{-2})$ partial gradient evaluations required in the non-convex setting. In addition, the batch sizes $\left| \mathcal{D}^{g_{y}} \right|$, $\left| \mathcal{D}^{f_{y}} \right|$, $\left| \mathcal{D}^{g_{yy}} \right|$ required in the strongly convex setting are smaller than those required in the non-convex setting.

Comment 4: The explanation in the experimental part is not clear. For example, table 2 shows the time to "converge" but there is no explanation of what the standard is (such as loss change is smaller than a threshold). Similarly, table 3 presents the number of hessian computing when reaching the same "training loss". Then how much? I suggest announcing those parts more clearly.

Our Response: Thank you for your suggestion. In Table 2, convergence is defined as the point where the training loss reaches its minimum and does not decrease further for 5 consecutive iterations. For Table 3, the training loss criteria for Task 1 is 2.35 and for Task 2 is 0.71. We will clarify these in the revised version of our paper.

Comment 5: In Figure 3, can the authors provide a plot with #hessian computing as the x-axis?

Our Response: Thank you for the suggestion. As per your suggestion we have added a figure using the number of Hessian computations as the x-axis in the revised version of our paper (see Figure 3(b)).

Weaknesses:

Comment 1: The paper does not seem very novel to me since the "lazy" strategies have already been tried though they may not be the lazy hessian in Bilevel optimization. I understood that no former works provided theoretical guarantees of Lazy hessian methods in BO. After I went through the proof details, it looked straightforward. However, I am open to changing my mind if I missed some important parts.

Our Response: Thank you for your comment. We would like to emphasize that the design of our LazyBLO algorithm is far from a straightforward adaptation of single-level optimization to the bilevel counterpart. Compared to single-level optimization, bilevel optimization is fundamentally different and has a far more complex problem structure. Specifically, in bilevel optimization, Hessian information is essential for hypergradient computation, rather than serving as an optional enhancement as in second-order methods for single-level optimization.

Moreover, due to the far more complex problem structure, analyzing the use of lazy Hessian in bilevel optimization is far more challenging than that in single-level optimization settings. There are two major challenges in our theoretical analysis:

1. Increased Stochastic Hypergradient Error: The use of lazy Hessian evaluation introduces additional errors in the stochastic hypergradient estimator for the upper-level function.
2. Error Propagation in the Coupled Structure: Due to the hierarchical and coupled structure of bilevel optimization problems, the error resulting from the stochastic hypergradient estimator with stale Hessian information further propagates to and increases the approximation error of $\mathrm{y}^*(\mathrm{x})$, which is the optimal solution of the lower-level problem, and the approximation error of $\mathrm{z}^*(\mathrm{x},\mathrm{y})$, which is the optimal Hessian-inverse estimation.

Note that $\mathrm{z}^*(\mathrm{x},\mathrm{y})$ is inherently linked to $\mathrm{y}^*(\mathrm{x})$, creating a tight coupling between these errors that significantly complicates the convergence analysis of LazyBLO. These complications in error propagation and approximation are unseen in single-level optimization convergence analysis and even in conventional bilevel optimization algorithm analysis, highlighting the novelty of our theoretical contributions.

Comment 2: The order selections in the corollaries are a little confusing to me. In Corollary 5.4, I assume that N is chosen as a constant, T is at the order of $\epsilon^{-1}$, and batch sizes are all $\epsilon^{-1}$. Then LazyBLO requires $\epsilon^{-2}$ partial gradients but why $\epsilon^{-2}$ for hessian evaluations? Is it $\mathcal{O}(T)=\mathcal{O}(\epsilon^{-1})$. The selection in Corollary 5.7 is more unclear to me. With the choice of a constant level of batch sizes of $D^{f_{y}}$ and $D^{g_{yy}}$, N cannot be constant and I did not get your results. Correct me if I misunderstand it.

Our Response: Thanks for your comment. In Corollary 5.4, the number of hessian evaluations is calculated as the product of the number of iterations $T$ and the batch size $D^{g_{yy}}$. Specifically:

• $\mathcal{O}(T)=\mathcal{O}(\epsilon^{-1})$, and
• $\left| \mathcal{D}^{g_{yy}} \right|=\Theta(\epsilon^{-1})$.

This results in $\mathcal{O}(\epsilon^{-2})$ Hessian evaluations.

In Corollary 5.7, we acknowledge a typo in the original paper. Specifically, we missed a $1/N$ term in the theorem's expression. Based on the corrected inequality: $\frac{1}{N}\sum_{n=0}^{N-1}\mathbb{E}\bigg[\ell\left(\mathrm{x}_t^n\right)-\ell^*\bigg] \leq \left(1-\mu_f \alpha\right)^t\hat{\Delta}_0$

$\qquad \qquad \qquad \qquad \qquad \quad +\mathcal{O}\bigg(\sigma_{g_{xy}}^2 + \sigma_{f_{x}}^2 + \frac{1}{N^{4}}\sigma_{g_{yy}}^2 +\frac{1}{N^{4}} \sigma_{f_{y}}^2 + \frac{1}{N}\sigma_{g_{y}}^2 \bigg),$ where $\hat{\Delta}_0=\frac{1}{N} \sum _{n=0}^{N-1}\left(\ell\left(\mathrm{x}_0^n\right)-\ell^*\right)+\frac{1}{N} \sum _{n=0}^{N-1}\left\|\mathrm{y}_0^n-\mathrm{y}^*\left(\mathrm{x}_0^n\right)\right\|^2+\frac{1}{N}\left\|\mathrm{z}_0-\mathrm{z}_0^*\right\|^2,$

we have $\mathcal{O}(T)=\mathcal{O}(\log\epsilon^{-1})$.

By choosing

• $\left| \mathcal{D}^{f_{x}} \right|,\left| \mathcal{D}^{g_{xy}} \right|=\Theta\left( \epsilon^{-1} \right)$,
• $\left| \mathcal{D}^{f_{y}} \right|,\left| \mathcal{D}^{g_{yy}} \right|=\Theta\left( N^{-4}\epsilon^{-1} \right)$, and
• $\left| \mathcal{D}^{g_{y}} \right|=\Theta\left( N^{-1}\epsilon^{-1} \right)$,

LazyBLO requires $\mathcal{O}(\epsilon^{-1}\log\epsilon^{-1})$ partial gradient evaluations and $\mathcal{O}(N^{-4}\epsilon^{-1}\log\epsilon^{-1})$ HVP evaluations to reach an $\epsilon$-optimal point which significantly improves the performance compared to the non-convex setting. We have corrected the typo in the revised version of the paper.

Questions:

Comment 1: While the problem of reducing Hessian computations is valid, the main technical novelty of LazyBLO is unclear. The approach reduces Hessian computation, which may seem like an incremental adaptation of existing methods rather than a novel algorithm. It would be helpful if you can specify them.

Our Response: Thank you for your comment. While reducing Hessian computations has been explored in single-level optimization to accelerate second-order methods, this is the first time such an approach has been applied to bilevel optimization. Before our work, it was unclear whether reducing Hessian computations was feasible or not in bilevel optimization due to the complexity of the bilevel problems.

Importantly, our lazy Hessian algorithm design is far from being a simple adaptation from single-level optimization to the bilevel counterpart. Compared to single-level optimization, bilevel optimization is fundamentally different and has a far more complex problem structure. Specifically, in bilevel optimization, Hessian information is essential for hypergradient computation, rather than being an optional enhancement as in second-order methods for single-level optimization.

Furthermore, due to the far more complex problem structure, analyzing the use of lazy Hessian in bilevel optimization introduces substantial challenges that are not present in single-level optimization settings. These challenges include:

1. Increased Stochastic Hypergradient Error: The use of lazy Hessian evaluation introduces additional errors into the stochastic hypergradient estimator for the upper-level function.
2. Error Propagation in the Coupled Structure: Due to the hierarchical and coupled structure of bilevel optimization problems, the error resulting from the stochastic hypergradient estimator with stale Hessian information further propagates to and increases the approximation error of $\mathrm{y}^*(\mathrm{x})$, which is the optimal solution of the lower-level problem, and the approximation error of $\mathrm{z}^*(\mathrm{x},\mathrm{y})$, which is the optimal Hessian-inverse estimation.

Note that $\mathrm{z}^*(\mathrm{x},\mathrm{y})$ is inherently linked to $\mathrm{y}^*(\mathrm{x})$, creating a tight coupling between these errors that significantly complicates the convergence analysis of LazyBLO. These complexities in error propagation and approximation are unseen in single-level optimization convergence analysis, even in the conventional bilevel optimization algorithm analysis.

Comment 2: Could you provide more insights into how LazyBLO performs on problems where Hessian is unstable?

Our Response: Thanks for your question. We assume your comment that "Hessian is unstable" refers to the situations where the Hessian matrix varies rapidly over time. If this is the case (e.g., in highly non-linear models), LazyBLO can be adapted by increasing the frequency of Hessian-vector product evaluations. This flexibility allows the algorithm to remain robust and perform effectively even in these challenging scenarios.

Questions:

Comment 1: Since the author assumes the strong convexity of the inner problem, why the challenge 2 is still challenging?

Our Response: Thank you for your question. While it is true that strongly convex inner problems are generally easier to solve compared to non-convex ones, the assumption of strong convexity does not imply that these problems can always be solved in closed form, except in a few very specific cases. In most scenarios, well-designed numerical procedures are still required to solve the inner problems approximately. In addition, the approximation to which the lower level problem is solved determines the bias in the implicit gradient computation of the objective function. As a consequence, this bias needs to be controlled carefully (both theoretically and numerically) to ensure convergence of the algorithm.

Comment 2: The method only considers the computation of HVP. Can the author further explain how to reduce the computation of JVPs as mentioned in line 308?

Our Response: Thanks for your question. Rewrite Equation (5) in our paper as $h_{t,n}^{f} = \nabla_{x} f \left( x_{t}^{n}, y_{t}^{n}; D_{t,n}^{f_{x}} \right) + H_{t}^{j}z_{t}$, where $H_{t}^{j}:=\nabla_{x y}^{2} g \left(x_{t}^{N}, y_{t}^{N}; D_{t}^{g_{xy}}\right)$. Then, it can be seen that JVP (i.e., $H_{t}^{j}z_{t}$) would be updated less frequently compared to HVP, since it is in the outer loop of the algorithm. If one wants to reduce the computation JVP, a similar lazy approach can also be applied on the JVP update in LazyBLO.

Comment 3: Theorem 5.3 indicates that N=1 results in the optimal convergence rate, which is SOBA. Then how can the proposed method improve the convergent in theory? Combining the weakness mentioned above, the proposed method has marginally improved.

Our Response: Thanks for your comment. We would like to clarify that when $N=1$, our proposed method recovers the result of AmIGO, not SOBA. This is because SOBA relies on an additional assumption—third-order Lipschitz continuity—which our method does not require. Furthermore, since we use a lazy approach where infrequent Hessian evaluations incurs additional errors, our goal is not to outperform the non-lazy approaches. Rather, our goal is to be as close to the non-lazy counterpart approaches as possible. Due to the additional errors accumulated because of the use of these stale Hessians, approximate HVP evaluations, and the coupling hierarchical structure of the bilevel problems, it is unclear whether LazyBLO will converge or not. In fact, establishing convergence for the lazy approach proposed in this work is itself surprising and a major contribution of our work.

Comment 4: Can we author explain what is the main reason for such a significant improvement in Table 2 since the theory can not show such an advantage?

Our Response: Thank you for your comment. The convergence rate alone can sometimes be misleading, as it is based on the number of iterations, without accounting for the computational cost of each iteration. To provide a more objective comparison, we measure performance using wall-clock time. LazyBLO achieves significant improvements in wall-clock time because it performs infrequent Hessian evaluations, which greatly reduce computational overhead. This efficiency explains the notable improvement observed in Table 2, even if the theoretical convergence rate does not directly reflect this advantage.

Weaknesses:

Comment 1: Theorem 5.3 looks like the algorithm has some trouble in convergence when the gradient estimation variance is not arbitrarily small. The only way to solve this problem is to take a sample size that is large enough. Previous methods allow such error terms to diminish with iteration. This is a big disadvantage of this algorithm. As a result, Corollary 5.4 shows that the proposed method has bare improvement in partial gradient and HVP evaluations.

Our Response: Thanks for your comment. We would like to address your concerns and provide additional context as follows:

1. The large batch size is required only for the theoretical convergence guarantees. In our experiments, we use a small batch size and still show fast-convergence of LazyBLO compared to the state-of-the-art baselines. This practical evidence indicates that the theoretical requirement of a large batch size does not negatively impact the real-world applicability of the proposed method.
2. The need for large batch sizes is not unique to LazyBLO. Many existing bilevel optimization algorithms, such as stocBiO, AmIGO, and iNEON [R1], also adopt large batch sizes in their theoretical analyses to control variance. Additionally, our approach employs SGD to approximate Hessian-vector product by solving a quadratic problem. Previous bilevel optimization algorithms using the same approximation are AmIGO and SOBA. Our approach shares a similarity with AmIGO, which avoids strong third-order Lipschitz continuity assumptions at the cost of requiring large batch sizes. In contrast, SOBA relies on these stronger assumptions to manage variance. In addition, previous methods, such as BSA, TTSA, SUSTAIN, MSTSA, VRBO, MRBO and ALSET require all the stochastic functions to satisfy Assumptions 3.1 and 3.2 stated in the manuscript to avoid large batch gradient computations. Note that this is certainly a very strong assumption and is not required in our analysis of LazyBLO. This indicates that large batch sizes are a common trade-off in bilevel optimization methods, rather than a disadvantage specific to LazyBLO.
3. Since we use a lazy approach where we perform infrequent Hessian evaluations that incur additional errors, our aim is not to outperform the non-lazy approaches. Rather, our goal is to be as close to the non-lazy counterpart approaches as possible. Due to the additional errors accumulated because of the use of these stale Hessians, approximate HVP evaluations, and the coupling hierarchical structure of the bilevel problems, it is unclear whether LazyBLO will converge or not. In fact, establishing convergence for the lazy approach proposed in this work is itself surprising and a major contribution of our work.

[R1] Huang, M., Ji, K., Ma, S., & Lai, L. (2022). Efficiently escaping saddle points in bilevel optimization. arXiv preprint arXiv:2202.03684.

Comment 2: It may be necessary to assume $l^*$ to be lower bounded for lines 372 and 428.

Our Response: Thank you for your comment. We will include the lower-bounded assumption of $l^*$ in the revised version of our paper.

Comment 3: The motivation for reducing only HVP computations, rather than both HVP and JVP computations, is unclear. In fully single-loop bilevel methods, both HVP and JVP are computed once per iteration, unlike the Neumann series approach discussed in this paper. There is no apparent reason to prioritize reducing only HVP computations. HVP and JVP should be considered equally important, particularly in single-loop methods, which form the basis of the proposed approach. Besides, it is better to provide the computation complexity comparison of the proposed method in terms of HVP and JVP.

Our Response: Thank you for your comment. We provide the following reasons for prioritizing the reduction of HVP computations:

1. Most bilevel optimization algorithms follow a double-loop structure and require multiple HVP computations per iteration.
2. While single-loop algorithms represent a subset of bilevel optimization methods and often compute HVPs only once per iteration, there are notable exceptions. Some single-loop algorithms, such as SUSTAIN, TTSA, BSA, and ALSET, require multiple HVP computations per iteration. In contrast, JVPs are computed only once in all these algorithms, regardless of the number of HVP computations per iteration.

Given these points, HVP computations are at least as frequent as JVP computations and incur a more significant computational bottleneck, especially in double-loop and certain single-loop bilevel algorithms. As such, reducing the number of HVP computations is more crucial than reducing JVP computations.

Questions:

Comment 1: It would be better to add the complexity of the proposed method in Table 1 for comparison as well.

Our Response: Thank you for your suggestion. We have included the complexity of our proposed method in Table 1 in the revised version of the paper.

Weaknesses:

Comment 1: The convergence rate of the proposed method depends on increasing batch sizes. However, the current state-of-the-art theoretical guarantees for stochastic bilevel optimization do not require this, as demonstrated by methods such as ALSET, SOBA, FSLA, F$^2$SA, and etc.

Our Response: The need for increasing batch sizes is not unique to LazyBLO. Many existing bilevel optimization algorithms, such as stocBiO, AmIGO, and iNEON [R1], also adopt increasing batch sizes in their theoretical analyses to control variance. Additionally, our approach employs SGD to approximate Hessian-vector product by solving a quadratic problem. Previous bilevel optimization algorithms using the same approximation are AmIGO and SOBA. Our approach shares a similarity with AmIGO, which avoids strong third-order Lipschitz continuity assumptions at the cost of requiring increasing batch sizes. In contrast, SOBA relies on these stronger assumptions to manage variance. In addition, previous methods, such as BSA, TTSA, SUSTAIN, MSTSA, VRBO, MRBO and ALSET require all the stochastic functions to satisfy Assumptions 3.1 and 3.2 stated in the manuscript to avoid increasing batch gradient computations. Note that this is certainly a very strong assumption and is not required in our analysis of LazyBLO. This indicates that increasing batch sizes are a common trade-off in bilevel optimization methods, rather than a disadvantage specific to LazyBLO.

It is also important to note that ALSET adopts a different strategy, using Neumann series to approximate the Hessian inverse, which distinguishes it from LazyBLO. While FSLA and F$^2$SA do not require increasing batch sizes, they introduce additional assumptions. Specifically, Assumptions A and B in FSLA are also hold for their stochastic version, which are not required in our analysis; and Assumptions 4 and 5 in F$^2$SA are not required for our LazyBLO.

[R1] Huang, M., Ji, K., Ma, S., & Lai, L. (2022). Efficiently escaping saddle points in bilevel optimization. arXiv preprint arXiv:2202.03684.

Comment 2: Although the proposed method achieves the same convergence rate to fully single-loop bilevel algorithms with respect to $\epsilon$ for general nonconvex–strongly convex bilevel problems, the HVP computational complexity appears to be $T=\mathcal{O}(N\epsilon^{-1})$, which is even higher than that of fully single-loop bilevel algorithms. This implies that the lazy strategy not only offers no theoretical improvement, but also is even worse than its non-lazy counterpart. Thus, the results presented here are not unexpected.

Our Response: Thank you for your comment. We would like to clarify that the order of the computational complexity of our LazyBLO is $-2$, which matches the state-of-the-art SGD-based bilevel optimization algorithms. In addition, since we use a lazy approach where infrequent Hessian evaluations incurs additional errors, our goal is not to outperform the non-lazy approaches. Rather, our goal is to be as close to the non-lazy counterpart approaches as possible. Due to the additional errors accumulated because of the use of these stale Hessians, approximate HVP evaluations, and the coupling hierarchical structure of the bilevel problems, it is unclear whether LazyBLO will converge or not. In fact, establishing convergence for the lazy approach proposed in this work is itself surprising and a major contribution of our work.

Thank you for your response. For the first question, I have reviewed the references you provided, and it appears that methods like TTSA, SUSTAIN, and ALSET do not impose any Lipschitz or strong convexity assumptions on the stochastic function, which aligns with your Assumptions 3.1–3.2. Furthermore, these methods also do not require an increasing batch size. For the third-order Lipschitz continuity assumption in SOBA, it is because they use fully single loop so that this assumption is needed to avoid the increasing batch size. However, the method you used is double loop, removing increasing batch size does not need the third order Lipschitz continuity; See TTSA, SUSTAIN, and ALSET. Therefore, I believe that eliminating the increasing batch size assumption would be necessary to achieve state-of-the-art theoretical performance.

On the second point, I still have some reservations. Fully single-loop algorithms such as SOBA and FLSA represent a recent trend in bilevel optimization as they avoid frequent HVP computations per iteration, aligning them with the direction of your work. Comparing against double-loop algorithms, which are not considered the current baselines, seems less convincing as a justification for focusing on the reduction of HVP rather than both HVP and JVP. Are there specific challenges or limitations that prevent skipping JVP entirely in your framework?

Based on the current response, I would like to keep my score.