Memory-Efficient Gradient Unrolling for Large-Scale Bi-level Optimization

Authors: qianli  shen
yezhen  wang
zhouhao  yang
xiang  li
haonan  wang
yang  zhang
jonathan  scarlett
zhanxing  zhu
kenji  kawaguchi
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 08:38:43

Abstract: Bi-level optimizaiton (BO) has become a fundamental mathematical framework for addressing hierarchical machine learning problems.
As deep learning models continue to grow in size, the demand for scalable bi-level optimization has become increasingly critical.
Traditional gradient-based bi-level optimizaiton algorithms, due to their inherent characteristics, are ill-suited to meet the demands of large-scale applications.
In this paper, we introduce **F**orward **G**radient **U**nrolling with **F**orward **G**radient, abbreviated as **$($FG$)^2$U**, which achieves an unbiased stochastic approximation of the meta gradient for bi-level optimizaiton.
$($FG$)^2$U circumvents the memory and approximation issues associated with classical bi-level optimizaiton approaches, and delivers significantly more accurate gradient estimates than existing large-scale bi-level optimizaiton approaches.
Additionally, $($FG$)^2$U is inherently designed to support parallel computing, enabling it to effectively leverage large-scale distributed computing systems to achieve significant computational efficiency.
In practice, $($FG$)^2$U and other methods can be strategically placed at different stages of the training process to achieve a more cost-effective two-phase paradigm.
Further, $($FG$)^2$U is easy to implement within popular deep learning frameworks, and can be conveniently adapted to address more challenging zeroth-order bi-level optimizaiton scenarios.
We provide a thorough convergence analysis and a comprehensive practical discussion for $($FG$)^2$U, complemented by extensive empirical evaluations, showcasing its superior performance in diverse large-scale bi-level optimizaiton tasks.

Main paper: 
 introduction: 
 Bi-level optimization is a mathematical framework with a long history of research [10, 65, 73] , dealing with hierarchical optimization problems where one problem is nested within the other. A bi-level optimization problem can be formulated as:EQUATIONwhere θ ∈ Θ ⊆ R M denotes the inner parameter, ϕ ∈ Φ ⊆ R N denotes the meta parameter, and f , g are called the meta objective function and inner objective function, respectively. Method/Unrolled Depth The performance of (FG) 2 U surpasses that of other large-scale bi-level optimization methods, owing to its accurate gradient approximation, while demonstrating better memory efficiency. Bottom Right: Efficiency tradeoff of (FG) 2 U on Data Condensation experiments. The efficiency of (FG) 2 U can be well enhanced via intra/inter-GPU parallelism.precipitating challenges associated with gradient-based bi-level optimization. Consequently, various gradient-based bi-level optimization algorithms have been developed [73] . These algorithms typically employ an iterative solution θ T obtained by executing multiple inner optimization steps to approximate the meta gradient, and provide different tradeoffs between computational costs and performance for meta gradient approximation.However, as the scale of deep learning models continues to expand, the requirements for scalability in bi-level optimization correspondingly increase. Existing gradient-based bi-level optimization algorithms, due to their inherent characteristics, are ill-suited to meet the demands of large-scale applications. Concretely, gradient unrolling (GU) methods [17, 16, 40, 60] are bottlenecked by the memory overhead associated with either the dimension of the inner parameter or the number of iterative steps for the inner problem. Implicit Function (IF) approaches [48, 19, 64, 76] are compromised by approximation errors, which stem from the iterative estimation of inner solutions and computations that involve the Hessian matrix. Value Function (VF) based strategies [39, 37, 61, 33] , although exhibit commendable theoretical properties [8] for deterministic bi-level optimization, have yet to gain traction in practical applications, predominantly due to their limitations in addressing largescale stochastic challenges [73] . Recent advancements in algorithms [60, 9] have been specifically tailored for large-scale bi-level optimization. Although these methodologies facilitate efficient gradient approximation by compromising accuracy, they may result in significantly suboptimal performance due to biased gradient approximations. Additionally, these methods struggle in more complex scenarios, such as when inner problems are addressed through black-box optimization.In this paper, we propose a novel method called Forward Gradient Unrolling with Forward Gradient, abbreviated as (FG) 2 U, which achieves an unbiased stochastic approximation of the meta gradient for bi-level optimization. (FG) 2 U circumvents the memory issues associated with GU-based approaches and approximation issues associated with IF-based approaches. Compared to recently developed large-scale bi-level optimization approaches, (FG) 2 U delivers significantly more accurate gradient estimates. Additionally, (FG) 2 U is inherently designed to support parallel computing, enabling it to effectively leverage large-scale distributed computing systems to achieve significant computational efficiency. In practice, a cost-effective two-phase paradigm can be achieved by strategically placing (FG) 2 U and other methods at different stages of the training process to balance efficiency and performance. Further, (FG) 2 U is easy to implement within popular deep learning frameworks, and can be conveniently adapted to address more challenging zeroth-order bi-level optimization scenarios.We provide an overview of (FG) 2 U in Figure 1 to illustrate its strengths and role in large-scale bi-level optimization. The rest of the paper is organized as follows. Firstly, in Section 2, we provide summaries of existing bi-level optimization algorithms and discuss their limitations in large-scale contexts. Next, in Section 3, we introduce the proposed method, (FG) 2 U, followed by a convergence analysis in Section 3.1 and a detailed discussion of the practical considerations in Section 3.2. Further, in Section 4, we conduct extensive empirical studies covering large-scale bi-level optimization in computer vision, natural language processing, and physics-informed machine learning to demonstrate the efficacy of (FG) 2 U in large-scale bi-level optimization scenarios.To address the memory issue of GU methods, truncated Reverse Gradient Unrolling (TRGU) [60] is proposed to reduce the memory usage by preserving only the last K steps of the inner optimization trajectory. However, this introduces a significant bias in large-scale scenarios, particularly when the permissible K is small. Recall ( 5) and (6) , where the conventional RGU method computes the hypergradient by fully unrolling the T -step inner optimization into a computational graph. Instead, TRGU performs s-step truncated back-propagation and approximates the gradient with the intermediate term cT -s:EQUATIONAccording to Proposition 3.1 in [60] , if the inner-level objective function g is L-smooth, twice-differentiable and globally α-strongly convex, and the gradient update rule writes θt = θt-1 -η∇ θ g(θt-1, ϕ), then the bias of s-step TRGU would be bounded byEQUATIONThe bound (22) demonstrates an exponentially decaying rate in s over the bias of s-step TRGU. However, when s gets smaller, which means that we truncate the computational graph heavier in pursuit of lower memory cost, the bias would grow exponentially. This would result in an inaccurate calculation of the hypergradient. Contrastively, our (FG) 2 U is an unbiased estimator of the hypergradient, while still keeping high memory efficiency with a small sample size of forward gradient as in (10) . 
background: 
 Gradient-based Bi-level Optimization. Within deep learning applications, the model concerned with optimizing over θ as presented in (1) typically constitutes deep neural networks. The optimal parameters of such networks are not explicitly accessible and are estimated through iterative procedures.Consequently, the primal problem of bi-level optimization in ( 1) is approximately reformulated as follows:EQUATIONwhere θ 0 (ϕ) = Ω 0 (ϕ), θ t (ϕ) = Ω t (θ t-1 (ϕ), ϕ) ∈ Θ, t = 1, . . . , T,where Φ ⊆ R N , Θ ⊆ R M are the parameter spaces; T , commonly called the unrolled depth, denotes the number of inner optimization steps for approximating θ * (ϕ); Ω 0 : R N → R M specifies the initialization of the inner optimization, and Ω t : Θ × Φ → Φ delineates the transition dynamics of the inner optimization at timestep t. In particular, for gradient descent,Ω t (θ t-1 (ϕ), ϕ) = θ t-1 -η t ∇ θ g(θ t-1 , ϕ), where η t denotes the step size at timestep t.To optimize ϕ using a first-order method, it is necessary to estimate the meta gradient ∇ ϕ h, which can be further decomposed according to the chain rule:EQUATIONThe computation of meta-gradient poses a significant challenge, primarily due to the need for efficient approximation of the implicit gradient. This task is complicated by the recursive dependency of θ T on ϕ. To surmount this challenge, a variety of gradient-based bi-level optimization algorithms have been developed, as extensively reviewed recently in [73] . These algorithms can be fundamentally categorized into three types based on their approach to meta-gradient approximation: Gradient Unrolling (GU), Implicit Function (IF), and Value Function (VF). Recent innovations such as truncated RGU (TRGU) [60] and Hessian-Free approaches [76, 75, 9] , which are predicated on GU and IF methodologies respectively, have introduced significant biases in their approximations to accommodate the computational constraints of large-scale scenarios. In the subsequent paragraph, we furnish a concise overview of GU-based approaches, addressing their non-constant memory issues in large-scale applications. Extended discussions on the remaining methods are reserved for Appendix B.Gradient Unrolling. The core idea behind GU [17, 16, 40, 60] entails unrolling the inner optimization into an expansive computational graph, followed by the employment of automatic differentiation (AD) techniques for the iterative computation of gradients.Forward Gradient Unrolling (FGU) [17, 16] computes the meta gradient using the following forward recursive formula, starting from Z 0 = dΩ0(ϕ) dϕ :EQUATIONReverse Gradient Unrolling (RGU) [46, 16] , instead of the employment of explict reccursive formulas of Z T , focuses on the implicit reccursive formulas of ∇ ϕ h:EQUATIONThe corresponding reverse recursive formulas can thus be summarized asEQUATIONWeakness (GU): Non-Constant Memory. Both GU approaches exhibit a non-constant memory overhead, which constrains their utility in large-scale scenarios. The forward reccursive formulas in (4) revolve around the Jacobian matrix product, demanding O(M N ) space consumption. The reverse recursive formulas in (6) necessitate the storage of the entire trajectory of the inner optimization θ 0:T for backward computation, thereby imposing a memory requirement of O(T M ). These requirements are often impractical for large-scale bi-level optimization, when ϕ and θ are of high dimension and a significant unrolled depth is required.Forward Gradient. Forward-mode automatic differentiation (forward-mode AD) has been applied to a variety of research fields, including the training of recurrent neural networks [70] , the computation of Hessian vector products [50] , etc. However, the computation of the true gradient via forward-mode AD requires the full Jacobian, which is typically too costly to compute.To solve this, forward gradient learning [69, 4, 63, 4, 56] , built upon forward-mode AD, was proposed. Forward gradient methods update parameters based on the directional gradient along a random perturbation direction for backpropagation-free training. More formally, given a differentiable function h : R N → R, the gradient for a given input ϕ ∈ R N can be approximated as ] = ∇h(ϕ)I = ∇h(ϕ). Despite the unbiasedness of ∇h, the dimension-dependent variance of the estimated gradient with a single direction impedes the scaling-up to high-dimensional problems. In practice, Monte Carlo gradient estimation can be used via averaged forward gradients over multiple random directions to reduce the variance.EQUATION 
(fg) 2 u: forward gradient unrolling with forward gradient: 
 We aim to circumvent the memory overhead issues associated with forward gradient unrolling (FGU) as discussed in Section 2. We begin by examining the forward gradient of h at ϕ,EQUATIONwhere v ∼ p(v) is a N -dimensional multivariate random variable, satisfying E[vv T ] = I. We follow the idea of FGU introduced in Section 2 to compute Z T v. By multiplying both sides of (4) by v on the right, we can obtain the recursive formulas for Z t v asEQUATIONThe revised recursive formulas in (9) facilitate the tracking of a M -dimensional vector Z t v, rather than full Jacobian Z t of size M × N , throughout the forward pass. The stochastic estimation in (8) is unbiased, adhering to the properties of forward gradient methods. To reduce the variance, we use Monte Carlo estimate via averaged forward gradients over b i.i.d. random directions:EQUATIONWe call this algorithm (FG) 2 U, as an abbreviation of Forward Gradient Unrolling with Forward Gradient. The algorithm is summarized in Appendix A as Algorithm 1.Compared to GU-based methods, as discussed in Section 2, (FG) 2 U eliminates the dependency on the meta parameter dimension N and the depth of unrolling T without introducing bias, significantly enhancing memory efficiency. Unlike IF-based methods, as discussed in Appendix B.2, (FG) 2 U overcomes the approximation issues associated with them while maintaining a constant memory overhead, thus providing superior gradient approximation. Compared to TRGU and Hessian-Free methods, which compromise approximation accuracy for efficiency, (FG) 2 U consistently delivers accurate gradient approximations. The computational efficiency of (FG) 2 U can be further enhanced by leveraging large-scale distributed computing resources, capitalizing on its inherently parallelizable formulation as presented in (10) . In practice, a more cost-effective two-phase paradigm can be achieved by strategically placing (FG) 2 U and other methods at different stages of the training process, as we will discuss in Section 3.2. For an illustration of the role of (FG) 2 U in large-scale bi-level optimization, please refer to Figure 1 .In this section, we provide a convergence analysis for (FG) 2 U. The proofs can be found in Appendix C.First, we establish a bound on the variance of the estimated gradient, when employing random vectors whose entries follow the Rademacher distribution. Lemma 3.1. For any ϕ ∈ Φ, if v i ∼ Unif({-1, 1} N ), the gradient estimation in (10) , satisfiesE∥ ∇h(ϕ) -∇h(ϕ)∥ 2 = 1 ρ ∥∇h(ϕ)∥ 2 ,where ρ := b N -1 ∈ (0, 1] as the sample size b is selected from 1,• • • , N -1.The resultant error is bounded byO N -1 b, where b represents the sample size used for computing the forward gradient, and N is the dimensionality of the gradient itself. This bound demonstrates how the error scales inversely with the sample size while also being influenced by the gradient's dimensionality.Next, we lay down the following assumptions, on which our main theorems are based. Let ψ = (θ, ϕ) ∈ Θ × Φ denote the combination of the lower-level parameter θ and the meta parameter ϕ. Following existing papers on the theory of bilevel optimization [45, 60, 28] , in Assumption 3.2, we adopt some standard assumptions over the smoothness of the objective functions f and g. Assumption 3.2. The meta objective function f (ψ) and the lower-level objective function g(ψ) are both C-Lipschitz and L-smooth, i.e., for any ψ, ψEQUATIONEQUATIONThe next assumption regulates that the transition functions Ω satisfy similar smoothness conditions. Assumption 3.3. The transition functions Ω 0:T are C Ω -Lipschitz and L Ω -smooth, i.e., for anyEQUATIONEQUATIONAssumption 3.3 is made to ensure the generality of our analysis over different optimizers. Note that Ω is scheme-dependent w.r.t. the gradient-based optimizer we adopt for lower-level problems. In many cases, such as gradient descent where Ω t (ψ t-1 ) = θ t-1 -η t ∇ θ g(ψ t-1 ), Assumption 3.3 is a direct consequence of Assumption 3.2.We propose the following theorem and remark for convergence analysis of (FG) 2 U on problem (2) . Notice the convergence result can be extended to the primal BO problem (1) with some further assumptions. We place a proof scratch and some discussions in Appendix C.3. Theorem 3.4 (Convergence). Suppose that Asumption 3.2 and Assumption 3.3 hold. Setting the learning rate β = ρ (ρ+1)L h for gradient descent over the hyperparameter ϕ, then there exists a constant L h (depending on C, L, C Ω , L Ω , and T , and defined formally in the proof) such thatEQUATION)Remark 3.5. Theorem 3.4 shows that Algorithm 1 converges to an ϵ-accurate stationary point with a convergence rate of = O(ϵ -1 ρ -1 ).Recall that ρ = b N -1 , which indicates that the convergence rate is linearly dependent on N , which poses a significant challenge when managing high-dimensional meta-parameters ϕ. However, it is important to note that the dimension-dependent convergence rate represents an upper bound, and scalability has been found to be feasible with several practical considerations, as discussed in the following subsection.Choice of b. According to the convergence analysis in Section 3.1, a sample size of b = O(N ) is required to achieve a convergence rate of O(ϵ -1 ). However, it has been widely observed that forward gradient and zeroth-order optimization, despite having dimension-dependent convergence rates, work well empirically with b = O(1) in large-scale scenarios, such as in LLM fine-tuning [47, 74] . In this paper, we select the largest possible b that does not exceed the GPU memory limit for our empirical study. Additionally, gradient accumulation is utilized to further control variance and stabilize the training process.Cost-Effective Two-Phase Paradigm. It is important to note that the upper bound delineated in (15) linearly depends on the performance discrepancy between the initialized meta parameter ϕ 0 and the optimal. This dependence motivates the adoption of a more cost-effective two-phase paradigm for large-scale bi-level optimization. In the initial phase, we utilize efficient yet less accurate gradient approximation methods, such as TRGU [60] and Hessian-Free [9] , to efficiently establish an initial ϕ 0 that surpasses random initialization, while keeping computational overhead manageable. Subsequently, in the second phase, (FG) 2 U is utilized for a more accurate, albeit less efficient, gradient approximation to further elevate the performance, leveraging extensive computational resources.Implementation. The technique employed in computing ∇h(ϕ)v is identified as forward-mode automatic differentiation (forward-mode AD). In advanced automatic differentiation libraries, such as JAX [5] and PyTorch [3] , forward-mode AD is efficiently implemented as Jacobian-vector product (jvp), without the necessity of explicitly computing the Jacobian matrix. The FLOP cost of jvp is approximately three times that of a standard forward pass, while the memory overhead is doubled. In practice, it is only necessary to define the forward computational graph of inner optimization and invoke forward-mode AD, which simplifies the implementation process significantly. Regarding distributed training, JAX offers the vmap interface for efficient intra-GPU parallelism and the pmap interface for effective inter-GPU parallelism.Zeroth-order Bi-level optimization. In certain applications of bi-level optimization, the inner problem is approached as a black box, where the gradient of Ω is inaccessible, rendering the analytical gradient unrolling unfeasible. For example, in PDE-constrained optimization [23, 62] , in which the inner problem entails solving a Partial Differential Equation (PDE) using a non-differentiable solver. In such scenarios, rather than employing forward-mode Automatic Differentiation (AD), one can resort to Finite Difference methods to approximate the directional gradient ∇h(ϕ)v byEQUATIONwith sufficiently small positive μ > 0. We refer to this zeroth-order variant of (FG) 2 U as (FG) 2 U-ZO, noting that the computation solely encompasses two forward passes and does not involve the utilization of any first-order information. The memory complexity is the same as forward-mode AD and the actual computation time will be slightly less than forward-mode AD, at the cost of introducing an approximation bias. We give a more detailed discussion within the context of zeroth-order optimization [41] in Appendix D, and empirically study a corresponding case in Section 4.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Justification: The assumptions are included in Section 3.1, and the complete proof is placed in Appendix C. Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and cross-referenced.• All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
experiments: 
 We conduct experiments across various contexts, as detailed in the respective subsections. Initially, we engage in an image data condensation task, where we focus on a comprehensive performance comparison between (FG) 2 U and both classical and large-scale bi-level optimization algorithms. Subsequently, we investigate meta-learning for the online adaptation of language models, employing a GPT model as the inner model, to illustrate how (FG) 2 U effectively circumvents the non-constant memory issue associated with RGU. Finally, we address a physics-informed bi-level optimization problem, where gradient-based inner solvers are ineffective, to demonstrate the efficacy of combining (FG) 2 U-ZO, the zeroth-order variant of (FG) 2 U discussed in Section 3.2, with non-differentiable numerical solvers.Data Condensation. To overcome the challenges posed by large-scale datasets, a line of works known as data condensation [68, 72] has been proposed. The main idea is to generate a compact, synthesized dataset, designed to elicit similar behaviors in machine learning models as those trained with the original, massive dataset. The objective of the mainstream principles [72] designed for data condensation can be naturally formulated as a bi-level optimization problem. We focus on the best-known principle performance matching [72] on classification tasks, which can be formulated asEQUATIONwhere D o , D c respectively denote the original and condensed dataset, θ denotes the model parameter, L denotes the cross-entropy loss function, and η represents the step-size for inner optimization. We conducted our experiments following the standard data condensation setting established by [68, 77, 67] . A more detailed task description is given in Appendix E.1 and implementation details are given in Appendix F.1.The condensed datasets are evaluated using 3-layer convolutional networks with randomly initialized parameters, and the average accuracies on test datasets are summarized in Table 1 . Compared to large-scale bi-level optimization methods like TRGU and Hessian-Free, which prioritize efficiency at the expense of approximation accuracy, (FG) The results of RGU, which represent the upper performance bound for both TRGU and (FG) 2 U, are provided for reference, along with the results from training on the entire dataset (denoted as WHOLE in Table 1 ), representing the upper performance bound for all approaches. However, it is crucial to acknowledge that RGU is not practical in large-scale bi-level optimization scenarios due to its non-constant memory requirements, as discussed in Section 2. This limitation will be further exemplified in the subsequent, where the inner model is significantly larger. In principle, the performance of (FG) 2 U can be further improved to approach that of RGU by increasing the number of random directions for gradient approximation.The memory and computational efficiencies of TRGU, Hessian-Free, and (FG) 2 U in the most challenging case (CIFAR-100, IPC=50) are reported in Figure 1 (Bottom Right), demonstrating that the efficiency of (FG) 2 U can be significantly enhanced through intra/inter-GPU parallelism.Meta Learning Online Adaptation of Language Models. The online adaptation of language models (LM) has been studied recently to keep the knowledge of LM current [34, 27] . However, trivial auto-regressive fine-tuning the LM, which applies uniform weights to all tokens, often results in suboptimal performance in downstream tasks. This issue stems from the default average negative log-likelihood (NLL) loss, which fails to capture the significance of tokens [25] . To overcome this limitation, [25] proposed Context-aware Meta-learned Loss Scaling (CaMeLS), a strategy that employs meta-learning to adjust token weights for more effective online adaptation. Specifically, they meta train a weight model to reweight the auto-regressive loss during online fine-tuning, aiming to enhance LM performance on downstream question-answering tasks. A comprehensive task description and the mathematical formulation of the objectives are detailed in Appendix E.2.The trained weight model is subsequently fine-tuned on unseen online documents and evaluated on corresponding question-answering tasks. In [25] , RGU is utilized for meta gradient approximation.To mitigate the non-constant memory issue associated with RGU, a DistilGPT2 model [59] is chosen as the surrogate base model for training the weight model, instead of larger models typically employed for online adaptation. Additionally, a very limited unrolled depth of 6 is utilized within a 40 GiB GPU memory budget. In our experiments, since (FG) 2 U has circumvented the non-constant memory issue associated with RGU, we are able to increase the unrolled depth and upscale the base model for training the weight model. Empirical evaluations are conducted on two datasets, StreamingQA [36] and SQuAD-Seq [54] .Model (# params) Method StreamingQA SQuAD-Seq EM (↑) F1 (↑) EM (↑) F1 (↑) DistilGPT2 (82M)CaMeLS + RGU [25, 66] 1.62 Table 2 : Comparison of the online adaptation performance. The reported evaluation metrics include the exact match (EM) and F1 scores. For vanilla CaMeLS [25] , RGU is conducted with unrolled depth 6, using DistilGPT2 as the base model. We present both the results reported by [66] and those from our implementation (denoted as impl.). For CaMeLS + (FG) 2 U, we select unrolled depths from {24, 48}, and the base model from {DistilGPT2, GPT2}. We report the results for the combination that yields the best F1 score. Additional details and ablation studies are documented in Appendix G.1.Firstly, we increased the unrolled depth while maintaining the base model as a DistilGPT2. We plotted the F1 scores and GPU memory usages for RGU with unrolled depths of {1, 2, 4, 6} and (FG) 2 U with unrolled depths of {24, 48} on StreamingQA in Figure 1 (Bottom Left). The performance of the weight model is positively correlated with the unrolled depth, substantiating the benefits of training with larger unrolled depths. The non-constant memory issue associated with RGU can be observed when the unrolled depth increases, while (FG) 2 U maintains constant memory even with large unrolled depth. Subsequently, we endeavored to upscale the base model to GPT2 to reduce the disparity between training and evaluation. The performances are summarized in EQUATIONwhere x denotes the space-time coordinate, Ψ denotes a bounded domain with boundary, u : [0, T ] × Ψ → R denotes the latent solution, u t represents the first-order derivative of u with respect to t, and N is a general differential operator parameterized by ϕ, acting on Ψ. This setup encompasses a broad spectrum of problems in physics. For example, the one-dimensional Burgers' equation is defined by Figure 2 : Left: Comparison of efficiency between the PINN solver and the numerical solver. We evaluated Adam [29] and SGD as the inner optimizers for the PINN solver, with steps ranging from 100 to 50,000. The results demonstrate that the numerical solver is significantly more efficient. Right: Comparison of relative L2 errors in the prediction of ϕ and u.N [u; ϕ] = µuu x -νu xx , where ϕ = (µ, ν) ∈ R 2 ,ϵ ϕ = ∥ϕ pred -ϕ∥ 2 /∥ϕ∥ 2 , ϵ u = ∥u pred -u∥ 2 /∥u∥ 2 .The problem of data-driven discovery of PDEs [52] can be framed as follows: given a set of scattered observations of the latent solution u(x), what are the parameters most accurately describing the observed data? The problem can be formulated as a PDE-constrained optimization problem (PDECO):EQUATIONwhere D = {(x i , u i )} 1:k denotes the observed data. In cases where the closed-form solutions of the nonlinear PDEs are intractable, parametric solutions u θ are used to approximate the latent solution u for given ϕ. The PDECO in (19) is then reformulated into a bi-level optimization problem:EQUATIONEmploying gradient-based PDE solvers, such as physics-informed neural networks (PINN) [52] , facilitates the direct application of (FG) 2 U. However, as demonstrated in Figure 2 (Left), the accuracy and efficiency of PINNs fall short of the rigorous demands of scientific computing. This limitation has prompted us to integrate faster and more accurate traditional solvers like the spectral method [1] (see also Appendix E.3.4) to tackle the inner problem. Given these solvers are non-differentiable, we employ (FG) 2 U-ZO, the zeroth-order variant of (FG) 2 U introduced in Section 3.2, to solve the problem.We conduct experiments on three non-linear PDEs: Burgers, Allen-Cahn, and KdV, with a more detailed task description available in Appendix E.3. The results are summarized in Figure 2 (Right). We can observe that the combination of (FG) 2 U-ZO and the numerical solver significantly outperforms (FG) 2 U and the PINN solver, in terms of both the prediction on ϕ and u. The implementation details are documented in Appendix F.3.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: To the best of our knowledge, all the information needed to reproduce the main experimental results are included in Section 4, Appendix E, and Appendix F. Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.• While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully. (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
conclusion: 
 In this work, we propose a novel algorithm Forward Gradient Unrolling with Forward Gradient, abbreviated as (FG) 2 U, designed to tackle the challenges associated with large-scale bi-level optimization. We conduct a convergence analysis of (FG) 2 U, perform extensive comparisons with existing methods, and provide detailed discussions on its practical applications. Additionally, we undertake an empirical evaluation across a series of large-scale bi-level optimization tasks. Our findings indicate that (FG) 2 U effectively complements existing bi-level optimization algorithms, addressing gaps in large-scale bi-level optimization scenarios.Limitations and future works. The experiments conducted in this paper are of relatively small scale, with the largest inner model being a GPT-2 model. We look forward to validating its effectiveness on larger-scale bi-level optimization tasks. Additionally, the application of black-box bi-level optimization and the potential of (FG) 2 U-ZO remain underexplored, considering the prevalent blackbox interaction between users and models today. We hope our work will inspire further development of large-scale bi-level optimization algorithms and their application in corresponding scenarios. Furthermore, we have not specifically addressed the efficiency issues inherited by (FG) 2 U from the forward gradient method. Enhancing the efficiency of (FG) 2 U while maintaining its gradient estimation accuracy will be an important direction for future research.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [No] Justification: The data used in this paper is public. We will release the code on the acceptance of the paper.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We have reported the error bars in Data Condensation experiments and Data-driven Discovery of PDE experiments. For Meta Learning Online Adaptation of LM experiments, we adhere to the standard evaluation protocol as employed in [25, 66] , specifically, using the identical evaluation script with the same random seed, to ensure a fair comparison, hence, without error bars. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We provide the information of computer resources in Appendix F.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA]• For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: We do not have plans to release any new assets or models in conjunction with this work. The code, which will be released upon acceptance, will be well documented, and the documentation will also be made publicly available.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Algorithm 1 (FG) 2 U: Forward Gradient Unrolling with Forward GradientRequire: Initial inner parameters θ 0 , initial meta parameter ϕ 0 , random direction distribution p, number of random directions b, total meta steps K, meta update mappings Ψ 1:K . 1: θ ← θ 0 , ϕ ← ϕ 0 2: for k = 1, . . . , K do for t = 1, . . . , T do 7: θ ← Ω t (θ, ϕ), A ← ∂Ωt(θ,ϕ) ∂θ , B ← ∂Ωt(θ,ϕ)w i ← ∂f (θ,ϕ) ∂θ y i + ∂f (θ,ϕ) ∂ϕ v i 14:end for 15:ϕ ← Ψ k (ϕ, 1 b b i=1 w i v T i ) 16: end for 17: return ϕ B Extended Discussion on Bi-level Optimization B.Another idea for computing the implicit gradient is to utilize the implicit function theorem (IFT) [30] . Suppose that the inner optimality ∇ θ g(θT (ϕ), ϕ) ≈ 0 is approximately achieved by sufficient inner optimization steps. If g is second-order differentiable, by applying the implicit function theorem and taking the first-order derivative of ϕ,EQUATIONThen, if the Hessian is further assumed to be invertible, the meta gradient can be approximated asEQUATIONEQUATIONThe main challenge lies in the computation of the inverse Hessian matrix H -1 , which is intractable when θ is of high dimensionality. Fortunately, several iterative inverse Hessian vector product (ihvp) approximators requiring only Hessian vector product (hvp) and O(M ) space can be employed to produce dH -1 for approximating dH -1 , based on Conjugate Gradient [48, 53] , Neumann Series [19, 28] and low-rank approximation [64, 24] .Neumann Series. The inverse Hessian vector product can be approximated with a truncated sum of Neumann series [19, 28] ,EQUATION)where α is a hyperparameter to ensure the convergence, and K is the number of truncated steps. Compared to other IF-based methods, the Neumann Series has demonstrated good empirical performance and stability [19] , and its stochastic variant has been well studied [28] .Weakness (IF): Approximation Errors. The errors of IF emanate from two distinct sources Firstly, IFT presupposes that the Karush-Kuhn-Tucker (KKT) conditions of the inner problem are satisfied, leading to an approximation error in (23) when iterative approximations of the inner solutions are used. Secondly, the singular nature of the Hessian within neural network training [57] leads to costly and unstable inverse Hessian approximation in practical applications, with a heavy reliance on engineering efforts [9] . More formally, recallEQUATIONThe approximation error can be decomposed intoEQUATION)To reduce the computational cost, Hessian-free approaches [76, 75, 9] propose approximating the Hessian as an identity matrix, incorporating additional assumptions about the inner model and objective.EQUATIONwhere α > 0 is a hyperparameter to control the magnitude. However, these numerous assumptions often diverge from practical scenarios, resulting in significant approximation errors and consequently inducing suboptimal outcomes. The VF-based methodology [39, 37, 61, 33] considers an equivalent reformulation of the original optimization problem as outlined in (2):EQUATIONThis reformulation casts the standard bi-level optimization challenge into a constrained single-level optimization framework. VF-based methods circumvent the need for second-order computations and have demonstrated near-optimal complexity, comparable to second-order methodologies in deterministic settings, as reported in [8] .Weakness (VF): stochastic optimization. VF-based strategies have yet to gain widespread acceptance in practical ML applications. This limited adoption is primarily attributed to the challenges these methods face in addressing large-scale stochastic problems, where the complexity significantly impedes their performance [73] .In this section, we detail the proofs of Lemma 3.1 and Theorem 3.4. Our approach to proving Theorem 3.4 follows a similar high-level approach as [19, 28, 60] with some important distinctions. Initially, in Lemma B.2, we extend the smoothness properties of the objective functions f and g, and the transition functions Ω, to the T -th iteration lower-level parameter θT . Following this, Lemma B.3 establishes the smoothness of the meta-learning objective f (θT , ϕ), incorporating results from the inner-loop computations. Building on the demonstrated smoothness of the meta objective function and the variance of the forward gradient method (as shown in Lemma 3.1), we then validate the convergence properties of Algorithm 1.The novelty in our proof of Theorem 3.4 lies in two primary aspects. Firstly, our analysis does not presume that the lower-level optimization yields an optimal solution θ * ; instead, it more realistically assumes the use of θT , which is derived from a finite number of iterations. This assumption aligns more closely with the computational constraints encountered in real-world scenarios. Secondly, our convergence analysis explicitly accounts for the variance of our unbiased gradient estimator, achieving a convergence rate of O(ϵ -1 ρ -1 ). This demonstrates that utilizing the forward gradient method, while significantly reducing memory requirements, does not adversely affect the algorithm's convergence rate, underscoring the practical viability and efficiency of our approach even with memory constraints.In Appendix C.3, we discuss how to extend the convergence of optimization problem (2) into (1), with additional assumptions.For convenience, the lemma is restated as follows. Lemma 3.1. For any ϕ ∈ Φ, the gradient estimation with forward gradient method:∇h(ϕ) = 1 b b i=1 ∇h(ϕ)vivi ⊤ ,where vi ∼ Unif({-1, 1} N ), and b denotes the sample size, satifiesE∥ ∇h(ϕ) -∇h(ϕ)∥ 2 = 1 ρ ∥∇h(ϕ)∥ 2 ,where ρ := b N -1 ∈ (0, 1] as b is selected from 1, . . . , N -1.Proof. We start by computing the variance of one-sample estimation, ∇h(ϕ) = ∇h(ϕ)vv ⊤ . Since E[vv ⊤ ] = I, we know that E[ ∇h(ϕ)] = ∇h(ϕ). Consequently,EQUATION)Since v is an N -dimensional Rademacher random variable, we have E∥∇h(ϕ)v∥ 2 = ∥∇h(ϕ)∥ 2 andEQUATION) = E∥∇h(ϕ)vv ⊤ ∥ 2 -∥∇h(ϕ)∥ 2 = (N -1)∥∇h(ϕ)∥ 2 .For the multi-sample estimation ∇h(ϕ) = 1 b b i=1 ∇h(ϕ)vivi ⊤ . Sincevi are i.i.d. sampled, and E[ ∇h(ϕ)] = ∇h(ϕ), we haveE∥ ∇h(ϕ) -E ∇h(ϕ)∥ 2 = E ∇h(ϕ) 1 b b i=1 (vivi ⊤ -I) 2 = 1 b 2 b i=1 E∥∇h(ϕ)(vivi ⊤ -I)∥ 2 = N -1 b ∥∇h(ϕ)∥ 2 .To prove our main result (Theorem 3.4), we first establish useful smoothness properties of the hyperparameter learned from solving the lower-level optimization problem. Subsequently, we establish the smoothness of the meta objective function examined at the approximated lower-level parameters in Lemma B.3.Regarding the lower-level parameter θ(ϕ), we present the following lemma, which is based on Assumption 3. ∈ Φ, ∥θT (ϕ) -θT (ϕ ′ )∥ ≤ C Z ∥ϕ -ϕ ′ ∥, ∥∇θT (ϕ) -∇θT (ϕ ′ )∥ ≤ L Z ∥ϕ -ϕ ′ ∥,whereC Z = C T +2 Ω -C Ω C Ω -1 and L Z = L Ω C T Ω + C T +2 Ω T C Ω -1 - C T Ω -1 (C Ω -1) 2 .Proof. We start with the proof of Lipshitz continuity. For any pair of ϕ, ϕ ′ ∈ Φ, using (2) and Assumption 3.3, we haveEQUATIONApplying (31) recursively over s = 1, . . . , t givesEQUATIONwhere the last inequality holds from the fact that θ0 = Ω0 as well as Assumption 3.3, and the subsequent equality follows from the geometric series summation formula.Therefore, θt(ϕ) is C Z (t)-Lipchitz, where C Z (t) :=C t+2 Ω -C Ω C Ω -1 . Substituting t = T , we get C Z = C Z (T ) = C T +2 Ω -C Ω C Ω -1 .We now proceed with the proof of L Z -smoothness. For simplicity of notation, we follow (4) and denoteZt(ϕ) = ∇θt(ϕ); At(ϕ) = ∂Ωt(θt-1(ϕ), ϕ) ∂θt-1 ; Bt(ϕ) = ∂Ωt(θt-1(ϕ), ϕ) ∂ϕ .Subsequently, considering the update rule Zt(ϕ) = At(ϕ)Zt-1(ϕ) + Bt(ϕ) of Forward Gradient Unrolling (4), we haveEQUATIONwhere the last inequality follows from Assumption 3.3 that Ω0(ϕ) and Ω1:T (ψ) are C Ω -Lipshitz and L Ωsmooth, and the previously proved result that θt(ϕ) is C Z (t)-Lipchitz.Noting that Z0(ϕ) = ∇θ0(ϕ) = ∇Ω0(ϕ), applying (33) recursively over t = 1, . . . , T gives∥∇θT (ϕ) -∇θT (ϕ ′ )∥ ≤ C T Ω ∥Z0(ϕ) -Z0(ϕ')∥ + T t=1 C T -t Ω (C Z (t) + 1)L Ω ∥ϕ -ϕ ′ ∥ = C T Ω ∥∇Ω0(ϕ) -∇Ω0(ϕ ′ )∥ + C T Ω T t=1 C Z (t) + 1 C t Ω L Ω ∥ϕ -ϕ ′ ∥ ≤ L Ω C T Ω ∥ϕ -ϕ ′ ∥ + C T Ω T t=1 C t+2 Ω -1 C t Ω (C Ω -1) L Ω ∥ϕ -ϕ ′ ∥ = L Ω C T Ω + C T Ω T t=1 C 2 Ω C Ω -1 - C T Ω C Ω -1 T t=1 1 C t Ω ∥ϕ -ϕ ′ ∥ = L Ω   C T Ω + C T +2 Ω T C Ω -1 - C T Ω C Ω -1 1 C Ω (1 -1 C T Ω ) 1 -1 C Ω   ∥ϕ -ϕ ′ ∥ = L Ω C T Ω + C T +2 Ω T C Ω -1 - C T Ω -1 (C Ω -1) 2 ∥ϕ -ϕ ′ ∥,where the third line follows from Assumption 3.3 that Ω0(ϕ) is L Ω -smooth and the choice C Z (t) =C t+2 Ω -C Ω C Ω -1 (which gives C Z (t) + 1 = C t+2 Ω -1C Ω -1 ), and the fifth line again uses the geometric series summation formula.Hence, θT (ϕ) is L Z -smooth with L Z = L Ω C T Ω + C T +2 Ω T C Ω -1 - C T Ω -1 (C Ω -1) 2 .Next, we provide a lemma establishing that the upper-level objective f , evaluated at the learned parameter (θT (ϕ), ϕ), also adheres to certain smoothness properties. Lemma B.3. Define h(ϕ) := f (θT (ϕ), ϕ). Under Assumptions 3.2 and 3.3, h(ϕ) is L h -smooth, i.e., for any ϕ, ϕ ′ ∈ Φ, ∥∇h(ϕ) -∇h(ϕ ′ )∥ ≤ L h ∥ϕ -ϕ ′ ∥, where L h = (C Z + 1) 2 L + CL Z , with C Z and L Z defined in Lemma B.2.Proof. For simplicity of notation, we follow (5) and denoteZt(ϕ) = ∇θt(ϕ); cT (ϕ) = ∂f (θT (ϕ), ϕ) ∂ϕ ; dT (ϕ) = ∂f (θT (ϕ), ϕ) ∂θT .For any ϕ, ϕ ′ ∈ Φ, following a similar proof as Lemma B.2, we have ∥∇h(ϕ) -∇h(ϕ ′ )∥= ∥dT (ϕ)ZT (ϕEQUATIONSubsequently, we deduce that(34) ≤ C Z ∥dT (ϕ) -dT (ϕ ′ )∥ + CL Z ∥ϕ -ϕ ′ ∥ + ∥cT (ϕ) -cT (ϕ ′ )∥ ≤ C Z ∂f (θT (ϕ), ϕ) ∂θT - ∂f (θT (ϕ ′ ), ϕ) ∂θT + ∂f (θT (ϕ ′ ), ϕ) ∂θT - ∂f (θT (ϕ ′ ), ϕ ′ ) ∂θT + ∂f (θT (ϕ), ϕ) ∂ϕ - ∂f (θT (ϕ ′ ), ϕ) ∂ϕ + ∂f (θT (ϕ ′ ), ϕ) ∂ϕ - ∂f (θT (ϕ ′ ), ϕ ′ ) ∂ϕ ′ + CL Z ∥ϕ -ϕ ′ ∥ ≤ C Z L∥θT (ϕ) -θT (ϕ ′ )∥ + C Z L∥ϕ -ϕ ′ ∥ + L∥θT (ϕ) -θT (ϕ ′ )∥ + L∥ϕ -ϕ ′ ∥ + CL Z ∥ϕ -ϕ ′ ∥ ≤ C Z LC Z ∥ϕ -ϕ ′ ∥ + C Z L∥ϕ -ϕ ′ ∥ + LC Z ∥ϕ -ϕ ′ ∥ + L∥ϕ -ϕ ′ ∥ + CL Z ∥ϕ -ϕ ′ ∥ = [(C Z + 1) 2 L + CL Z ]∥ϕ -ϕ ′ ∥,where the first, third, and fourth lines all follow directly from Lemma B.2 and Assumption 3.3 (recall that the latter states that f is L-Lipschitz and C-smooth).Therefore, h(ϕ) is L h -smooth with L h = (C Z + 1) 2 L + CL Z .Now based on the aforementioned lemmas, we put forward the proof of our main theorem: the convergence analysis for our bilevel optimization method (FG) 2 U. Theorem 3.4 (Convergence). Suppose that Asumption 3.2 and Assumption 3.3 hold. Setting the learning rate β = ρ (ρ+1)L h for gradient descent over the hyperparameter ϕ, we haveEQUATION)Proof. We haveEQUATION)where the second line is a well-known inequality for smooth functions with the L h -smoothness itself following from Lemma B.3, and the third line uses the gradient descent rule ϕ k+1 = ϕ k -β ∇h(ϕ k ).By Lemma 3.1 and the fact that ∇h is unbiased (see the proof of Lemma 3.1), we know thatE[⟨∇h(ϕ k ), ∇h(ϕ k )⟩|ϕ k ] = ∥∇h(ϕ k )∥ 2 ; E[∥∇h(ϕ k ) -∇h(ϕ k )∥ 2 |ϕ k ] = 1 ρ ∥∇h(ϕ k )∥ 2 .Therefore, taking the conditional expectation E[ • | ϕ k ] over (36) givesEQUATION)Furthermore, taking the full expectation and telescoping (37) over k form 0 to K -1 yieldsEQUATIONChoosing β = ρ (ρ+1)L h , we haveEQUATIONHence, Algorithm 1 requires O(ϵ -1 ρ -1 ) steps to attain an ϵ-accurate stationary point.Convergence of Problem (1). To extend the convergence of optimization problem (2) into (1), we need to assume that the lower-level objective function g is strongly convex w.r.t. θ as commonly done by previous works [60, 28] . From the strong convexity and first-order smoothness (Assumption 3.2) of g, we have 1) the zeroth and first-order smoothness of θ * (ϕ); 2) ∥θT (ϕ ′ ) -θ * (ϕ ′ )∥ → 0 as T → +∞. Then the inequalityEQUATIONimplies the the zeroth and first-order smoothness of θT (ϕ). Following the same line of proof as presented in our paper, we derive the smoothness of f (θ * (ϕ), ϕ) and f (θT (ϕ), ϕ), and subsequently the convergence of either problem (1) or (2).However, as discussed in Section 2, given that the scope of this paper is large-scale BO, the inner optimization typically involves deep neural networks. Therefore, the optimal parameters are not explicitly accessible and can only be estimated through iterative procedures. Most related works [9, 60, 28] are implicitly or explicitly solving (2) instead of (1). Additionally, it is important to acknowledge that achieving strong convexity is often unfeasible in practical applications. Consequently, we focus on (2), aiming to present a more practical convergence theory that proves the effectiveness of our method.In this section, we give a more detailed introduction to zeroth-order (ZO) derivative estimators. These estimators are pivotal in scenarios where the computation of exact derivatives is either infeasible due to memory constraints or computationally prohibitive. Apart from the forward gradient method employed in (FG) 2 U, randomized smoothing (RS) is another widely-used derivative estimator, both in Reinforcement Learning [22, 32] and Large Language Models [18, 47, 74] .For a function F : R n → R, gradient estimation via RS can be mathematically formulated as:EQUATION)where b is the number of random samples, ϵ is the smoothing parameter, vi are samples drawn from a standard Gaussian distribution. Regarding the accuracy of estimation, it has been shown in [13, 42] that the variance of RS is roughly in the order of O(N/b), which is the same as FG as proved in Lemma 3.1.RS stands out particularly in its ability to estimate gradients of functions evaluated through black-box systems, where internal operations are inaccessible or highly complex. This characteristic makes RS exceptionally valuable in practical applications such as adversarial robustness and black-box optimization, where obtaining direct gradients might not be possible. Another advantage of RS is its robustness against noise and discontinuities in the function landscape. Unlike deterministic methods, the stochastic nature of RS allows it to approximate the gradient over a smoothed version of the function, providing stability in scenarios where slight perturbations can lead to substantial changes in the output.While RS provides robust gradient estimates across various scenarios, it is critical to recognize that RS inherently introduces bias if the expectation is not computed during inference. In many CV and NLP applications, the computational expense of Monte Carlo sampling at the evaluation stage is prohibitive, leading to a biased estimation when using RS. However, in the context of inverse PDE problems, where the inner-loop solvers are non-differentiable numerical solvers, we employ RS as a zeroth-order derivative estimator.In the era of rapid advancement in machine learning, a multitude of foundation models [11, 6, 55] has benefited from training on large-scale datasets, exhibiting formidable performance that models trained on small-scale data cannot match. However, the exponential growth of data also presents challenges: (1) Models updated with only new data are prone to catastrophic forgetting [20] while retaining all historical data for subsequent training imposes significant storage and computational burdens. ( 2) Applications within the realm of meta-learning, such as hyperparameter tuning [46, 43] and neural architecture search [78, 38] , necessitate multiple training iterations over datasets. The computational cost of these operations scales dramatically with the size of the datasets, posing a bottleneck for efficiency and scalability. (3) The widespread dissemination and utilization of datasets have raised significant concerns regarding privacy and copyright [12] .To overcome the challenges posed by large-scale datasets, a line of work known as data condensation [68, 72] has been proposed, with the idea to generate a compact, synthesized dataset, designed to elicit similar behaviors in machine learning models as those trained with the original, massive dataset. The objectives of the mainstream principles [72] designed for data condensation can be naturally formulated as a bi-level optimization problem. We focus on the best-known principle performance matching [72] on classification task, which can be formulated as, min Do L(θT ; Do), where θt = θt-1 -η∇L(θt-1; Dc), t = 1, . . . , T,We conduct our experiments to condense the following image datasets:• MNIST [35] : a handwritten digits dataset containing 60, 000 training images and 10, 000 testing images with the size of 28 × 28 from 10 categories.• CIFAR 10/100 [31] : colored natural images datasets contraining 50, 000 training images and 10, 000 testing images from 10/100 categories, respectively.The scale of the condensed dataset will fundamentally impact the results. Therefore, we consider different scales for each dataset, with images per class set to 1, 10, and 50. The condensed dataset will be used to train random initialized models, and evaluated on a test dataset.The online adaptation of language models (LM) [34, 27] has been studied recently to keep the knowledge of LM updated to date. However, trivial auto-regressive fine-tuning the LM with uniform weights for all tokens results in poor performance in downstream tasks, as the default average negative log-likelihood (NLL) loss does not accurately reflect the importance of tokens [25] . To address the issue, [25] proposed Context-aware Meta-learned Loss Scaling (CaMeLS) to meta-learning the weights of tokens for effective online adaption. More formally, let θ denote the parameter of the base model for adaptation, ϕ denote the parameter of a parametric weight model to assign weights for each token, the meta-learning online adaption of LM can be formulated as the following bi-level optimization,EQUATIONWe follow the setting studied by [25] , where the downstream task is question-answering. The meta-objective consists of a question-answering term measuring the performance gained from adaptation, and a locality term that prevents the updated base model parameters from excessively changing the base model's behavior. Let DQA denotes the question-answering dataset, D loc denotes the locality dataset, and c ∈ R + denotes the weight of the locality term, then the meta objective is formally defined asEQUATIONThe inner objective is defined as a weighted NLL loss, where the weights are determined by the weight model The trained weight model is then fine-tuned on unseen online documents and evaluated on corresponding question-answering tasks. The nonlinear viscous Burgers equation is a pivotal partial differential equation arising in diverse domains of applied mathematics such as fluid mechanics, nonlinear acoustics, and traffic flow. This equation can be deduced from the Navier-Stokes equations for the velocity field by omitting the pressure gradient term. In our experiment, the equation along with Dirichlet boundary conditions, is expressed as follows:EQUATIONwith actual viscosity ν = 0.01 π ≈ 0.0031831. For PINN, following [44] , we enforced the initial condition into the output by choosing a surrogate model of the solution asû(x) = (1 -exp(-t)) NN(x; θ) -sin(πx),where NN(x; θ) is a neural network.The Allen-Cahn equation is a reaction-diffusion equation of mathematical physics describing the process of phase separation in multi-component alloy systems, including order-disorder transitions. In our experiment, it is expressed as follows:EQUATIONwith the actual diffusion coefficient ν = 0.001. For PINN, we enforced the initial condition into the output by choosing a surrogate model of the solution asû(x) = (1 -exp(-t)) NN(x; θ) + x 2 cos(πx),where NN(x; θ) is a neural network.The Korteweg-de Vries (KdV) equation serves as a mathematical model for waves on shallow water surfaces. This equation is distinguished as a prototypical example of an integrable PDE. It is characterized by features typical of integrable systems, including a plethora of explicit solutions, notably soliton solutions, and an infinite number of conserved quantities. These properties are particularly noteworthy given the inherent nonlinearity of the equation, which generally complicates the solvability of PDEs. In specific, we consider:EQUATIONwith the actual coefficient of dispersion ν equal to 0.0025. For PINN, we enforced the initial condition into the output by choosing a surrogate model of the solution asû(x) = (1 -exp(-t)) NN(x; θ) + cos(πx),where NN(x; θ) is a neural network.Numerical solvers play a critical role in the study and application of PDEs, enabling the simulation and analysis of complex physical phenomena that cannot be addressed analytically [51] . These solvers convert PDEs into a form that can be handled computationally, typically by discretizing the domain into a finite set of points or elements and approximating the derivatives. Conventional numerical methods include finite difference methods, finite element methods, and spectral methods [1] .Among the various numerical methods for solving PDEs, spectral methods stand out for their ability to deliver highly accurate solutions, particularly for problems with smooth solutions [21, 7] . Spectral methods involve representing the solution to a PDE as a sum of basis functions, such as trigonometric polynomials, which are globally defined over the domain. This approach contrasts with finite difference or finite element methods, where the solution is localized to the grid points or elements. In this paper, we mainly adopt spectral methods, as we focus on the Burgers, Allen-Cahn, and KdV equations. All these three equations can be efficiently and accurately resolved by spectral techniques.We conducted our experiments following the standard data condensation setting established by [68, 77, 67] . The condensation and evaluation are both performed on a depth-3 convolutional neural network [58] . The hyperparameters we used for (FG) We adhered to the standard settings of CaMeLS [25] and adapted their official code for our implementation. The only modification made was replacing the meta gradient approximation module with (FG) 2 U. It is important to note that the base models used for meta-learning were initially pre-trained on a split QA-paired set. While the official codebase provided the script for pretraining, it did not include the exact base model (weights) they used. We executed the official script to generate the pre-trained base models and observed that meta-learning performance is sensitive to the choice of base models. For a fair comparison, we reported both the results from [66] (where CaMeLS [25] presented performance improvements over baselines using bar plots without specific metric values) and the results with our best custom pre-trained base models. Following the two-phase training paradigm introduced in Section 3.2, we performed training of (FG) 2 U on RGU (DistilGPT2, unrolled depth 6) results. The hyperparameters we used for (FG) 2 U are summarized in Appendix F.2, while all remaining hyperparameters were kept the same as in [25] . All experiments are conducted on one NVIDIA A100 GPU (80G • The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.Justification: To the best of our knowledge, all necessary information has been included in Section 4, Appendix E, and Appendix F.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material.• The answer NA means that the paper does not include experiments. • The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.Justification: The research conducted in the paper conforms, in every respect, with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


Decision: This paper introduces a Forward Gradient Unrolling with Forward Gradient, abbreviated as (FG)2U, achieving an unbiased stochastic approximation of the meta gradient. All the reviewers agree to accept this paper based on reasons of (1) simple and easy to implement, (2) well linked theoretical contribution, (3) comprehensive experiments, etc. I mostly agree.

Some concerns are also raised: (1) Theorem 3.4 indicates that the convergence rate depends linearly on N and states the various useful implications about it; (2) the computational cost takes the order of O(N) and it is suggested to add discussion on the zero-order BO.

Thus, I recommend an accept subject to addressing the mentioned issues in the final version

Reviews: review: paper_summary: paper_summary: This paper presents a novel gradient unrolling algorithm for bi-level optimization. The authors highlight that existing methods for calculating meta gradients in the literature are not memory efficient. They propose a sampling-based method, (FG)^2U, to estimate the meta gradient. This approach approximates the meta gradient by multiplying it with a random rank-1 matrix, thereby simplifying the Forward Gradient Unrolling (FGU) scheme. The paper includes discussions on sampling efficiency, convergence analysis, and numerical experiments for (FG)^2U.
strength_weakness: strength_weakness: A major concern is the theoretical dependence on the parameter dimension $N$. Theorem 3.4 indicates that the convergence rate depends linearly on N. It means that
- With a fixed batch size, convergence on large-scale applications would be slow (vanilla (FG)^2U without additional techniques).
- With an $O(N)$ batch size, convergence is satisfactory, but the memory cost grows as $O(MN)$, similar to FGU.
- While the authors mention that (FG)^2U allows parallelization to mitigate computational overhead, it seems that the calculations of FGU (4) also permit parallelization, correct? Specifically, $d \theta / d \phi$ can be computed in parallel as $d \theta[1] / d \phi, d \theta[2] / d \phi, d \theta[3] / d \phi, \cdots$
questions: questions: Regarding the formula for data condensation (17) and related formulas in the appendix: should it be minimizing over ${\mathcal{D}_c}$ instead of minimizing over ${\mathcal{D}_o}$?
limitations: limitations: The paper appears to have no potential negative societal impact. The authors discussed the limitations in the Appendix H.
paper_summary: paper_summary: The paper introduces a method called Forward Gradient Unrolling with Forward Gradient, abbreviated as (FG)²U, which is designed to address the large memory requirements of forward method in bi-level optimization in large-scale machine learning model
strength_weakness: strength_weakness: 1. The proposed method introduces additional computational complexity. Can the author give an analysis of the complexity?
2. The convergences analysis is for the problem (2) not the original problem (1). 
3. More large scale datasets are needed.
questions: questions: see weakness
limitations: limitations: NA
paper_summary: paper_summary: This paper studies scalable bi-level optimization problems and points out the limitations of most traditional methods. Among all, to mitigate the high memory cost issue of the GU method, it proposes (FG)2U such that the space consumption is reduced to $\mathcal{O}(M)$ from $\mathcal{O}(MN)$. Convergence analysis is provided and extensive experiments show the benefit of the algorithm.
strength_weakness: strength_weakness: 1.	The memory cost is reduced but the computational cost seems to increase significantly. The computational cost takes the order of $\mathcal{O}(KTN)$ since $b=\mathcal{O}(N)$. This feels unrealistic as well in the large-scale application when $N$ is large. 
2.	Based on the previous question, I checked the choice of $b$ in the experiment, which is small. How do authors tune this parameter? Is there any ablation study on it?
3.	In terms of the discussion in Appendix B.2, IF methods also require $\mathcal{O}(M)$ space using some approximation tricks. There indeed exist approximation errors in IF methods. However, normally, the approximation errors can be controlled to be very small based on the hyperparameter selections. Thus, IF methods seem suitable for scalable bi-level optimization problems to some extent. Do authors have time to include IF methods in Table 2?
4.	The algorithm is deterministic without consideration of the batch data. Can authors provide some insights in stochastic settings? For example, can the unbiasedness of the hypergradient estimator still be satisfied, or is there any requirement for the batch size? 
5.	Can authors provide real space consumption in the comparison (Table 1, Table 2)?
questions: questions: See weakness part.
limitations: limitations: See weakness part.

Comments: title: Thanks again for your valuable feedback
Thanks again for your valuable feedback
Response to authors
Thanks again
Thanks for the response
Response to author
comment: Thank the authors for their further clarifications. 

Best,

Reviewer
Thank you for your response and for the recognition of our work.

Regarding the convergence guarantee dependent on the problem dimension, we acknowledge that further improvement of the theoretical results would require assumptions that may be impractical. In the field of backpropagation-free optimization (FG/ZO), the gap between dimension-dependent theoretical results and the significantly positive empirical outcomes in large-scale cases remains an open question. We hope that future research will help bridge this gap.

Regarding the parallelization plan, the essential information we intended to convey in our rebuttal response is that "a special parallelization plan is covered by $($FG$)^2$U as a universal framework." We acknowledge that we may have overstated the advantages of $($FG$)^2$U over FGU and RGU in terms of memory management, as pointed out by the reviewer. We will exercise caution with related statements in our revised paper.

We would like to once again extend our sincere thanks for your valuable feedback. We remain open to answering any further questions during the remaining time of the discussion period.
Thank you for your response. We regret to hear that some concerns remain.

Regarding the convergence guarantee dependent on the problem dimension, we acknowledge that further improvement of the theoretical results would require assumptions that may be impractical. In the field of backpropagation-free optimization (FG/ZO), the gap between dimension-dependent theoretical results and the significantly positive empirical outcomes in large-scale cases remains an open question. We hope that future research will help bridge this gap.

Regarding BO + ZO, we followed the reviewer’s suggestions and identified several related works [1][2][3]. Given the limited time remaining in the discussion period, we are unable to conduct a comprehensive empirical study of these methods. However, we offer some preliminary comments here and will seriously consider incorporating comparisons in our revised paper.

1. [1][3] propose utilizing zeroth-order Hessian/Jacobian approximations for IF-based methods, whereas our work focuses on GU-based methods. It is important to note that zeroth-order approximation cannot eliminate the inherent bias introduced by IF-based methods, and the theoretical guarantees provided in these works are also dimension-dependent.

2. [2] employs zeroth-order optimization in a GU manner. However, [2] is limited to Neural Architecture Search rather than universal BO. Additionally, the inner problem considered in [2] is differentiable (white box), while our exploration of $($FG$)^2$U-ZO is in the more challenging black-box setting.

We would like to once again extend our sincere thanks for your valuable feedback. We remain open to answering any further questions during the remaining time of the discussion period.

$ $

[1] On the Convergence Theory for Hessian-Free Bilevel Algorithms, https://arxiv.org/abs/2110.07004

[2] ZARTS: On Zero-order Optimization for Neural Architecture Search, https://arxiv.org/abs/2110.04743

[3] Fully Zeroth-Order Bilevel Programming via Gaussian Smoothing, https://arxiv.org/pdf/2404.00158
I greatly appreciate the efforts made by the authors during the rebuttal phase, including their responses to my questions and the inclusion of additional experimental results.

However, my primary theoretical concern, which I raised in the initial review, remains only partially addressed.

The statement "Plan A requires exactly $M$ threads" may not be entirely accurate. It is possible to adaptively allocate the $M$ calculation tasks based on the constraints of the hardware. For instance, if $M=5$, one could opt for 3 threads and distribute the tasks as $[1,1,3]$, rather than being limited to $[1,1,1,1,1]$. This flexibility also casts doubt on the assertion that "the memory cost of $\text{(FG)}^2\text{U}$ can be more manageable than FGU and RGU," as one could potentially implement simple variants of FGU or RGU.

Of course, as the authors have noted, these analyses are merely upper bounds. The actual performance in the experiments presented by the authors appears promising. Therefore, I will keep my positive score.
We are pleased to hear that our response has addressed your concerns. We would like to extend our sincere thanks for your valuable feedback again.
I thank the authors for their response. Some of my concerns have been addressed. However, I am still slightly concerned about the theoretical guarantee on the problem dimension. There are also some works on bilevel optimization using zeroth-order types of methods. Maybe the authors would like to include them and provide a comparison. I keep my current score. 

Best,
Reviewer
Thank you for the rebuttal, the authors have addressed all my concerns. I will increase my score.