Non-asymptotic Analysis of Biased Adaptive Stochastic Approximation

Authors: sobihan  surendran
adeline  fermanian
antoine  godichon-baggioni
sylvain  le corff
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 18:50:22

Abstract: Stochastic Gradient Descent (SGD) with adaptive steps is widely used to train deep neural networks and generative models. Most theoretical results assume that it is possible to obtain unbiased gradient estimators, which is not the case in several recent deep learning and reinforcement learning applications that use Monte Carlo methods.
This paper provides a comprehensive non-asymptotic analysis of SGD with biased gradients and adaptive steps for non-convex smooth functions. Our study incorporates time-dependent bias and emphasizes the importance of controlling the bias of the gradient estimator. 
In particular, we establish that Adagrad, RMSProp, and AMSGRAD, an exponential moving average variant of Adam, with biased gradients, converge to critical points for smooth non-convex functions at a rate similar to existing results in the literature for the unbiased case. Finally, we provide experimental results using Variational Autoenconders (VAE) and applications to several learning frameworks that illustrate our convergence results and show how the effect of bias can be reduced by appropriate hyperparameter tuning.

Main paper: 
 introduction: 
 Stochastic Gradient Descent (SGD) algorithms are standard methods to train statistical models based on deep architectures. Consider a general optimization problem:EQUATIONwhere V is the objective function. Then, Gradient Descent methods produce a sequence of parameter estimates as follows: θ 0 ∈ R d and for all n ∈ N,θ n+1 = θ n -γ n+1 ∇V (θ n ) ,where ∇V denotes the gradient of V and for all n ≥ 1, γ n > 0 is the learning rate. In many cases, it is not possible to compute the exact gradient of the objective function, hence the introduction of vanilla Stochastic Gradient Descent, defined for all n ∈ N by:θ n+1 = θ n -γ n+1 ∇V (θ n ) ,where ∇V (θ n ) is an estimator of ∇V (θ n ). For example, in deep learning, stochasticity emerges with the use of mini-batches. While these algorithms have been extensively studied, both theoretically and practically, see, e.g., [10] , many questions remain open. In particular, most results are based on the case where the estimator ∇V is unbiased. Although this assumption is valid in the case of vanilla SGD, it breaks down in many common applications. For example, zeroth-order methods used to optimize black-box functions [61] in generative adversarial networks [58, 16] have access only to noisy biased realizations of the objective functions.Furthermore, in reinforcement learning algorithms such as Q-learning [42] , policy gradient [5] , and temporal difference learning [8, 52, 18] , gradient estimators are often obtained using a Markov chain with state-dependent transition probability. These estimators are then biased [69, 23] . Other examples of biased gradients can be found in the field of generative modeling with Markov Chain Monte Carlo (MCMC) and Sequential Monte Carlo (SMC) [34, 13] . In particular, the Importance Weighted Autoencoder (IWAE) proposed by [12] , which is an extension of the standard Variational Autoencoder (VAE) [48] , yields biased estimators. Finally, this is also the case in Bilevel Optimization [43, 36, 41] and Conditional Stochastic Optimization [40, 39] .Moreover, in practical applications, vanilla SGD shows difficulties in calibrating the step sequences. Therefore, modern variants of SGD employ adaptive steps that use past stochastic gradients or Hessians to avoid saddle points and deal with ill-conditioned problems. The idea of adaptive steps was first proposed in the online learning literature by [4] and later adopted in stochastic optimization, with the Adagrad algorithm of [27] .In this paper, we give non-asymptotic convergence guarantees for modern variants of SGD where both the estimators are biased and the steps are adaptive. To our knowledge, existing results consider either adaptive steps but unbiased estimators [27, 77, 67, 74, 19] , or biased estimators with non-adaptive steps [70, 44, 2, 22, 21] .More precisely, our contributions are summarized as follows.• We provide convergence guarantees for the Biased Adaptive Stochastic Approximation framework, under weak assumptions on the bias. To the best of our knowledge, these are the first convergence results to incorporate adaptive steps in biased Stochastic Approximation.• In particular, we establish that Adagrad, RMSProp, and AMSGRAD, an exponential moving average variant of Adam, with a biased gradient, converge to a critical point for non-convex smooth functions with a convergence rate of O(log n/ √ n + b n ), where b n is related to the bias at iteration n. However, we achieve an improved linear convergence rate with the Polyak-Łojasiewicz (PL) condition.• Finally, we show how our theoretical results apply to several applications with biased gradients. In particular, we show that our hypotheses hold for Stochastic Bilevel Optimization and Conditional Stochastic Optimization, but also for Self-Normalized Importance Sampling estimators or Coordinate Sampling. We also propose a first non-asymptotic bound on the bias of IWAE, which allows us to illustrate through several experiments the effect of bias on the convergence of the optimization, and to show how this effect can be reduced by an appropriate choice of hyperparameters.Organization of the paper. In Section 2, we introduce the setting of the paper and relevant related works. In Section 3, we present the Adaptive Stochastic Approximation framework and the main assumptions. In Section 4, we present our principal results, i.e., convergence rates for the risk when the PL condition is assumed, and on the gradient norm without this hypothesis. We illustrate our results in Section 5. All proofs are postponed to the appendix.2 Setting and Related Works Stochastic Approximation. Stochastic Approximation (SA) methods go far beyond SGD. They consist of sequential algorithms designed to find the zeros of a function when only noisy observations are available. Indeed, [68] introduced the Stochastic Approximation algorithm as an iterative recursive algorithm to solve the following integration equation:EQUATION)where h is the mean field function, X is a random variable taking values in a measurable space (X, X ), and E π is the expectation under the distribution π. In this context, H θ can be any arbitrary function.If H θ (X) is an unbiased estimator of the gradient of the objective function, then h(θ) = ∇V (θ). As a result, the minimization problem (1) is then equivalent to solving problem (2), and we can note that SGD is a specific instance of SA. SA methods are then defined as follows:θ n+1 = θ n -γ n+1 H θn (X n+1 ) ,where the term H θn (X n+1 ) is the n-th stochastic update, also known as the drift term, and is a potentially biased estimator of ∇V (θ n ). It depends on a random variable X n+1 which takes its values in (X, X ). In machine learning, V typically represents the theoretical risk, θ the model parameters, and X n+1 the data.Adaptive Stochastic Gradient Descent. SGD can be traced back to [68] , and its averaged counterpart was proposed by [66] . The non-asymptotic analysis of SGD in both convex and strong convex cases can be found in [59] . [32] prove the convergence of a random iterate of SGD for nonconvex smooth functions, which was already suggested by the results of [9] . They show that SGD with constant or decreasing stepsize γ n = 1/ √ n converges to a stationary point of a non-convex smooth function V at a rate of O(1/ √ n) where n is the number of iterations.Most adaptive first-order methods, such as Adam [47] , Adadelta [78] , RMSProp [72] , and NADA [25] , are based on the blueprint provided by the Adagrad family of algorithms. The first known work on adaptive steps for non-convex stochastic optimization, in the asymptotic case, was presented by [50] . [74] proved that Adagrad converges to a critical point for non-convex objectives at a rate of O(log n/ √ n) when using a scalar adaptive step. In addition, [79] extended this proof to multidimensional settings. More recently, [19] focused on the convergence rates for Adagrad and Adam. Furthermore, several modified versions of Adam have been proposed, such as AMSGRAD [77] and YOGI [67] .Biased Stochastic Approximation. The asymptotic results of Biased SA have been studied by [70] . The non-asymptotic analysis can be found in the reinforcement learning literature, especially in the context of temporal difference (TD) learning, as explored by [8, 52, 18] . The case of non-convex smooth functions has been studied by [44] . The authors establish convergence results for the mean field function at a rate of O(log n/ √ n + b), where b corresponds to the bias and n to the number of iterations. For strongly convex functions, the convergence of SGD with biased gradients can be found in [2] , specifically addressing the case of Martingale noise with a constant step size. [46, 21] introduce a novel assumption, known as "Expected Smoothness", which is the weakest assumption compared to the existing literature on biased SGD that we extend to cover the adaptive case. The authors provide convergence results in the case of non-convex smooth functions. Convergence results with assumptions on the control of bias and MSE can be found in [56, 22] . Applications of biased gradients can be found in Bilevel Optimization [43, 36, 41] and Conditional Stochastic Optimization [40, 39] . Moreover, biased gradients are also used in various other applications [38, 54, 6, 56] . Finally, [3] studied convergence results of biased gradients with Adagrad in the Markov chain case, focusing on the norm of the gradient of the Moreau envelope while assuming the boundedness of the objective function.Our analysis provides non-asymptotic results in a more general setting, for a wide variety of objective functions and adaptive algorithms and treating both the Martingale and Markov chain cases. 
limitations: 
 Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: In our paper, we discuss the limitations of our results. Furthermore, for each assumption, we also discuss where it can be verified and provide the limitations of these assumptions.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• 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. 
framework: 
 Consider the optimization problem (1) where the objective function V is assumed to be differentiable. In this paper, we focus on the following SA algorithm with adaptive steps: θ 0 ∈ R d and for all n ∈ N,EQUATIONwhere γ n+1 > 0 and A n is a sequence of symmetric and positive definite matrices. In a context of biased gradient estimates, choosingA n = δI d + 1 n + 1 n k=0 H θ k (X k+1 )H θ k (X k+1 ) ⊤ -1/2can be assimilated to the full Adagrad algorithm [27] . However, computing the square root of the inverse becomes expensive in high dimensions, so in practice, Adagrad is often used with diagonal matrices. This approach has been shown to be particularly effective in sparse optimization settings. Denoting by Diag(A) the matrix formed with the diagonal terms of A and setting all other terms to 0, Adagrad with diagonal matrices is defined in our context as:EQUATIONwhere Hn (X 1:n+1 , θ 0:n ) = 1 n + 1 n k=0 H θ k (X k+1 )H θ k (X k+1 ) ⊤ .In RMSProp [72] , Hn (X 1:n+1 , θ 0:n ) in ( 4) is an exponential moving average of the past squared gradients, defined by:Hn (X 1:n+1 , θ 0:n ) = (1 -ρ) n k=0 ρ n-k H θ k (X k+1 )H θ k (X k+1 ) ⊤ ,where ρ is the moving average parameter. Furthermore, when A n is a recursive estimate of the inverse Hessian, it corresponds to the Stochastic Newton algorithm [11] .H1 There exists a constant µ > 0 such that for allθ ∈ R d , 2µ V (θ) -V (θ * ) ≤ ∥∇V (θ)∥ 2 .H1 corresponds to the Polyak-Łojasiewicz condition, which is weaker than strong convexity and remains satisfied even when the function is non-convex. It ensures uniqueness of the minimizer θ * . The PL condition has been extensively studied theoretically [45] and has been verified empirically in many applications, such as over-parameterized deep networks [26] and Linear Quadratic Regulator models [29] .H2 The objective function V is L-smooth. For all (θ, θ ′ ) ∈ R d × R d , ∥∇V (θ) -∇V (θ ′ )∥ ≤ L ∥θ -θ ′ ∥ .This assumption is crucial to obtain our convergence rate and is very common see, e.g., [59, 10] . Under this assumption, for all (θ, θEQUATIONH3 (i) Biased Gradients: There exist two non-increasing positive sequences (λ n ) n≥1 and (r n ) n≥1 such that for all n ∈ N,E ⟨∇V (θ n ) , A n H θn (X n+1 )⟩ ≥ λ n+1 E ∥∇V (θ n )∥ 2 -r n+1 .(ii) Expected Smoothness: there exists a non-increasing non-negative sequence (σ 2 n ) n≥1 , and positive constants σ1 , σ2 such that for all n ∈ N,E ∥H θn (X n+1 )∥ 2 ≤ σ 2 n + σ1 E ∥∇V (θ n )∥ 2 + σ2 E V (θ n ) -V (θ * ) .In this assumption, r n+1 represents an additive bias term, generally of the order of the square of the bias, and λ n+1 may depend on the minimum eigenvalue of A n . In [21, Theorem 2] , it has been demonstrated that this assumption is weaker than the alternatives used in the literature on biased SGD. We have extended these assumptions to the adaptive case. It is important to note that the first point of H3 depends on the application (objective function V ) and on the adaptive algorithm (matrix A n ) that we want to use. The purpose of this assumption is to provide a more general framework that covers all possible applications and adaptive algorithms. In the biased SGD setting, if the bias term ∥E[H θn (X n+1 ) | F n ] -∇V (θ n )∥ is bounded by bn+1 , we can easily verify the first point of H3 by considering λ n+1 = 1/2 and r n+1 = b2 n+1 . We show in Section 4.3 that this assumption is also verified in algorithms such as Adagrad and RMSProp. The second point of H3 is a weaker assumption compared to bounding the variance of the noise term. Applications where we can verify these assumptions are discussed in Appendix D.We finally consider an additional assumption on A n . Let ∥A∥ be the spectral norm of a matrix A.H4 There exists (β n ) n≥1 such that for all n ∈ N, ∥A n ∥ := λ max (A n ) ≤ β n+1 .In our setting, since A n is assumed to be a symmetric matrix, the spectral norm is equal to the largest eigenvalue. H4 plays a crucial role, as the estimates may diverge when this assumption is not satisfied. Given a sequence (β n ) n≥1 , one way to ensure that H4 is satisfied is to replace the random matricesEQUATIONIt is then clear that ∥ Ãn ∥ ≤ β n+1 . Furthermore, in most cases, especially for Adagrad, RMSProp and Stochastic Newton, control of λ max (A n ) in H4 is satisfied. For example, in Adagrad and RMSProp, in (4), we haveλ max (A n ) ≤ δ -1/2 .4 Convergence ResultsQuestion: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Justification: The paper provides a comprehensive set of assumptions with discussions, along with complete proofs for each theoretical result in the Appendix.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 crossreferenced. • 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. 
convergence under the pl condition: 
 In this section, we study the convergence rate of SGD with biased gradients and adaptive steps under the PL condition. We give below a simplified version of the bound we obtain on the risk and refer to Theorem A.2 in the appendix for a formal statement with explicit constants. Theorem 4.1. Assume that H1 -H4 hold. Let θ n ∈ R d be the n-th iterate of the recursion (3) andEQUATIONThe rate obtained is classical and shows the tradeoff between a term coming from the adaptive steps (with a dependence on γ, β, λ) and a term r n which depends on the control of the bias. To minimize the right hand-side of (7), we would like to have β = λ = 0. For example, it is verified in the case of Adagrad and RMSProp if the gradients are bounded, as will be discussed later.We stress that Theorem 4.1 applies to any adaptive algorithm of the form (3), with the only assumption being H4. Without any information on these eigenvalues, the choice that β n ∝ n β and λ n ∝ n -λ allows us to remain very general, which can even be seen as a worst-case scenario. Finally, note that non-adaptive SGD is a particular case of Theorem 4.1. Thus, our theorem gives new results also in the non-adaptive case with generic step sizes and biased gradients with decreasing bias.In the non-convex smooth case, theoretical results are generally based on a randomized version of SA, as described in [60, 32, 44] . Instead of considering the final parameter θ n , we introduce a random variable R, which takes its values in {1, . . . , n}, and the quantity of interest becomes θ R . Note that this procedure is a technical tool, in practical applications we use classical SA. The following theorem provides a bound in expectation on the gradient of the objective function V , which is the best we can have given that no assumption is made about the existence of a global minimum of V . Theorem 4.2. Assume that H2 -H4 hold. Assume also that for all k ∈ N, we have γ k+1 ≤ λ k+1 /(σ 1 Lβ 2 k+1 ). For any n ≥ 1, let R ∈ {0, . . . , n} be a discrete random variable such that:P(R = k) := w k+1 γ k+1 λ k+1 n j=0 w j+1 γ j+1 λ j+1,where w k+1 = k+1 j=1 (1 + σ2 δ j ) -1 with δ j = Lγ 2 j β 2 j /2. Then, E ∥∇V (θ R )∥ 2 ≤ 2 V * + α 1,n + α 2,n n j=0 w j+1 γ j+1 λ j+1, whereα 1,n = n k=0 w k+1 γ k+1 λ k+1 r k+1 , α 2,n = n k=0 w k+1 δ k+1 σ 2 k , and V * = E[V (θ 0 ) -V (θ * )] .If σ2 = 0, Theorem 4.2 recovers the asymptotic convergence rate obtained by [44] with respect to the hyperparameters γ, β, and λ, and to the bias. We can observe that if γ ≤ λ + 2β, the condition on (γ k ) k≥1 can be met simply by tuning C γ . In particular, if A n = I d , the requirement on the step sizes can be expressed as γ k+1 ≤ 1/(σ 1 L).We give below the convergence rates obtained from Theorem 4.2 under the same assumptions on γ n , β n , and λ n as in Theorem 4.1. Corollary 4.3. Assume that H2-H4 hold. Letγ n = C γ n -γ , β n = C β n β , λ n = C λ n -λ with C γ > 0, C β > 0, and C λ > 0.Assume that γ, β, λ ≥ 0 and γ + λ < 1. Then, if σ2 = 0, we have:E ∥∇V (θ R )∥ 2 =    O n -γ+λ+2β + b n if γ -β < 1/2 , O n γ+λ-1 + b n if γ -β > 1/2 , O n γ+λ-1 log n + b n if γ -β = 1/2 ,where the bias term b n can be constant or decreasing. In the latter case, writing r n = C r n -r , we have:b n =    O (n -r ) if r + λ + γ < 1 , O n γ+λ-1 if r + λ + γ > 1 , O n γ+λ-1 log n if r + λ + γ = 1 .In practice, the value of r is known in advance while the other parameters can be tuned to achieve the optimal rate of convergence. In any scenario, we can never achieve a bound of O(1/ √ n + b n ), and the best rate we can reach is O(log n/ √ n + b n ) when γ = 1/2, β = 0, and λ = 0. In this case, all eigenvalues of A n must be bounded from both below and above. Note that we could also have obtained such a rate by taking λ n = n -1/2 and β n = n -1/2 while keeping γ n constant. However, the assumption that β n = n -1/2 is too strong (fast decay of the eigenvalues of A n ), hence our choice of β n = C β n β . Finally, for a decreasing bias, if r ≥ 1/2, the bias term contributes to the convergence rate of the algorithm. Otherwise, the other term is the leading term of the upper bound. In both cases, the best achievable bound isO(log n/ √ n) if r ≥ 1/2.Bounded Gradient Case. Now, we analyze the convergence of Randomized Adaptive Stochastic Approximation when the stochastic updates are bounded, as given by the following assumption.H5 There exists M ≥ 0 such that for all n ∈ N, ∥H θn (X n+1 ) ∥ ≤ M .Boundedness of the stochastic gradient of the objective function is a classical assumption in adaptive stochastic optimization [67, 74, 19, 73] . Corollary 4.4. Assume that H2-H5 hold. Letγ n = C γ n -γ , β n = C β n β , λ n = C λ n -λ with C γ > 0, C β > 0, and C λ > 0.Assume that γ, β, λ ≥ 0 and γ + λ < 1. For any n ≥ 1, let R ∈ {0, . . . , n} be a uniformly distributed random variable. Then,E ∥∇V (θ R )∥ 2 ≤ V * + α ′ 1,n + LM 2 α ′ 2,n /2 √ n , where α ′ 1,n = n k=0 γ k+1 λ k+1 r k+1 , α ′ 2,n = n k=0 γ 2 k+1 β 2 k+1, andV * = E[V (θ 0 ) -V (θ * )].Importantly, in Corollary 4.4, there are no assumptions on the step sizes, and we obtain a better bound than in Theorem 4.2.We give a convergence analysis of Adagrad and RMSProp with a biased gradient estimator. First, note that, under H5, for all eigenvalues λ of A n , the adaptive matrix in Adagrad or RMSProp, it holds that (M 2 + δ) -1/2 ≤ λ ≤ δ -1/2 , i.e., H4 is satisfied with λ = 0 and β = 0.Corollary 4.5. Assume that H2 and H5 hold. Let γ n = c γ n -1/2 and A n denote the adaptive matrix in Adagrad or RMSProp. For any n ≥ 1, let R ∈ {0, . . . , n} be a uniformly distributed random variable. Suppose that for any n ≥ 1, there exist positive constants α and C α such that:EQUATIONThen,E ∥∇V (θ R )∥ 2 = O log n √ n + b n ,where the bias b n is explicitly given in Appendix A.5.In the case of an unbiased gradient, we obtain the same bound of O(log n/ √ n) as in [74, 79, 19 ] under the same assumptions. If the bias is of the order O(n -1/4 ), the algorithm achieves the same convergence rate as in the case of an unbiased gradient.Algorithm 1 AMSGRAD with Biased Gradients Input: Initial point θ 0 , maximum number of iterations n, step sizes {γ k } k≥1 , momentum parameters ρ 1 , ρ 2 ∈ [0, 1) and regularization parameter δ ≥ 0. Set m 0 = 0, V 0 = 0 and V0 = 0for k = 0 to n -1 do Compute the stochastic update H θ k (X k+1 ) m k = ρ 1 m k-1 + (1 -ρ 1 )H θ k (X k+1 ) V k = ρ 2 V k-1 +(1-ρ 2 )H θ k (X k+1 )H θ k (X k+1 ) ⊤ Vk = max Vk-1 , Diag(V k ) A k = δI d + Vk -1/2 θ k+1 = θ k -γ k+1 A k m k end for Output: (θ k ) 0≤k≤nFinally, we show the convergence of AMS-GRAD [67] with a biased gradient estimator. At each iteration, AMSGRAD uses an exponential moving average of past gradients instead of the current gradient as in Equation (3), which is detailed in Algorithm 1. The key difference between Adam and AMSGRAD lies in their handling of the second moment estimate. Specifically, AMSGRAD uses the updated term Vk = max( Vk-1 , Diag(V k )) instead of directly using V k , with the maximum taken coordinate-wise. This approach is crucial, as it ensures that the eigenvalues of A n decrease at each iteration. The following theorem provides a bound in expectation on the gradient of the objective function V using randomized iterations with AMSGRAD. Theorem 4.6. Assume that H2, H3 (i), and H5 hold. Let γ n = c γ n -1/2 , A n denote the adaptive matrix of AMSGRAD in Algorithm 1, and ρ 1 , ρ 2 ∈ [0, 1). For any n ≥ 1, let R ∈ {0, . . . , n} be a uniformly distributed random variable. Then,E ∥∇V (θ R )∥ 2 = O log n √ n + b n ,where b n corresponds to the bias which comes from r n in H3(i). Choosing r n = C r n -r , we get:b n =    O (n -r ) if r < 1/2 , O n -1/2 if r > 1/2 , O n -1/2 log n if r = 1/2 .If the bias is of the order O(n -1/4 ), we achieve a convergence rate of O(log n/ √ n), which is the same as that of an unbiased gradient [19] and similar to that of Adagrad and RMSProp. It is worth noting that our results are also applicable to SGD momentum by taking A n = I d in Algorithm 1.For illustrative purposes, in this subsection we give the form of the bias of the gradient estimator, denoted by bn , in two simple scenarios, i.e., when {X n , n ∈ N} is either an i.i.d. sequence or a Markov chain. For Adagrad, RMSProp, and AMSGRAD, bounding the bias of the gradient estimator is a sufficient condition for verifying H3(i), which in turn enables us to derive convergence results in each scenario.I.i.d. case. Assume that {X n , n ∈ N} are i.i.d. random variables. If the mean field function h (θ n ) = E [H θn (X n+1 ) | F n ]aligns with the true gradient, then the estimator is unbiased. Otherwise, the bias of the gradient estimator is bn+1 = ∥h(θ n ) -∇V (θ n )∥ .Markov Chain case. Assume now that {X n , n ∈ N} is a Markov Chain. The bias consists of two parts: the difference between the mean field function and the true gradient, and a term due to the Markov chain dynamics. For all T ≥ 0, we define the stochastic update as follows:H θ k (X k+1 ) = 1 T T i=1 H θ k X (i) k+1 ,whereX (i)k+1 represents the i-th sample generated at iteration k + 1. This multi-sample estimator is commonly used in applications such as Reinforcement Learning, Markov Chain Monte Carlo, and Sequential Monte Carlo methods, effectively reducing the variance of the gradient estimator. The mixing time τ mix of a Markov chain with stationary distribution π and transition kernel P is characterized as:τ mix := inf t ; sup x D TV (P t (x, •), π) ≤ 1 4 ,where D TV denotes the total variation distance. For an ergodic Markov chain with stationary distribution π, the bias of this gradient estimator when using T samples per step isbn+1 = ∥h(θ n ) -∇V (θ n )∥ + M τ mix /T ,where h(θ) = H θ (x)π(dx). If the general optimization problem reduces to the following stochastic optimization problem with Markov noise, as considered in most of the literature [28, 24, 7] :min θ∈R d V (θ) := E x∼π [f (θ; x)],where θ → f (θ; x) is a loss function, and π is some stationary data distribution of the Markov Chain and H θ k (XEQUATION), then bn+1 = M τ mix /T , similar to SGD with Markov Noise [24] .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: Our work is primarily theoretical, and the experiments are conducted to illustrate our results. The paper provides detailed descriptions of the experimental setup, parameters, and methodologies, ensuring that the main results can be reproduced and verified.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. 
bilevel and conditional stochastic optimization: 
 We can now apply our theoretical results in various settings where biased gradients are involved. In particular, they apply to the fields of Stochastic Bilevel Optimization and Conditional Stochastic Optimization. Stochastic Bilevel Optimization consists of minimizing an objective function V with respect to θ, where V is itself a function of ϕ * (θ) and ϕ * (θ) is obtained by solving another minimization problem. Conditional Stochastic Optimization focuses on optimizing the expected value of a function that contains a nested conditional expectation on a random variable η. We provide in Table 1 a summary of the assumptions satisfied in these settings, which allow to apply the results of Section 4 and to obtain a O(log n/ √ n + b n ) convergence rate in both cases, and explicit forms for b n . To our knowledge, these are the first convergence rates obtained in these settings.We refer to Appendix D for other examples in which the bias of the estimator can be controlled, in particular Self-Normalized Importance Sampling (Appendix D.1), Sequential Monte Carlo Methods (Appendix D.2), Policy Gradient (Appendix D.3), Zeroth-Order Gradient (Appendix D.4), and Coordinate Sampling (Appendix D.5).In this section, we illustrate our theoretical results in the context of deep VAE. The experiments were conducted using PyTorch [65] , and the source code can be found here 2 . In generative models, θ∈R d E ξ [f (θ, ϕ * (θ); ξ)] min θ∈R d E ξ f ξ E η|ξ [g η (θ, ξ)] s.t. ϕ * (θ) ∈ argmin ϕ∈R q E ζ [g(θ, ϕ; ζ)] Lipchitz Constant H2 Lemma C.2 Lemma C.6 Bias Control H3 Lemma C.3 Lemma C.5 Gradient Bound H5 Lemma C.3 Lemma C.6 Convergence Theorem C.4 Theorem C.7the objective is to maximize the marginal likelihood log p θ (x), which is the marginalization of (x, z) → p θ (x, z), where x represents the observations and z is the latent variable. Under some simple technical assumptions, by Fisher's identity, we have:EQUATIONHowever, in most cases, the conditional density z → p θ (z | x) is intractable and can only be sampled. Variational Autoencoders introduce an additional parameter ϕ and a family of variational distributions z → q ϕ (z | x) to approximate the true posterior distribution. Parameters are estimated by maximizing the Evidence Lower Bound (ELBO):log p θ (x) ≥ E q ϕ (•|x) log p θ (x, Z) q ϕ (Z | x) =: L ELBO (θ, ϕ; x) .The Importance Weighted Autoencoder (IWAE) [12] is a variant of the VAE that incorporates importance weighting to obtain a tighter ELBO. The IWAE objective can be written as follows:L IWAE k (θ, ϕ; x) = E q ⊗k ϕ (•|x) log 1 k k ℓ=1 p θ (x, Z (ℓ) ) q ϕ (Z (ℓ) | x) ,where k corresponds to the number of samples drawn from the variational posterior distribution. The estimator of the gradient of ELBO in IWAE is a biased estimator of ∇ θ log p θ (x). In Theorem B.1, we establish that the bias of this estimator is of order O(1/k), thereby allowing us to derive a convergence rate for IWAE. Since bias has an impact on convergence rates, we propose to use one of the bias reduction techniques, the Biased Reduced Importance Weighted Autoencoder (BR-IWAE) [14] , which is detailed in Appendix B.Dataset and Model. We conduct our experiments on the CIFAR-10 dataset [51] and use a Convolutional Neural Network (CNN) architecture with the Rectified Linear Unit (ReLU) activation function for both the encoder and the decoder. The latent space dimension is set to 100. We estimate the log-likelihood using VAE, IWAE, and BR-IWAE models, all of which are trained for 100 epochs. Then, we illustrate empirically the convergence rates obtained in Corollary 4.5 and Theorem 4.6 for IWAE. Since the bias of the estimator of the gradient in IWAE is of the order O(1/k), choosing a bias of order O(n -α ) is equivalent to using n α samples at iteration n to estimate the gradient. We plot in Figure 2 the gradient squared norm ∥∇V (θ n )∥ 2 and the Negative Log-Likelihood is given in Appendix E.2. Note that all figures are with respect to epochs, whereas here, n represents the number of updates of the gradient. The dashed curves correspond to the expected convergence rate O(n -1/4 ) for α = 1/8 and O(log n/ √ n) for α = 1/4 and α = 1/2. We observe that the algorithms converge at the expected theoretical rates, and even faster. In Appendix E.2, we have included an additional experiment on the FashionMNIST dataset [76] , which shows similar behavior, but the convergence is closer to the expected rates, suggesting that our upper bounds may be tight. We see similar convergence rates for Adagrad, RMSProp, and Adam, although, as expected, Adam performs slightly better. Moreover, it is clear that convergence is faster with a larger α but beyond a certain threshold for α the rate of convergence does not change significantly. Since choosing a larger α induces an additional computational cost, it is crucial to choose an appropriate value that achieves fast convergence without being too computationally expensive. Choosing an optimal number of samples at each iteration remains an open problem depending on the chosen generative model.V( n) 2 = 1/8 (RMSProp) = 1/4 (RMSProp) = 1/2 (RMSProp) = 1/8 (Adam) = 1/4 (Adam) = 1/2 (Adam)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: [Yes] Justification: The code supporting our experiments is available on GitHub, as referenced in the Experiments section. All data used in our study is publicly accessible (CIFAR dataset). Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
discussion: 
 This paper provides a non-asymptotic analysis of Biased Adaptive Stochastic Approximation with and without the PL condition in the non-convex smooth setting. We derive a convergence rate of O(log n/ √ n + b n ) for non-convex smooth functions, where b n corresponds to the time-dependent decreasing bias, and an improved linear convergence rate with the Polyak-Łojasiewicz (PL) condition. We also establish that Adagrad, RMSProp, and AMSGRAD with biased gradients converge to critical points for non-convex smooth functions. Our results provide insights on hyper-parameters tuning to achieve fast convergence and reduce computational time. A natural extension of this work is the analysis of the assumptions, the bias and convergence rates for specific deep learning architectures. A theoretical analysis of the Monte Carlo effort required at each iteration to obtain an optimal convergence rate is another interesting perspective.[75] Xidong Wu, Jianhui Sun, Zhengmian Hu, Junyi Li, Aidong Zhang, and Heng Huang. Federated conditional stochastic optimization. In Advances in Neural Information Processing Systems, volume 36, 2024.[76] Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv:1708.07747, 2017.[77] Manzil Zaheer, Sashank Reddi, Devendra Sachan, Satyen Kale, and Sanjiv Kumar. Adaptive methods for nonconvex optimization. In Advances in neural information processing systems, volume 31, 2018.[78] Matthew D Zeiler. Adadelta: an adaptive learning rate method. arXiv preprint arXiv:1212.5701, 2012.[79] Fangyu Zou, Li Shen, Zequn Jie, Ju Sun, and Wei Liu. Weighted adagrad with unified momentum. arXiv preprint arXiv:1808.03408, 2018.Supplementary Material for "Non-asymptotic Analysis of Biased Adaptive Stochastic Approximation" We first establish a technical lemma which is essential for the proof.Lemma A.1. Let (δ n ) n≥0 , (γ n ) n≥1 , (η n ) n≥1, and (v n ) n≥1 be some positive sequences satisfying the following assumptions.• The sequence δ n follows the recursive relation:δ n ≤ (1 -2ωγ n + η n γ n ) δ n-1 + v n γ n ,with δ 0 ≥ 0 and ω > 0.• Let n 0 = inf {n ≥ 1 : η n ≤ ω}, then for all n ≥ n 0 + 1, we assume that ωγ n ≤ 1.Then, for all n ∈ N,δ n ≤ exp -ω n k=n/2 γ k exp 2 n k=1 η k γ k δ 0 + 2 max 1≤k≤n v k η k + 1 ω max n/2≤k≤n v k .The proof is given in [35, Proposition 6.1] Theorem A.2. Assume that H1 -H4 hold. Let θ n ∈ R d be the n-th iterate of the recursion (3). Then,E [V (θ n ) -V (θ * )] ≤ E [V (θ 0 ) -V (θ * )] + 2 σ max 1≤k≤n λ k+1 v k β 2 k+1 γ k+1 exp - µ 2 n k=n/2 λ k+1 γ k+1 × exp 2 n k=1 C k β 2 k+1 γ 2 k+1 + 2 µ max n/2≤k≤n v k ,whereσ = σ2 L 2 + σ1 L 2 , C k = max σ, µ 2 λ 2 k+1 4β 2 k+1 and v k = r k+1 + Lσ 2 k 2 β 2 k+1 λ k+1 γ k+1 .with the conventionC k = 1 if σ1 = σ2 = 0.Proof. Since V is L-smooth (Assumption H2) and using the recursion (3) of Adaptive SA, we obtain:V (θ n+1 ) ≤ V (θ n ) + ⟨∇V (θ n ) , θ n+1 -θ n ⟩ + L 2 ∥θ n+1 -θ n ∥ 2 ≤ V (θ n ) -γ n+1 ⟨∇V (θ n ) , A n H θn (X n+1 )⟩ + Lγ 2 n+1 2 ∥A n ∥ 2 ∥H θn (X n+1 )∥ 2 . Writing V n = V (θ n ) -V (θ * ), we get V n+1 ≤ V n -γ n+1 ⟨∇V (θ n ) , A n H θn (X n+1 )⟩ + L 2 γ 2 n+1 β 2 n+1 ∥H θn (X n+1 )∥ 2 .Then, using H3,E [V n+1 ] ≤ E [V n ] -γ n+1 E ⟨∇V (θ n ) , A n H θn (X n+1 )⟩ + L 2 γ 2 n+1 β 2 n+1 σ 2 n + L 2 γ 2 n+1 β 2 n+1 σ1 E[∥∇V (θ n )∥ 2 ] + σ2 E[V n ] ≤ 1 + σ2 L 2 β 2 n+1 γ 2 n+1 E[V n ] -γ n+1 λ n+1 - σ1 L 2 γ n+1 β 2 n+1 E[∥∇V (θ n )∥ 2 ] + γ n+1 λ n+1 r n+1 + Lσ 2 n 2 γ 2 n+1 β 2 n+1 .Given the smoothness condition (5) , with θ ′ = θ n and θ = θ n -1 L ∇V (θ n ), we derive:V (θ * ) ≤ V (θ n ) - 1 L ∥∇V (θ n )∥ 2 + 1 2L ∥∇V (θ n )∥ 2 , ∥∇V (θ n )∥ 2 ≤ 2LV n .Using the above inequality and the Polyak-Łojasiewicz condition (H1), we obtain:E [V n+1 ] ≤ 1 -µλ n+1 γ n+1 + σ2 L 2 + σ1 L 2 β 2 n+1 γ 2 n+1 E [V n ] + γ n+1 λ n+1 r n+1 + Lσ 2 n 2 γ 2 n+1 β 2 n+1 .By choosing γn+1 = λ n+1 γ n+1 , we get:E [V n+1 ] ≤ 1 -µγ n+1 + σ2 L 2 + σ1 L 2 β 2 n+1 λ 2 n+1 γ2 n+1 E [V n ] + γn+1 r n+1 + Lσ 2 n 2 β 2 n+1 λ n+1 γn+1 γ n+1 .In order to satisfy the assumptions of Lemma A.1, considerC n = max σ, µ 2 λ 2 n+1 4β 2 n+1 with σ = σ2 L 2 + σ1 L 2 .Then, since C n ≥ σ, we have:E [V n+1 ] ≤ 1 -µγ n+1 + C n β 2 n+1 λ 2 n+1 γ2 n+1 E [V n ] + γn+1 r n+1 + Lσ 2 n 2 β 2 n+1 λ n+1 γn+1 γ n+1 .Now, using lemma A.1 by choosing:δ n = E [V n ] , η n = C n β 2 n+1 λ 2 n+1 γn+1 , ω = µ 2 , v n = r n+1 + Lσ 2 n 2 β 2 n+1 λ n+1 γ n+1 ,we have:E [V (θ n ) -V (θ * )] ≤ E [V (θ 0 ) -V (θ * )] + 2 σ max 1≤k≤n v k λ 2 k+1 β 2 k+1 γk+1 exp - µ 2 n k=n/2 γk+1 × exp 2 n k=1 C k β 2 k+1 γ2 k+1 /λ 2 k+1 + 2 µ max n/2≤k≤n {v k } ,which concludes the proof by taking γn+1 = λ n+1 γ n+1 .A.2 Proof of Theorem 4.2 By H2, using (5), we obtain:V (θ k+1 ) ≤ V (θ k ) + ⟨∇V (θ k ) , θ k+1 -θ k ⟩ + L 2 ∥θ k+1 -θ k ∥ 2 ,which, using the recursion (3) of Adaptive SA and H4, yields:V (θ k+1 ) ≤ V (θ k ) -γ k+1 ⟨∇V (θ k ) , A k H θ k (X k+1 )⟩ + δ k+1 ∥H θ k (X k+1 )∥ 2 ,with δ k+1 = Lγ 2 k+1 β 2 k+1 /2. Using Assumption H3, we have:E[V (θ k+1 )] ≤ E[V (θ k )] -γ k+1 λ k+1 E ∥∇V (θ k )∥ 2 + γ k+1 λ k+1 r k+1 + δ k+1 σ 2 k + δ k+1 σ1 E[∥∇V (θ k )∥ 2 ] + σ2 E [V (θ k ) -V (θ * )] .Therefore,γ k+1 λ k+1 - Lσ 1 2 γ k+1 β 2 k+1 E ∥∇V (θ k )∥ 2 ≤ (1 + σ2 δ k+1 ) E [V (θ k )] -V (θ * ) -E [V (θ k+1 )] -V (θ * ) + γ k+1 λ k+1 r k+1 + δ k+1 σ 2 k .Let us now consider the sequence of weights w k defined by w 0 = 1 andw k = k j=1 (1 + σ2 δ j ) -1 .Then,w k+1 γ k+1 λ k+1 - Lσ 1 2 γ k+1 β 2 k+1 E ∥∇V (θ k )∥ 2 ≤ w k E[V (θ k )] -V (θ * ) -w k+1 E [V (θ k+1 )] -V (θ * ) + w k+1 γ k+1 λ k+1 r k+1 + w k+1 δ k+1 σ 2 k .In the sequel, let us denoteV n = V (θ n ) -V (θ * ), so that n k=0 w k+1 γ k+1 λ k+1 1 - Lσ 1 2λ k+1 γ k+1 β 2 k+1 E ∥∇V (θ k )∥ 2 ≤ w 0 E [V 0 ] -w n+1 E [V n+1 ] + 1 2 n k=0 w k+1 γ k+1 λ k+1 r k+1 + n k=0 w k+1 δ k+1 σ 2 k .Then, given that γ k+1 ≤ λ k+1 /(Lσ 1 β 2 k+1 ), we have1 2 E n k=0 w k+1 γ k+1 λ k+1 ∥∇V (θ k )∥ 2 ≤ w 0 E [V 0 ] -w n+1 E [V n+1 ] + 1 2 n k=0 w k+1 γ k+1 λ k+1 r k+1 + n k=0 w k+1 δ k+1 σ 2 k .Consequently, by definition of the discrete random variable R,E ∥∇V (θ R )∥ 2 = n k=0 w k+1 γ k+1 λ k+1 n j=0 w j+1 γ j+1 λ j+1 E ∥∇V (θ k )∥ 2 ≤ 2 E [V 0 ] -w n+1 E [V n+1 ] + n k=0 w k+1 γ k+1 r k+1 + n k=0 w k+1 δ k+1 σ 2 k n j=0 w j+1 γ j+1 λ j+1, which concludes the proof by noting that V (θ n+1 ) ≥ V (θ * ).The proof is a direct consequence of the fact that for a sufficiently large n:n k=1 1 k s =    O n -s+1 if 0 ≤ s < 1 , O (1) if s > 1 , O (log n) if s = 1 .A.4 Proof of Corollary 4.4By H2, using (5), we obtain:V (θ k+1 ) ≤ V (θ k ) + ⟨∇V (θ k ) , θ k+1 -θ k ⟩ + L 2 ∥θ k+1 -θ k ∥ 2 ≤ V (θ k ) -γ k+1 ⟨∇V (θ k ) , A k H θ k (X k+1 )⟩ + Lγ 2 k+1 2 ∥A k ∥ 2 ∥H θ k (X k+1 )∥ 2 ,which, using H4 and H5 yields:V (θ k+1 ) ≤ V (θ k ) -γ k+1 ⟨∇V (θ k ) , A k H θ k (X k+1 )⟩ + L 2 γ 2 k+1 β 2 k+1 M 2 .Using H3,E[V (θ k+1 )|F k ] ≤ V (θ k ) -γ k+1 λ n+1 ∥∇V (θ k )∥ 2 + γ k+1 λ k+1 r k+1 + LM 2 2 γ 2 k+1 β 2 k+1 .Therefore,γ k+1 λ k+1 ∥∇V (θ k )∥ 2 ≤ V (θ k ) -E [V (θ k+1 ) |F k ] + γ k+1 λ k+1 r k+1 + LM 2 2 γ 2 k+1 β 2 k+1 ,andn k=0 γ k+1 λ k+1 E ∥∇V (θ k )∥ 2 ≤ E [V (θ 0 ) -V (θ n+1 )] + n k=0 γ k+1 λ k+1 r k+1 + LM 2 2 n k=0 γ 2 k+1 β 2 k+1 .Consequently, by definition of the discrete random variable R,E ∥∇V (θ R )∥ 2 = 1 n n k=0 E ∥∇V (θ k )∥ 2 ≤ n k=0 γ k+1 λ k+1 √ n E ∥∇V (θ k )∥ 2 ≤ V 0,n + n k=0 γ k+1 λ k+1 r k+1 + LM 2 n k=0 γ 2 k+1 β 2 k+1 /2 √ n , where V 0,n = E[V (θ 0 ) -V (θ n+1 )],which conclude the proof by noting that V (θ n+1 ) ≥ V (θ * ).A.5 Proof of Corollary 4.5Here, we consider the case where the regularization is non-increasing, i.e., where δ = β -2 n+1 . The constant case is strictly analogous. To verify H4, we demonstrate that the control of the maximum and minimum eigenvalues is satisfied for Adagrad and RMSProp.• Lower bound for the smallest eigenvalue of A n . By assumption H5, we have:1 n + 1 n k=0 H θ k (X k+1 )H θ k (X k+1 ) ⊤ ≤ M 2 .This implies that:λ min (A n ) = λ max β -2 n+1 I d + Diag 1 n + 1 n k=0 H θ k (X k+1 )H θ k (X k+1 ) ⊤ -1/2 ≥ (β -2 1 + M 2 ) -1/2 .• Upper bound for the largest eigenvalue of A n .λ max (A n ) = λ min β -2 n+1 I d + Diag 1 n + 1 n k=0 H θ k (X k+1 )H θ k (X k+1 ) ⊤ -1/2 ≤ β n+1 .Therefore, by setting λ n+1 = (β -2 1 + M 2 ) -1/2 and β n = C β n β , we have λ = 0 and one can arbitrarily choose β (one can take β = 0 for the constant regularization case).• Lower bound for the smallest eigenvalue of A n . By assumption H5, we have:∥V n ∥ ≤ (1 -ρ) n k=1 ρ n-k ∥H θ k (X k+1 )∥ 2 ≤ M 2 (1 -ρ) n k=1 ρ n-k ≤ M 2 ,where we used the fact thatn k=1 ρ n-k ≤ (1 -ρ) -1. This implies that:λ min (A n ) = λ max β -2 n+1 I d + Diag (V n ) -1/2 ≥ (β -2 1 + M 2 ) -1/2 . • Upper bound for the largest eigenvalue of A n . Note that λ max (A n ) = λ min β -2 n+1 I d + Diag (V n ) -1/2 ≤ β n+1 .Therefore, under H2, H3(i), and H5, we can conclude thatE ∥∇V (θ R )∥ 2 = O log n √ n + b n ,where b n corresponds to the bias which comes from r n in H3(i). Choosing r n = C r n -r , we get:b n =    O (n -r ) if r < 1/2 , O n -1/2 if r > 1/2 , O n -1/2 log n if r = 1/2 .Now, we show that under the control of bias, i.e.,∥E [H θn (X n+1 ) |F n ] -∇V (θ n )∥ ≤ C α n -α ,we can verify H3(i) with a similar bound on the bias, where r = 2α. This yields the bias term b n as follows:b n =    O n -2α if α < 1/4 , O n -1/2 if α > 1/4 , O n -1/2 log n if α = 1/4 .Verifying Assumption H3 (i) for Adagrad. Using the tower property, we have:E [⟨∇V (θ n ) , A n H θn (X n+1 )⟩] = E [E [⟨∇V (θ n ) , A n H θn (X n+1 )⟩ |F n ]] ,where (F n ) n≥0 represents the filtration generated by the random variables (θ 0 , {X k } k≤n ). Let Ãn be an adaptive F n -measurable matrix. Then,E [⟨∇V (θ n ) , A n H θn (X n+1 )⟩ |F n ] = ∇V (θ n ) , Ãn E [H θn (X n+1 ) |F n ]Treated as in SGD but with λmin( Ãn)+ E ∇V (θ n ) , (A n -Ãn )H θn (X n+1 ) |F nControl error between An and Ãn .We only verify Assumption H3(i) for Adagrad algorithm since it is analogous to RMSProp. Consider A n given by:A n = diag β -2 n+1 I d + 1 n + 1 n k=0 H θ k (X k+1 ) H θ k (X k+1 ) ⊤ -1/2 .Ãn = diag β -2 n+1 I d + 1 n + 1 n-1 k=0 H θ k (X k+1 ) H θ k (X k+1 ) ⊤ -1/2and denoting by A[i] the i-th element of the diagonal of a matrix A, we haveA n [i] -Ãn [i] = u -1/2 n v 1/2 n -u 1/2 n v -1/2 n ≤ 0 , where u n = β -2 n+1 + 1 n + 1 n k=0 (H θ k (X k+1 ) [i]) 2 and v n = β -2 n+1 + 1 n + 1 n-1 k=0 (H θ k (X k+1 ) [i]) 2 .Then, sinceu n ≥ v n , A n [i] -Ãn [i] = v n -u n √ u n v n √ u n + √ v n ≥ - 1 n + 1 (H θn (X n+1 ) [i]) 2 1 2v n 3/2 ≥ - β 3 n+1 n + 1 (H θn (X n+1 ) [i]) 2 .Since the bias ofH θn (X n+1 ) is bounded by bn := C α n -α , E [⟨∇V (θ n ) , A n H θn (X n+1 )⟩ |F n ] = ∇V (θ n ) , Ãn E [H θn (X n+1 ) |F n ] + E ∇V (θ n ) , (A n -Ãn )H θn (X n+1 ) |F n ≥ λ min Ãn ∥∇V (θ n )∥ 2 -λ max Ãn ∥∇V (θ n )∥ bn -∥∇V (θ n )∥ β 3 n+1 n + 1 E ∥H θn (X n+1 )∥ 3 |F n .As H θn (X n+1 ) and the gradient of V are uniformly bounded by M , λ min ( Ãn )≥ (β -2 1 + M 2 ) -1/2 , so that E [⟨∇V (θ n ) , A n H θn (X n+1 )⟩ |F n ] ≥ 1 β -2 1 + M 2 ∥∇V (θ n )∥ 2 -β n+1 M bn -M 4 β 3 n+1 n + 1 ,and Assumption H3(i) is satisfied withλ n+1 = (β -2 1 + M 2 ) -1/2 and r n+1 = M β 2 n+1 b2 n /λ n+1 + M 4 β 3 n+1 /(n + 1).A.6 Proof of Theorem 4.6The proof of this theorem is inspired by [67] and [73] , considering biased gradient estimators and decreasing step sizes. We define the operation max(D 1 , D 2 ) for diagonal matrices D 1 and D 2 as the matrix formed by taking the maximum between the diagonal elements of D 1 and D 2 . We say that the sequence (A n ) n≥1 of diagonal matrices is decreasing if all diagonal terms are decreasing, in other words, if all eigenvalues are decreasing.Let θk+1 = θ k+1 + κ (θ k+1 -θ k ), for k ≥ 1, κ ∈ [0, 1) and m k = ρ 1 m k-1 + (1 -ρ 1 )g k with g k = H θ k (X k+1 ).Using the recursion of AMSGRAD, we have:θk+1 -θk = (1 + κ)θ k+1 -(1 + 2κ)θ k + κθ k-1 = (1 + κ) (θ k+1 -θ k ) -κ (θ k -θ k-1 ) = -(1 + κ)γ k+1 A k m k + κγ k A k-1 m k-1 . Choosing κ = ρ 1 /(1 -ρ 1 ), we can rewrite it as:θk+1 -θk = κ (γ k A k-1 -γ k+1 A k ) m k-1 -γ k+1 A k g k .By Assumption H2, V is L-smooth, using the recursion of AMSGRAD together with a Taylor expansion with θk , we obtain:V ( θk+1 ) ≤ V θk + ∇V θk , θk+1 -θk + L 2 θk+1 -θk 2 ≤ V θk -γ k+1 ∇V θk , A k g k + κ ∇V θk , (γ k A k-1 -γ k+1 A k ) m k-1 + Lγ 2 k+1 ∥A k g k ∥ 2 + Lκ 2 ∥(γ k A k-1 -γ k+1 A k ) m k-1 ∥ 2 ≤ V θk + T 1,k + T 2,k + T 3,k + T 4,k ,whereT 1,k = -γ k+1 ⟨∇V (θ k ) , A k g k ⟩ + Lγ 2 k+1 ∥A k g k ∥ 2 , T 2,k = -γ k+1 ∇V θk -∇V (θ k ) , A k g k , T 3,k = κ ∇V θk , (γ k A k-1 -γ k+1 A k ) m k-1 , T 4,k = Lκ 2 ∥(γ k A k-1 -γ k+1 A k ) m k-1 ∥ 2 .Note first thatn k=1 E [T 1,k ] = - n k=1 γ k+1 E [⟨∇V (θ k ) , A k g k ⟩] + L n k=1 γ 2 k+1 E ∥A k g k ∥ 2 ≤ -C λ n k=1 γ k+1 E ∥∇V (θ k )∥ 2 + C λ n k=1 γ k+1 r k+1 + L n k=1 γ 2 k+1 E ∥A k g k ∥ 2 ,whereC λ = (δ + M 2 ) -1/2 .For the second term, using the inequality xy ≤ x 2 /2 + y 2 /2 for all x, y, and the smoothness of V , we get:n k=1 E [T 2,k ] = - n k=1 E ∇V θk -∇V (θ k ) , γ k+1 A k g k ≤ 1 2 n k=1 E ∇V θk -∇V (θ k ) 2 + 1 2 n k=1 E ∥γ k+1 A k g k ∥ 2 ≤ L 2 2 n k=1 E θk -θ k 2 + n k=1 γ 2 k+1 2 E ∥A k g k ∥ 2 ≤ κ 2 L 2 2 n k=1 E ∥θ k -θ k-1 ∥ 2 + n k=1 γ 2 k+1 2 E ∥A k g k ∥ 2 ≤ κ 2 L 2 2 n k=1 γ 2 k E ∥A k-1 m k-1 ∥ 2 + n k=1 γ 2 k+1 2 E ∥A k g k ∥ 2 .For the third term, using the boundedness of the gradient of V and the fact that ∥m k ∥ ≤ M by Lemma A.3, we have:n k=1 E [T 3,k ] = κ n k=1 E ∇V θk , (γ k A k-1 -γ k+1 A k ) m k-1 ≤ κM 2 d i=1 n k=1 E [γ k A k-1 [i] -γ k+1 A k [i]] ≤ κM 2 d i=1 E [γ 1 A 0 [i] -γ n+1 A n [i]] ≤ κM 2 dC γ ,where in the second inequality, we used the fact that γ k and A k are decreasing since we use Vk = max( Vk-1 , V k ). For the last term, using the boundedness of the gradient of V yields:n k=1 E [T 4,k ] = Lκ 2 n k=1 E ∥(γ k A k-1 -γ k+1 A k ) m k-1 ∥ 2 ≤ Lκ 2 M 2 d i=1 n k=1 E (γ k A k-1 [i] -γ k+1 A k [i]) 2 ≤ Lκ 2 M 2 d i=1 n k=1 E (γ k A k-1 [i]) 2 -(γ k+1 A k [i]) 2 ≤ Lκ 2 M 2 dC 2 γ ,where we used the inequality (x -y) 2 ≤ x 2 -y 2 when x ≥ y in the second last inequality.Combining all these terms, we finally obtain:C λ n k=1 γ k+1 E ∥∇V (θ k )∥ 2 ≤ V * + C λ n k=1 γ k+1 r k+1 + L n k=1 γ 2 k+1 E ∥A k g k ∥ 2 + n k=1 γ 2 k+1 2 E ∥A k g k ∥ 2 + κ 2 L 2 2 n k=1 γ 2 k E ∥A k-1 m k-1 ∥ 2 + κM 2 dC γ + Lκ 2 M 2 dC 2 γ ,whereV * = E[V (θ 0 ) -V (θ * )] ≥ E[V (θ 0 ) -V ( θn+1 )].Choosing γ n = n -1/2 and using Lemma A.3 and [15, Lemma 24] yieldsn k=1 γ 2 k+1 E ∥A k m k ∥ 2 ≤ (1 -ρ 1 ) n k=1 γ 2 k+1 E ∥A k g k ∥ 2 ≤ (1 -ρ 1 )dC 2 γ log 1 + nM 2 δ = O (d log n) .Therefore, by dividing both sides by C λ n -1/2 , we obtain1 n n k=1 E ∥∇V (θ k )∥ 2 = O 1 √ n + d log n √ n + d √ n + b n ,which concludes the proof.Lemma A.3. Let γ k+1 ≤ γ k for all k ∈ N, and let A k be the adaptive matrix defined in Algorithm 1.Assume that ρ 1 ∈ [0, 1). Then, for all k ∈ N:∥m k ∥ ≤ M and n k=1 γ 2 k+1 E ∥A k m k ∥ 2 ≤ (1 -ρ 1 ) n k=1 γ 2 k+1 ∥A k g k ∥ 2 .Proof. For the first inequality, we have:∥m k ∥ = (1 -ρ 1 ) k ℓ=1 ρ k-ℓ 1 g ℓ ≤ (1 -ρ 1 ) k ℓ=1 ρ k-ℓ 1 ∥g ℓ ∥ ≤ M (1 -ρ 1 ) ℓ≥0 ρ ℓ 1 ≤ M ,where we used the fact that ℓ≥0 ρ ℓ 1 = 1/(1 -ρ 1 ). For the second inequality, using the fact that γ k and A k are decreasing (in the sense that all eigenvalues of A k are decreasing), since we use Vk = max( Vk-1 , V k ), we can write:n k=1 γ 2 k+1 ∥A k m k ∥ 2 = n k=1 γ 2 k+1 A k (1 -ρ 1 ) k ℓ=1 ρ k-ℓ 1 g ℓ 2 ≤ (1 -ρ 1 ) 2 n k=1 γ 2 k+1 k ℓ=1 ρ k-ℓ 1 ∥A ℓ g ℓ ∥ 2 ≤ (1 -ρ 1 ) 2 n k=1 k ℓ=1 ρ k-ℓ 1 γ 2 ℓ+1 ∥A ℓ g ℓ ∥ 2 ≤ (1 -ρ 1 ) 2 n ℓ=1 n k=ℓ ρ k-ℓ 1 γ 2 ℓ+1 ∥A ℓ g ℓ ∥ 2 ,which concludes the proof.In our case, we have a dependence on δ in the logarithm, which is common for adaptive algorithms. The regularization parameter δ, originally introduced to avoid the zero denominator issue when V k approaches 0, is often overlooked. However, it has been empirically observed that the performance of adaptive methods can be sensitive to the choice of this parameter, especially when a very small δ is used, which has resulted in performance issues in some applications.In practice, δ is typically chosen as 10 -8 . In our convergence rate analysis, even though the logarithm of δ -1 is small, it still impacts the convergence rate. A larger δ will lead to a better convergence rate, while a smaller δ will preserve stronger adaptivity. We need to find a better compromise between the convergence rate and the adaptivity to choose δ. In [77, 67, 73] , it was shown that by choosing δ between 10 -3 and 10 -1 , better results were obtained in some applications of deep learning.Furthermore, several modified versions of Adam have been proposed, such as AMSGRAD [77] and YOGI [67] with the discussion of the regularization parameter δ. The authors of [73] proposed a new modified version of Adam called SADAM to represent the calibrated ADAM using the softplus function. In this algorithm, they define Vk = softplus √ V k while other terms remain unchanged. Since we have:In this section, we elaborate on the IWAE procedure within our framework to illustrate its convergence rate. The IWAE objective function is defined as:L IWAE k (θ, ϕ; x) = E q ⊗k ϕ (•|x) log 1 k k ℓ=1 p θ (x, Z (ℓ) ) q ϕ (Z (ℓ) | x) ,where k corresponds to the number of samples drawn from the encoder's approximate posterior distribution. Denoting V as the objective function, i.e., V (θ) = log p θ (x), the gradient of V and the estimator of the gradient of the ELBO of the IWAE objective are given by:EQUATIONwhere w (ℓ) = p θ (x, z (ℓ) )/q ϕ (z (ℓ) |x) the unnormalized importance weights. Theorem B.1 provides an upper bound for the bias of this estimator. Theorem B.1. Let X ⊆ R dx and Z ⊆ R dz denote the data space and the latent space, respectively. Assume that there exists M such that for all θ ∈ Θ ⊂ R d , x ∈ X and z ∈ Z, ∥∇ θ log p θ (x, z)∥ ≤ M (x). Then, there exists a constant C > 0 such that for all θ ∈ Θ, ϕ ∈ Φ and x ∈ X,E q ⊗k ϕ (•|x) ∇ θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) ≤ C k ,where ∇ θ V (θ) and ∇ θ L IWAE k (θ, ϕ; x) are defined in (10) .Proof. The proof is adapted from [1, Theorem 2.1]. By definition,∇ θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) = k ℓ=1 w (ℓ) ∇ θ log p θ (x, z (ℓ) ) -E p θ (•|x) [∇ θ log p θ (x, z)] k ℓ=1 w (ℓ).Writing H(x, z (ℓ) ) = ∇ θ log p θ (x, z (ℓ) ) -E p θ (•|x) [∇ θ log p θ (x, z)], yields ∇θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) = k ℓ=1 w (ℓ) H(x, z (ℓ) ) k ℓ=1 w (ℓ).Since E q ϕ [w H(x, z)] = 0, we have:∇ θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) = 1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) -E q ϕ w H(x, z) 1 k k ℓ=1 w (ℓ).As k ℓ=1 w (ℓ) H(x, z (ℓ) )/k is an unbiased estimator of E q ϕ w H(x, z) , E q ϕ ∇θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) = E q ϕ 1 1 k k ℓ=1 w (ℓ)-1E q ϕ [w] 1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) -E q ϕ w H(x, z) ,so thatE q ϕ ∇θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) = E q ϕ   1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) -E q ϕ w H(x, z) E q ϕ [w] -1 k k ℓ=1 w (ℓ) E q ϕ [w] 1 k k ℓ=1 w (ℓ)   .Therefore,E q ϕ ∇θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) ≤ A 1 + A 2 ,whereA 1 = E q ϕ ∇θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) 1 { 2 k k ℓ=1 w (ℓ) >Eq ϕ [w]} , A 2 = E q ϕ ∇ θ L IWAE k (θ, ϕ; x) -∇ θ V (θ) 1 { 2 k k ℓ=1 w (ℓ) ≤Eq ϕ [w]} .Note thatA 1 ≤ E q ϕ 2 E q ϕ [w] 2 1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) -E q ϕ w H(x, z) E q ϕ [w] - 1 k k ℓ=1 w (ℓ)≤ 2E q ϕ [w] 2 E q ϕ 1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) -E q ϕ w H(x, z) 1 k k ℓ=1 w (ℓ) -E q ϕ [w] ≤ 2 E q ϕ [w] 2 E q ϕ   1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) -E q ϕ w H(x, z) 2   1/2 × E q ϕ   1 k k ℓ=1 w (ℓ) -E q ϕ [w] 2   1/2, where we used Cauchy-Schwarz inequality in the last inequality. On the other hand,EQUATIONE q ϕ   1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) -E q ϕ w H(x, z) 2   = Tr V 1 k k ℓ=1 w (ℓ) H(x, z (ℓ) ) ≤ 4dM 2 E q ϕ w 2 k .Finally, we deduce thatA 1 ≤ 2 E q ϕ [w] 2 1 √ k E q ϕ w 2 1/2 2 √ dM √ k E q ϕ w 2 1/2 = E q ϕ w 2 E q ϕ [w] 2 4 √ dM k .Using the assumption on the boundedness of ∥∇ θ log p θ (x, z)∥ and the Markov inequality, we obtain:A 2 ≤ 2M P 2 1 k k ℓ=1 w (ℓ) ≤ E q ϕ [w] ≤ 2M P 2 1 k k ℓ=1 w (ℓ) -E q ϕ [w] ≤ -E q ϕ [w] ≤ 2M P 1 k k ℓ=1 w (ℓ) -E q ϕ [w] ≥ E q ϕ [w] 2 ≤ E q ϕ w 2 E q ϕ [w] 2 8M k ,which concludes the proof.Input: Initial point θ 0 , maximum number of iterations n, step sizes {γ k } k≥1 and a hyperparameter α ≥ 0 to control the bias and MSE. for k = 0 to n -1 do Compute the stochastic update ∇ θ,ϕ L IWAE k α (θ k , ϕ k ; X k+1 ) using k α samples from the variational posterior distribution and adaptive stepsA k . Set θ k+1 = θ k -γ k+1 A k ∇ θ L IWAE k α (θ k , ϕ k ; X k+1 ). Set ϕ k+1 = ϕ k -γ k+1 A k ∇ ϕ L IW AE k α (θ k , ϕ k ; X k+1 ). end for Output: (θ k ) 0≤k≤nIn this section, we provide additional details on the Biased Reduced Importance Weighted Autoencoder (BR-IWAE). In IWAE, instead of estimating the gradient of the ELBO with respect to θ via the Monte Carlo method, we estimate the gradient of the true objective function E p θ (•|x) [∇ θ log p θ (x, z)] using the BR-SNIS estimator [14] . This estimator aims to reduce the bias of self-normalized importance sampling estimators without increasing the variance.Input: Maximum number of iterations t max for MCMC and number of samples k from the variational distribution q ϕ (• | x). Initialization: Draw z0 from the variational distribution q ϕ (• | x). The BR-SNIS estimator of E p θ (•|x) [∇ θ log p θ (x, z)] is given by:∇ θ log p θ (x, z 1:k t0:tmax ) = 1 t max -t 0 tmax t=t0+1 k ℓ=1 ω (ℓ) t ∇ θ log p θ (x, z ℓ t ) ,where t 0 corresponds to a burn-in period. By [14, Theorem 4] the bias of this estimator decreases exponentially with t 0 . The BR-IWAE algorithm proceeds in two steps, which are repeated during optimization:• Update the parameter ϕ as in the IWAE algorithm, that is, for all n ≥ 1:ϕ n+1 = ϕ n -γ n+1 A n ∇ ϕ L IWAE k (θ n , ϕ n ; X n+1 ) .• Update the parameter θ by estimating (9) using BR-SNIS as detailed in Algorithm 3:ϕ n+1 = ϕ n -γ n+1 A n ∇ θ log p θ (X n+1 , z 1:k t0:tmax ) .In the previous section, we discussed one technique for reducing bias, BR-IWAE. Here, we provide an overview of some other bias reduction techniques within our context. First, the jackknife biascorrected estimator [62] is defined as:L Jackknife (θ, ϕ; x) = kL IWAE k (θ, ϕ; x) -(k -1)L IWAE k-1 (θ, ϕ; x) ,which achieves a reduced bias of O(k -2 ). This can also be generalized to have a bias of order O(k -m ) for some m ≥ 1 by considering:L Jackknife k,m = m j=0 c(k, m, j)L IWAE k-j ,where the coefficients c(k, m, j) are given asc(k, m, j) = (-1) j (k -j) m (m -j)!j! .The Delta method Variational Inference (DVI) [71] is defined by:L DV I k = E q ⊗k ϕ (•|x) log 1 k k ℓ=1 w (ℓ) + s2 k 2k wk ,wherew (ℓ) = p θ (x, z (ℓ) ) q ϕ (z (ℓ) | x) , wk = 1 k k ℓ=1 w (ℓ) and s2 k = 1 k -1 k ℓ=1 (w (ℓ) -wk ) 2 .The Monte Carlo estimator of the Delta method Variational Inference objective achieves a reduced bias of O(k -2 ). Some other techniques for reducing bias include the iterated bootstrap for bias correction, the debiasing lemma [57] , and Multi-Level Monte Carlo and its variants [39] .We consider the Stochastic Bilevel Optimization problem given by: minθ∈R d V (θ) = E ξ [f (θ, ϕ * (θ); ξ)] (upper-level) (11) subject to ϕ * (θ) ∈ argmin ϕ∈R q E ζ [g(θ, ϕ; ζ)] (lower-level)where the upper and inner level functions f and g are both jointly continuously differentiable and ξ and ζ are random variables. The goal of equation ( 11) is to minimize the objective function V with respect to θ, where ϕ * (θ) is obtained by solving the lower-level minimization problem. This bilevel problem involves many machine learning problems with a hierarchical structure, which include hyper-parameter optimization [31] , metalearning [30] , policy optimization [37] and neural network architecture search [55] . The gradient of the objective function V is given by:∇V (θ) = ∇ θ f (θ, ϕ * (θ)) -∇ θϕ g(θ, ϕ * (θ))v *, where v * is the solution of the following linear system:∇ 2 ϕ g(θ, ϕ * (θ))v = ∇ ϕ f (θ, ϕ * (θ)) . Instead of computing v * , the solution of the linear system above, [43, 17] proposes a method to estimate v * . This estimation introduces bias in the gradient of the objective function.Consider the following assumptions.H6 For all θ ∈ R d , g(θ, ϕ) is strongly convex with respect to ϕ with parameter µ g > 0.H7 (Regularity Lipschitz condition) Assume that f , ∇f , ∇g, ∇ 2 g are respectively Lipschitz continuous with Lipschitz constants ℓ f,0 , ℓ f,1 , ℓ g,1 and ℓ g,2 .Assumptions H6 and H7 are the same assumptions used in [17] to obtain the convergence results with SGD. Furthermore, these two assumptions ensure that the first-and second-order derivatives of f and g, as well as the solution mapping ϕ * (θ), are well-behaved. Proposition C.1. ([17, Lemma 2.2]) Under Assumption 6, we have:EQUATION)Due to the dependence of the minimizer of the lower-level problem ϕ * (θ), obtaining an unbiased estimate of ∇V (θ) is challenging. To address this, we replace ϕ * (θ) in the gradient with ϕ and define ∇θ f (θ, ϕ):= ∇ θ f (θ, ϕ) -∇ 2 θϕ g(θ, ϕ) ∇ 2 ϕ g(θ, ϕ) -1 ∇ ϕ f (θ, ϕ) .Furthermore, by estimating ∇ 2 ϕ g(θ, ϕ)-1, we define the stochastic update H k [17] as follows:EQUATIONwhere G = N ℓg,1 N ′ i=1 I -1 ℓg,1 ∇ 2 ϕ g θ k , ϕ k+1 ; ζ (i) k with N ′ is drawn from {1, . . . ,) ∈ R d × R d , ∥∇V (θ) -∇V (θ ′ )∥ ≤ L V ∥θ -θ ′ ∥ , with the constant L V is given by L V = ℓ f,1 + ℓ g,1 (ℓ f,1 + L f ) µ g + ℓ f,0 µ g ℓ g,2 + ℓ g,1 ℓ g,2 µ g ,and L f is defined asL f = ℓ f,1 + ℓg,1ℓ f,1 µg + ℓ f,0 µg ℓ g,2 +ℓg,1ℓg,2 µg .Input: Initial points θ 0 , ϕ 0 , maximum number of iterations for the upper-level n and for the lower-level T , step sizes {γ k , γk } k≥1 , momentum parameters ρ 1 , ρ 2 ∈ [0, 1) and regularization parameter δ ≥ 0.Set m 0 = 0, V 0 = 0 and V0 = 0for k = 0 to n -1 do Set ϕ k,0 = ϕ k . for t = 0 to T -1 do ϕ k,t+1 = ϕ k,t -γk+1 ∇ ϕ g (θ k , ϕ k,t ; ζ k,t ) end for Set ϕ k+1 = ϕ k,T . Compute the stochastic update H k using ϕ k+1 . m k = ρ 1 m k-1 + (1 -ρ 1 )H k V k = ρ 2 V k-1 + (1 -ρ 2 )H k H ⊤ k Vk = max Vk-1 , Diag(V k ) A k = δI d + Vk -1/2 θ k+1 = θ k -γ k+1 A k m k end for Output: (θ k , ϕ k ) 0≤k≤nLemma C.3. Under Assumptions H6 and H7, the following inequalities hold:∥∇V (θ k ) -E [H k | F k ]∥ 2 ≤ 2L 2 f ∥ϕ k+1 -ϕ * (θ k )∥ 2 + 2 b2 k , ∇θ f (θ, ϕ) ≤ ℓ f,0 + ℓ g,1 ℓ f,0 µ g , where L f = ℓ f,1 + ℓg,1ℓ f,1 µg + ℓ f,0 µg ℓ g,2 +ℓg,1ℓg,2 µg and bk = ℓ g,1 ℓ f,11 µg 1 - µg ℓg,1 N .Proof. For the bias term, since ∇V (θ k ) = ∇f (θ k , ϕ * (θ k )), we have:∥∇V (θ k ) -E [H k | F k ]∥ 2 = ∇f (θ k , ϕ * (θ k )) -∇f (θ k , ϕ k+1 ) + ∇f (θ k , ϕ k+1 ) -E [H k | F k ] 2 ≤ 2 ∇f (θ k , ϕ * (θ k )) -∇f (θ k , ϕ k+1 ) 2 + 2 ∇f (θ k , ϕ k+1 ) -E [H k | F k ] 2 ≤ 2L 2 f ∥ϕ k+1 -ϕ * (θ k )∥ 2 + 2 b2 k ,where we used [33, Lemma 2.2] for the first term and [37, Lemma 11] for the second term.For the second inequality, we have:∇θ f (θ, ϕ) = ∇ θ f (θ, ϕ) -∇ 2 θϕ g(θ, ϕ) ∇ 2 ϕ g(θ, ϕ) -1 ∇ ϕ f (θ, ϕ) ≤ ∥∇ θ f (θ, ϕ)∥ + ∇ 2 θϕ g(θ, ϕ) ∇ 2 ϕ g(θ, ϕ) -1 ∥∇ ϕ f (θ, ϕ)∥ ≤ ℓ f,0 + ℓ g,1 ℓ f,0 µ g .Theorem C.4. Assume that H6 and H7 hold. Let θ n ∈ R d be the n-th iterate of Algorithm 4, γ n = c γ n -1/2 and γn = c γ n -1/2 /T . For any n ≥ 1, let R ∈ {0, . . . , n} be a uniformly distributed random variable. Assume the boundedness of the variance of the estimators of ∇f , ∇g, and ∇ 2 g.Then,E ∥∇V (θ R )∥ 2 = O log n √ n + b n .Proof. By using Lemma C.3, V is smooth and Lemma C.3, the bias and the gradient of V are bounded. Using our Corollary 4.5, we obtain:E ∥∇V (θ R )∥ 2 = O log n √ n + b n , where b n = O n k=0 γ k+1 b2 k + n k=0 γ k+1 ∥ϕ k+1 -ϕ * (θ k )∥ 2 √ n .Then, with [33, Lemma 2.3 ] and [17, Lemma 3], we derive:b n = O n k=0 γ k+1 b2 k √ n + 1 √ n .We achieve a classical convergence rate of O(log n/ √ n) for Stochastic Bilevel Optimization problems. Two types of bias emerge in this context: firstly, the challenge of directly computing ϕ * (θ), and secondly, the necessity of estimating [∇ 2 ϕ g(θ, ϕ)] -1 . Our results extend those of [17] to the adaptive case, particularly Adagrad, RMSProp, and AMS-GRAD. This provides convergence guarantees for the Alternating Stochastic Gradient Descent (ALSET) method. We can apply our convergence analysis to Stochastic Min-Max and Compositional Problems, as well as to the Actor-Critic method with linear value function approximation [49] , which can be viewed as a special case of the Stochastic Bilevel algorithm.We now consider a class of Conditional Stochastic Optimization:EQUATIONwhere f ξ (•) : R q → R depends on the random vector ξ and g η (•, ξ) : R d → R q is a vector-valued function dependent on both random vectors ξ and η. The inner expectation is taken with respect to the conditional distribution of η given ξ. Given certain conditions on the regularity of these functions, the gradient of V as defined in ( 14) can be expressed as:EQUATIONConstructing an unbiased stochastic estimator of this gradient can be both costly and, in some cases, impractical. Instead, we opt for a biased estimator of ∇V (θ), using just one sample ξ and m i.i.d. samples {η j } m j=1 from the conditional distribution of η given ξ:EQUATION)H8 For all ξ and η, assume that f ξ (•), ∇f ξ (•), g η (•, ξ), and ∇g η (•, ξ) are respectively Lipschitz continuous with Lipschitz constants ℓ f,0 , ℓ f,1 , ℓ g,0 and ℓ g,1 .H9 For all θ and ξ, we assume thatE η|ξ g η (θ, ξ) -E η|ξ [g η (θ, ξ)] 2 ≤ σ 2 g . Lemma C.5. ([40, Lemma 2.2]) Under Assumptions H8 and H9, the following holds:E ∇V (θ; ξ, {η j } m j=1 ) -∇V (θ) 2 ≤ ℓ 2 g,0 ℓ 2 f,1 σ 2 g m .Lemma C.6. Under Assumption H8, we have:∥∇V (θ) -∇V (θ ′ )∥ ≤ ℓ g,1 ℓ f,0 + ℓ 2 g,0 ℓ f,1 ∥θ -θ ′ ∥ , ∥∇V (θ)∥ ≤ ℓ g,0 ℓ f,0 . Proof. Denoting G θ = E η|ξ [g η (θ, ξ)] and ∇G θ = E η|ξ [∇g η (θ, ξ)],we establish the smoothness of V and Boundness of ∇V . Smoothness of V :∥∇V (θ) -∇V (θ ′ )∥ = E ξ ∇G T θ ∇f ξ (G θ ) -E ξ ∇G T θ ′ ∇f ξ (G θ ′ ) ≤ E ξ ∇G T θ ∇f ξ (G θ ) -E ξ ∇G T θ ′ ∇f ξ (G θ ) + E ξ ∇G T θ ′ ∇f ξ (G θ ) -E ξ [∇G θ ′ ∇f ξ (G θ ′ )] ≤ E ξ (∇G θ -∇G θ ′ ) T ∇f ξ (G θ ) + E ξ ∇G T θ ′ (∇f ξ (G θ ) -∇f ξ (G θ ′ )) ≤ E ξ [∥∇G θ -∇G θ ′ ∥ ∥∇f ξ (G θ )∥] + E ξ [∥∇G θ ′ ∥ ∥∇f ξ (G θ ) -∇f ξ (G θ ′ )∥] ≤ ℓ g,1 ℓ f,0 ∥θ -θ ′ ∥ + ℓ g,0 ℓ f,1 E ξ [∥G θ -G θ ′ ∥] ≤ ℓ g,1 ℓ f,0 ∥θ -θ ′ ∥ + ℓ 2 g,0 ℓ f,1 ∥θ -θ ′ ∥ .Boundness of ∇V :∥∇V (θ)∥ = E ξ ∇G T θ ∇f ξ (G θ ) ≤ E ξ [∥∇G θ ∥ ∥∇f ξ (G θ )∥] ≤ ℓ g,0 ℓ f,0 .Theorem C.7. Assume that H8 and H9 hold. Let γ n = c γ n -1/2 , A n denote the adaptive matrix in AMSGRAD and ρ 1 , ρ 2 ∈ [0, 1). For any n ≥ 1, let R ∈ {0, . . . , n} be a uniformly distributed random variable. Then,E ∥∇V (θ R )∥ 2 = O log n √ n + b n ,where b n is defined by writing m k as the number of conditional samples at iteration k:b n = O n k=0 m k √ k √ n .Proof. This is an immediate implication of Theorem 4.6 using Lemmas C.5 and C.6.These results can also be extended to the Federated Conditional Stochastic Optimization problem [75], which is defined by:min θ∈R d V (θ) = 1 L L ℓ=1 E ξ ℓ f ℓ ξ ℓ E η l |ξ ℓ g ℓ η ℓ (θ, ξ ℓ ) ,where E ξ ℓ f ℓ ξ ℓ (•) : R q → R is the outer-layer function on the ℓ-th device with the randomness ξ ℓ , and E η ℓ |ξ ℓ g ℓ η ℓ (•, ξ ℓ ) : R d → R q is the inner-layer function on the ℓ-th device with respect to the conditional distribution of η ℓ given ξ ℓ . If the functions f ℓ ξ ℓ (•) and g ℓ η ℓ (•, ξ ℓ ) for all L devices verify Assumptions H8 and H9, we obtain the same convergence rate.The following Table 2 provides a comprehensive summary of the key points, including the verification of our assumptions and the convergence results obtained in both Stochastic Bilevel Optimization and Conditional Stochastic Optimization. E ξ f (θ, ϕ * (θ); ξ) min θ∈R d E ξ f ξ E η|ξ gη (θ, ξ) s.t. ϕ * (θ) ∈ argmin ϕ∈R q E ζ [g(θ, ϕ; ζ)] Gradient ∇ θ f (θ, ϕ * (θ)) -∇ θϕ g(θ, ϕ * (θ))v * E ξ E η|ξ ∇gη (θ, ξ) T ∇f ξ E η|ξ gη (θ, ξ) Lipchitz Constant H2 ℓ f,1 + ℓ g,1 ℓ f,1 +L f µg + ℓ f,0 µg ℓ g,2 + ℓ g,1 ℓ g,2 µg ℓ g,1 ℓ f,0 + ℓ 2 g,0 ℓ f,1Bias Control H3 ℓ g,1 ℓ f,11 µg 1 - µg ℓ g,1 N ℓ 2 g,0 ℓ 2 f,1 σ 2 g m Gradient Bound H5 ℓ f,0 + ℓ g,1 ℓ f,0 µg ℓ g,0 ℓ f,0 Convergence O log n √ n + bn O log n √ n + bnIn this section, we explore examples of applications using biased gradient estimators while having control over the bias.Let π be a probability measure on a measurable space (X, X ). The objective is to estimateπ(f ) = E π [f (X)] for a measurable function f : X → R d such that π(|f |) < ∞. Assume that π(dx) ∝ w(x)λ(dx), where w is a positive weight function and λ is a proposal probability distribution, and that λ(wEQUATIONleads to the Self-Normalized Importance Sampling (SNIS) estimator:Π N f X 1:N = N i=1 ω i N f X i , ω i N = w X i N ℓ=1 w (X ℓ ), where X 1:N = X 1 , . . . , X N are independent draws from λ and the ω i N are called the normalized weights. [1] shows that the bias of the SNIS estimator can be expressed as:E Π N f X 1:N -π(f ) ≤ 12 N λ ω 2 λ(ω) 2 .This particular type of estimator can be found in the domain of Monte Carlo methods, particularly in the context of Bayesian inference and Sequential Monte Carlo methods.We focus here in the task of estimating the parameters, denoted as θ, in Hidden Markov Models. In this context, the hidden Markov chain is denoted by (X t ) t≥0 . The distribution of X 0 has density χ with respect to the Lebesgue measure µ and for all t ≥ 0, the conditional distribution of X t+1 given X 0:t has density m θ (X t , •). It is assumed that this state is partially observed through an observation process (Y t ) 0≤t≤T . The observations Y 0:t are assumed to be independent conditionally on X 0:t and, for all 0 ≤ t ≤ T , the distribution of Y t given X 0:t depends on X t only and has density g θ (X t , •) with respect to the Lebesgue measure. The joint distribution of hidden states and observations is given byp θ (x 0:T , y 0:T ) = χ(x 0 )g θ (x 0 , y 0 ) T -1 t=0 m θ (x t , x t+1 )g θ (x t+1 , y t+1 ) .Our objective is to maximize the likelihood of the model: p θ (y 0:T ) = p θ (x 0:T , y 0:T ) dx 0:T .To use a gradient-based method for this maximization problem, we need to compute the gradient of the objective function. Under simple technical assumptions, by Fisher's identity, ∇ θ log p θ (y 0:T ) = ∇ θ log p θ (x 0:T , y 0:T )p θ (x 0:T |y 0:T )dx 0:T = E x 0:T ∼p θ (.|y 0:T ) [∇ θ log p θ (x 0:T , y 0:T )]= E x 0:T ∼p θ (.|y 0:T ) T -1 t=0 s t,θ (x t , x t+1 ) ,where s t,θ (x, x ′ ) = ∇ θ log{m θ (x, x ′ ) g θ (x, y t+1 )} for t > 0 and by convention s 0,θ (x, x ′ ) = ∇ θ log g θ (x, y 0 ). Given that the gradient of the log-likelihood represents the smoothed expectation of an additive functional, one may opt for Online Smoothing algorithms to mitigate computational costs. The estimation of the gradient ∇ θ log p θ (y 0:T ) is given by:H θ (y 0:T ) = N i=1 ω i T Ω T τ i T,θ ,where {τ i T,θ } N i=1 are particle approximations obtained using particles { ξ i T , ω i T } N i=1 targeting the filtering distribution ϕ T , i.e. the conditional distribution of x T given y 0:T . In the Forward-only implementation of FFBSm [20] , the particle approximations {τ i T,θ } N i=1 are computed using the following formula, with an initialization of τ i 0 = 0 for all i ∈ 1, N :τ i t+1,θ = N j=1 ω j t m θ (ξ j t , ξ i t+1 ) ℓ=1 ω ℓ t m θ (ξ ℓ t , ξ i t+1 ) τ j t,θ + s t,θ (ξ j t , ξ i t+1 ) , t ∈ N .The estimator of the gradient H θ (y 0:T ) computed by the Forward-only implementation of FFBSm is biased. The bias and MSE of this estimator are of order O (1/N ) [20] , where N corresponds to the number of particles used to estimate it. Using alternative recursion methods to compute {τ i T,θ } N i=1 results in different algorithms, such as the particle-based rapid incremental smoother (PARIS) [64] and its pseudo-marginal extension [34] and Parisian particle Gibbs (PPG) [13] . In such cases, one can also control the bias and MSE of the estimator.Consider a finite Markov Decision Process (MDP) denoted as (S, A, R, P ), where S represents the state space, A denotes the action space, R : S × A → [0, R max ] is a reward function, and P is the transition model. The agent's decision-making process is characterized by a parametric family of policies {π θ } θ∈R d , employing the soft-max parameterization. The reward function is given by:V (θ) := E (S,A)∼v θ [R(S, A)] = (s,a)∈S×A v θ (s, a)R(s, a) ,where v θ represents the unique stationary distribution of the state-action Markov Chain sequence {(S t , A t )} t≥1 generated by the policy π θ . Let λ ∈ (0, 1) be a discount factor and T be sufficiently large, the estimator of the gradient of the objective function V is given by:H θ (S 1:T , A 1:T ) = R (S T , A T ) T -1 i=0 λ i ∇ log π θ (A T -i ; S T -i ) ,where (S 1:T , A 1:T ) := (S 1 , A 1 , . . . , S T , A T ) is a realization of state-action sequence generated by the policy π θ . It's important to note that this gradient estimator is biased, and the bias is of order O(1 -λ) [44] .Consider the problem of minimizing the objective function V . The zeroth-order gradient method is particularly valuable in scenarios where direct access to the gradient of the objective function is challenging or computationally expensive. The zeroth-order gradient oracle obtained by Gaussian smoothing [61] is given by:EQUATIONwhere τ > 0 is a smoothing parameter and X ∼ N (0, I d ) a random Gaussian vector. [61, Lemma 3] provide the bias of this estimator:EQUATION)The application of these zeroth-order gradient methods can be found in generative adversarial networks [58, 16] .The coordinate descent method is based on the iteration:θ n+1 = θ n -γ n+1 H θn (X n+1 ) jn e jn ,where {e 1 , . . . , e d } is the canonical basis of R d and H θn (X n+1 ) j is the j-th coordinate of the gradient. The randomized coordinate selection rule chooses j n uniformly from the set {1, 2, . . . , d}.Alternatively, the Gauss-Southwell selection rule [63] uses:j n+1 := argmax j∈{1,...,d} |H θn (X n+1 ) j | .This corresponds to a greedy selection procedure since at each iteration we choose the coordinate with the largest directional derivative. Another approach to choosing j n is Coordinate Sampling [53] , a variant of the stochastic gradient descent algorithm that incorporates a selection step by sampling to perform random coordinate descent. The distribution of ζ n+1 , which selects the coordinate, is characterized by the probability weights vector (wn , . . . , wn ) defined as: w (j) n = P(ζ n+1 = j|F n ), j ∈ {1, . . . , d} . This distribution of ζ n+1 is referred to as the coordinate sampling policy. The Stochastic Coordinate Gradient Descent algorithm is defined by:θ n+1 = θ n -γ n+1 D(ζ n+1 )H θn (X n+1 ) ,where D(k) = e k e ⊤ k ∈ R d×d has its entries equal to 0 except for the (k, k) entry, which is 1. Observe that the distribution of the random matrix D(ζ n+1 ) is fully characterized by the matrixEQUATIONn ). In this context, A n represents a diagonal matrix D n where the diagonal terms characterize the probability weights for sampling each coordinate. These weights typically depend on preceding iterations and even on current gradients. In this case, we always have β n+1 ≤ 1 and to control the minimum eigenvalue, we only require a lower bound on the probability weights. This method can be easily extended to incorporate biased gradients and adaptive steps by introducing Ān = D n A n , where A n represents the adaptive matrix as before, and D n is the matrix of probability weights. and Adam algorithms with a decaying learning rate, as mentioned before. For SGD, we employ the clipping method to clip the gradients to prevent excessively large steps.For this experiment, we set k = 5 samples in both IWAE and BR-IWAE, while restricting the maximum iteration of the MCMC algorithm to 5 and the burn-in period to 2 for BR-IWAE. For comparison, we estimate the Negative Log-Likelihood using these three models with SGD, SGD with momentum, Adagrad, RMSProp, and Adam, and the results are presented in Table 3 . Similar to the case of CIFAR-10, we observe that IWAE outperforms VAE, while BR-IWAE outperforms IWAE by reducing bias in all cases. The adaptive methods surpass SGD, and momentum further improves their performances. Consequently, Adam excels among all algorithms due to its adaptive steps and momentum. Similarly, as we did in the case of CIFAR-10, we incorporate a time-dependent bias that decreases by choosing a bias of order O(n -α ) at iteration n. We vary the value of α for both FashionMNIST and CIFAR-100. All figures are plotted on a logarithmic scale for better visualization and with respect to the number of epochs. The dashed curve corresponds to the expected convergence rate O(n -1/4 ) for α = 1/8, and O(log n/ √ n) for α = 1/4, as well as for α = 1/2, just as in the case of CIFAR-10. We can clearly observe that for all cases, convergence is achieved when n is sufficiently large. In the case of the FashionMNIST dataset, the bound seems tight, and the convergence rate of O(n -1/2 ) does not seem to be possible to reach, in contrast to the case of CIFAR-10 where the curves corresponding to α = 1/4 and α = 1/2 approach the O(n -1/2 ) convergence rate. For all figures, with a larger α, the convergence in both the squared gradient norm and negative log-likelihood occurs more rapidly.Additional Experiments on CIFAR-10 Dataset. The effect of C γ .Figure 8 illustrates the convergence in both the squared gradient norm and the negative log-likelihood for C γ = 0.001 and C γ = 0.01 in Adagrad. In the case of the squared gradient norm, we have only plotted the results for C γ = 0.001 for better visualization, and the plot for C γ = 0.01 was already presented in Figure 2 . It is clear that when C γ is set to 0.001, the convergence of the negative log-likelihood is slower. Similarly, the convergence in the squared gradient norm for C γ = 0.001 achieves convergence, but it is slower compared to the case of C γ = 0.01. The Impact of regularization parameter δ.In Section A.7, we discussed the impact of the regularization parameter δ in Adam. It has been empirically observed that the performance of adaptive methods can be sensitive to the choice of this parameter. Here, we illustrate the impact of this regularization parameter in IWAE. To achieve this, we plot the test loss for different sets of values for δ ∈ {10 -8 , 10 -5 , 10 -3 , 10 -2 , 5 × 10 -2 , 10 -1 } in Figure 9 . Our experimental results align with prior work [77, 67, 73] , affirming the consistent impact of δ. Notably, we find that employing δ = 5 × 10 -2 yields improved performance in IWAE.The Impact of Bias over Time.Our experiments illustrate the negative log-likelihood with respect to epochs, and we observed that a higher value of α leads to faster convergence. The key point to consider when tuning α is that while convergence may be faster in terms of iterations, it may lead to higher computational costs. To illustrate this, we set a fixed time limit of 1000 seconds and tested different values of α, plotting the test loss as a function of time in Figure 10 . It is clear that with α = 1/8, the convergence is always slower, whereas choosing α = 1/4 achieves faster convergence than α = 1/2. While the difference may seem small here, with more complex models, the disparity becomes significant. Therefore, it is essential to tune the value of α to attain fast convergence and reduce computational time.In this paper, all simulations were conducted using the Nvidia Tesla T4 GPU. The total computing hours required for the results presented in this paper are estimated to be around 100 to 200 hours of GPU usage. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.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? Answer: [Yes] Justification: The paper provides comprehensive details on the experimental setup, including hyperparameters, all optimizer algorithms, and other relevant algorithms with pseudo code, ensuring the reproducibility of the results. 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. 
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: All algorithms used to illustrate the results are run 5 times, and the mean and significance of all results are plotted, ensuring an appropriate representation of statistical significance.Guidelines:• 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. 
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: The paper provides detailed information on the computer resources required for all experiments, including the type of compute workers (CPU or GPU), and approximate time of execution.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: 
 • Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [NA] Justification: The paper does not use existing assets; therefore, no creators or original owners need to be credited, and no license or terms of use need to be mentioned or respected. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • 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. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [NA] Justification: The paper does not introduce new assets; therefore, there is no documentation provided alongside any assets. 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. https://github.com/SobihanSurendran/Adaptive-SAV k = V k ,where b is the parameter to control for achieving a better convergence rate. In this case, we have λ max (A k ) ≤ b/ log 2, which is similar to δ -1/2 in Adagrad and Adam. Additionally, they demonstrate that b ≈ 50 appears to be a good choice based on the empirical observations. 
 
 Appendix:  


Decision: This paper studies the convergence rate for a set of biased adaptive stochastic estimations (e.g., Adagrad, RmsProp), which is different from the prior works on the unbiased estimations. Most reviewers acknowledge the theoretical contributions of this paper and suggest acceptance. The authors should take the reviewers' comments (especially Reviewer nF8X) into consideration in the revision.

Reviews: review: paper_summary: paper_summary: The authors provide convergence guarantees for the biased adaptive stochastic approximation framework and establishes convergence to a critical point for non-convex setting of Adagrad type methods. The authors illustrate their results in the setting of IWAE.
strength_weakness: strength_weakness: Checking assumption H3 in its current form is not obvious for different types of noise, since it involves taking expectation w.r.t. dynamics of the process $\{\theta_n\}$ itself. With this assumption the proof of the main results is rather classical, following standard techniques in non-convex optimization. At the same time, it is not clear how H3 can be checked even for independent sequence $\{X_n\}$. Checking the same assumption for $\{X_n\}$ being a Markov chain is even more tricky due to correlations between $X_n$ and $\theta_n$. It is typical that previous papers on the subject (see e.g. [Karimi et al, 2019]) state the assumptions, that do not involve expectations over $\theta_n$. I suggest the authors to clarify H3 for independent/martingale/Markov noise separately and to clarify respective assumption (when $X_n$ is a Markov chain - probably in terms of mixing time of this sequence). Provided that my concerns regarding H3 are resolved, I will be happy to increase my score.

Also the bibliography can be complemented by the paper [Dorfman and Levy, 2022], where the authors consider adaptivity of Adagrad-type algorithms to mixing time. 

References:
[Dorfman and Levy, 2022] Dorfman, R. and Levy, K.Y. Adapting to mixing time in stochastic optimization with markovian data. In International Conference on Machine Learning 2022, pp. 5429-5446. PMLR.
[Karimi et al, 2019] Karimi, B., Miasojedow, B., Moulines, E. and Wai, H.T. Non-asymptotic analysis of biased stochastic approximation scheme. In Conference on Learning Theory, 2019, pp. 1944-1974.
questions: questions: See the weakness section.
limitations: limitations: The authors should discuss in more details the types of considered data stream $(X_n)$ and its relations with H3.
paper_summary: paper_summary: Stochastic and adaptive optimization algorithms are commonly used in advanced machine learning techniques. However, the analysis of non-convex optimization with biased gradients is lacking in the literature. This work considers the general scenario of optimizing machine learning objectives with practical optimizers given biased gradient estimators. Convergence bounds for different optimizers are proved with or without the Polyak-Łojasiewicz condition, and developed under an adaptive (preconditioning) matrix, which indicates the possibility of achieving the same convergence speed as with unbiased estimators by proper hyperparameter tuning.

The results are theoretically applied to bi-level optimization and others, empirically applied to importance weighted variational-autoencoder (IWAE) optimization. It is shown that the proved bounds apply to several applications. The experiments on Cifar-10 and FasionMNIST conclude that controlling the bias is beneficial to the convergence speed when the bias term dominates the optimizer term and the benefits become marginal after a certain threshold.
strength_weakness: strength_weakness: Despite having wide applications, the experiments only include one case that verifies the proved bounds. The paper spent some efforts in bi-level optimization, but did not provide empirical analysis. IWAE has a remote connection to bi-level optimization when q is flexible enough, which could be utilized with a different optimization scheduling. 

minor: line 135, there seems to be a redundant factor of 1/2 for $r_{n+1}$
questions: questions: The variational distribution q may also affect the convergence. In the experiments, is q fixed or trained? I think the analysis applies to fixed q, but would like to learn about assumptions and implementation details here.

The introduction also talked about Markov chain based optimization schemes, but the analysis later in the paper seems to be far from it. How is assumption H3 achieved in such case? How are the theorems connected to the Martingale and Markov chain cases?

The main analysis in the paper is limited to optimization with vanishing bias, which does not directly apply to some practical methods (e.g. relaxation with Gumbel-softmax/concrete distribution, variational sequential Monte Carlo that drops part of the gradient). Can you provide a useful bound when r_n is a (or only bounded by a) constant from Theorem 4.2?
limitations: limitations: I do not see a limitation section and it would be beneficial to sketch future works by discussing limitations.
paper_summary: paper_summary: This paper studies the non-asymptotic convergence guarantees of SGD with adaptive step sizes and (time-dependent) biased gradient estimators for nonconvex smooth functions. Applications to AdaGrad, RMSProp and Adam are developed. Numerical experiments on bilevel and conditional stochastic optimization, as well as deep VAE are used to illustrate the established theoretical results.
strength_weakness: strength_weakness: The bounded spectrum assumption of $A_n$ might not hold for every optimizer considered, and it remains to see how to relax it.
questions: questions: Sometimes “adaptive step sizes” is referred to as “adaptive steps” in the paper. I guess using “adaptive step sizes” for all instances will avoid possible confusion. Can you also address the above point in weaknesses?
limitations: limitations: Yes.
paper_summary: paper_summary: This paper proposes non-asymptotic convergence guarantuees on several gradient-based optimization methods, ranging from the SGD to AdaGrad-like methods, when the estimation of the gradients is biased. More specifically, the convergence results include a theorem with the Polyak-Lojasiewicz condition and another one without the PL condition, both in a non-convex smooth setting.
strength_weakness: strength_weakness: # Signifiance

How limitating are hypotheses H1-5? It is difficult to check when the proposed models fulfill them.
questions: questions: How limitating are hypotheses H1-5? It is difficult to check when the proposed models fulfill them.
limitations: limitations: -
paper_summary: paper_summary: This paper aims to analyze SGD with biased gradients in a non-asymptotic manner, where the steps  are also adaptive. In particular, under certain assumptions and conditions, it provides convergence guarantees and establishes convergence rates for a critical point of non-convex smooth functions for various adaptive methods.  

Overall, I think the theoretical results of this paper are interesting and insightful.
strength_weakness: strength_weakness: One thing unclear to me is the role of parameters that are related to sampling noise (e.g. batch size) in the theorems. These parameters are supposed to be crucial since they are actually closely related to the training process and generalization. How do theorems in this paper characterize the effects of such parameters (other than some simple bounds from the assumptions)?
questions: questions: - Please see the weaknesses part.

- In addition, although I understand that the settings considered in this paper are novel, there are already a lot of works in closely related topics. Does there exist any technical difficulty that cannot be solved by techniques from previous works? If it does, what are the technical difficulties compared to previous works?
limitations: limitations: The authors have adequately addressed the limitations in the main part of the paper.
paper_summary: paper_summary: This paper studied convergence rate of general adaptive stochastic approximation with biased gradient. This is an important field since only biased gradient is accessible in many practical machine learning problem. The authors establied non-asymptotic bound based on various assumption, some of which are quite strong.
strength_weakness: strength_weakness: 1. The techinque is not novel compared with existing works on biased SGD, e.g. [1]. The authors directly impose an upper and lower bound of preconditioning matrix in H3, in which case the analysis is not very different from that of SGD. A better way is to establish the boundedness of preconditioning matrix along the tracjectory as in e.g., [2], instead of taking it as prior assumption.

2. The authors claimed that they analyzed Adam throught out the paper, but in fact they only did AMSGRAD. This is misleading since they are two totally different algorithms with distinct practical performance. I suggest the authors correct this point.

### Typos:

line 230: $\mathcal{O}(n^{-1/4})$->$\mathcal{O}(\log n/n^{1/2})$

[1] Liu, Yin, and Sam Davanloo Tajbakhsh. "Adaptive stochastic optimization algorithms for problems with biased oracles." arXiv preprint arXiv:2306.07810 (2023).

[2] Li, Haochuan, Alexander Rakhlin, and Ali Jadbabaie. "Convergence of adam under relaxed assumptions." Advances in Neural Information Processing Systems 36 (2024).
questions: questions: In Theorem 4.2, what does $\delta_j$ stand for? Does it correspond to $\delta$ in eq (4)? What is the dependence of $\delta$ in the final convergence rate of Adagrad and other algorithms?
limitations: limitations: The limitations are pointed out in the weakness part.

Comments: title: Response to Rebuttal
Thanks for the clarification
comment: Thanks so much to the authors for your rebuttal. Your rebuttal has addressed my concern.
Thank you for your comments. Regarding Adagrad and RMSProp, although the analysis may appear similar to that of biased SGD under Assumption H3(i), verifying H3(i) requires managing the form $A_n$, which differs from the approach used in biased SGD, which is detailed in line 573. Additionally, the analysis of AMSGrad differs significantly from both the general case and the specific applications to Adagrad and RMSProp.

The proof for AMSGrad is derived from scratch and does not rely on assumptions about the boundedness of eigenvalues. However, we do assume that gradients are bounded (as per Assumption H5). Even with this assumption, the proof for AMSGrad is different from the analysis of biased SGD, as detailed in Appendix A.6.

While Assumption H5 is common in adaptive algorithm literature [8, 64, 70, 71], recent works have provided theoretical guarantees for unbiased gradients without relying on the boundedness of gradients (e.g., [82]).
Our work focuses on adaptive algorithms with biased gradients, including scenarios with decreasing bias. We have shown that in various applications, such as Stochastic Bilevel Optimization and Conditional Stochastic Optimization, all our assumptions, including H3 and H5, can be verified. 
Most results that do not rely on Assumption H5 are provided with high probability. Furthermore, the obtained bounds generally have a rate of $\mathcal{O}\left(\text{poly}(\log n) / \sqrt{n} \right)$, with a loss of a logarithmic factor in the polynomial.
However, even some recent works on SGD and adaptive methods continue to rely on the assumption of bounded gradients (e.g., see [1, 2]).

One approach to avoid assumption H5 is to use the truncation method. To control the smallest eigenvalue in Adagrad, we can consider a truncation that leads to

$$
A_{n} = \left( \text{diag} \left( \beta_{n+1}^{-2}I_{d} + \frac{1}{n+1}\sum_{k=0}^{n}H_{\theta_{k}} \left( X_{k+1} \right) H_{\theta_{k}}\left( X_{k+1} \right)^{T} \mathbf{1}\_{|| H_{\theta_{k}} \left( X\_{k+1} \right) ||^{2} \leq v\_{k+1}} \right)   \right)^{-1/2} ,
$$
so that 
$$
\lambda_{min} \left( A_{n} \right) \geq  \left(\beta_{1}^{-1} + \frac{1}{\sqrt{n+1}}\sqrt{\sum_{k=0}^{n}v_{k+1}} \right)^{-1} .
$$ 
Then, choosing $v_{k} = C_{v}k^{-v}$ with $v \geq 0$ leads to 
$$
\lambda_{\min} \left( A_{n} \right) \geq \frac{1}{\beta_{1}^{-1} + \sqrt{C_{v}} (n+1)^{v/2}} =: \lambda_{n+1}.
$$
The value of $v$ can be chosen arbitrarily to achieve the desired convergence rate. However, the larger the chosen $v$, the less truncation is applied, which may result in a slower convergence rate. It is also possible to adjust other parameters, such as $\gamma_{n+1}$ and $\beta_{n+1}$.

While there are a few other works on the scalar version of Adagrad with biased gradients, their assumptions are stronger than ours.
Our work is the first step to providing the general convergence results for the biased gradients and adaptive algorithms for any form of $A_n$, which applies to Adagrad and RMSProp and extends to AMSGRAD under the assumption of bounded gradients.

Additionally, at the suggestion of Reviewer ioBT and Reviewer BcqM, we have provided convergence results for the Martingale and Markov Chain cases. Our work focuses on deriving convergence results for various adaptive algorithms based on commonly used assumptions in the literature, applicable to many different scenarios. We believe that our paper is a valuable contribution to the community, and that exploring convergence results for biased gradients without relying on Assumption H5 may require a separate paper.
We plan to discuss the boundedness of gradients along the trajectory and will address this in future work.


[1] Xidong Wu, Jianhui Sun, Zhengmian Hu, Junyi Li, Aidong Zhang, and Heng Huang. Federated conditional stochastic optimization. In Advances in Neural Information Processing Systems, volume 36, 2024.

[2] Xidong Wu, Feihu Huang, Zhengmian Hu, and Heng Huang. Faster adaptive federated learning. In Proceedings of the AAAI conference on artificial intelligence, volume 37, pages 10379–10387, 2023.
Thank you for detailed answer. I will increase my score to 6.
I thank the authors for the detailed response which addressed my concerns and questions. I keep my positive score unchanged.
Thanks for the reply. I would like to clarify that regarding w1, I notice in the analysis of Adagrad, RMSProp and AMSGrad, the Assumption H5 is used. This assumption, along with the existence of $\delta$, directly imposes an upper and lower bound on the preconditioning matrix globally. Hence I don't think the analysis will be too much different from biasd SGD. This is different from [2] which establishes the boundedness along the tracjectory.
Thank you for the rebuttal which addresses my questions. After reading the other reviews and responses, I do not have any additional concerns.