Conceptio › Archive › arXiv CS

arXiv CSopen access

The Bias of Nonlinear Two-Time-scale Stochastic Approximation under Constant Step-Sizes

· arxiv_cs

arXiv CS · Papers · License: Open Access

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

The Bias of Nonlinear Two-Time-scale Stochastic Approximation under Constant Step-Sizes

arXiv:2609.20409v1 [cs.LG] 17 Sep 2026

Djamel Rassem Lamouri

Dorian Baudry

Nicolas Gast

Univ. Grenoble Alpes, CNRS, Inria Grenoble INP, LIG, 38000 Grenoble France

Abstract Two-timescale stochastic approximation (TTSA) is a fundamental tool for analyzing coupled iterative algorithms in reinforcement learning, optimization, and stochastic control. However, finite-time guarantees for nonlinear two-timescale schemes remain difficult to obtain, especially under constant step-sizes. In this paper, we study nonlinear TTSA with step-sizes α ≫ β. Under standard stability, regularity, and Markovian noise assumptions, we upper bound the mean-squared error and the bias of both iterates around their limiting equilibria. Our bounds scale as O(α + β 2 /α2 ), which we prove to be tight when β ≤ α3/2 . The analysis separates the contributions of initial conditions, fast-timescale tracking error, Markovian dependence, and timescale coupling, thereby clarifying the origin of the β 2 /α2 term. Our results reveal qualitative differences from the linear TTSA setting previously studied, showing that nonlinear dynamics introduce additional finite-time effects that are absent in the linear case.

Introduction

A broad family of machine learning and optimization algorithms can be formulated within the framework of Stochastic Approximation (SA) [4, 25, 3]. In the two-timescale setting [24, 29], two sequences of parameter vectors, {xk }k≥0 and {yk }k≥0 , are generated from the following recursive update rule: xk+1 = xk − αk f (xk , yk , ξk+1 ) (1) yk+1 = yk − βk g(xk , yk , ξk+1 ), where, for a measurable space (Ξ, S), the S-measurable functions f : Rdx × Rdy × Ξ → Rdx and g : Rdx × Rdy × Ξ → Rdy are the noisy dynamics of the system, {ξk }k≥1 is the source of noise modeled as a stochastic process defined on Ξ and {αk }k≥0 , {βk }k≥0 are positive sequences representing the step-sizes. Some notable examples of algorithms that naturally exhibit this scheme of recursive updates include AdaGrad [13], Generative adversarial networks [16], Bi-level Optimization [20] and Reinforcement Learning [34, 23]. To study the dynamics of SA, a classical tool is to use the ODE method [5]. It consists in viewing (1) as a noisy Euler discretization of the singularly perturbed ODE 1 ẋ(t) = − f (x(t), y(t)) ϵ ẏ(t) = −g(x(t), y(t)) ,

(2)

with ϵ = βk /αk ≈ 0 [4, Chapter 6]. The small parameter ϵ induces a separation between the two timescales: The variable x evolves on the fast timescale, whereas y evolves slowly and is Preprint.

therefore seen by the x-dynamics as quasi-static. This motivates the standard decoupling of the two ODEs: for each fixed y, the fast dynamics are described by ẋ(t) = −f (x(t), y). Under suitable stability assumptions, the trajectory solutions of this fast ODE converge to the asymptotically stable equilibrium h(y) satisfying f (h(y), y) = 0. Once this tracking property is established, the slow ODE is ẏ(t) = −g(h(y(t)), y(t)). Under suitable assumptions, the solution of the ODE converges to a point y ∗ . The ODE method consists in assuming that xk ≈ h(yk ) and yk ≈ y ∗ . To guarantee the accuracy of the ODE method, it is therefore important to characterize the error that is made when approximating (1) by the solutions of the ODE (2). The classical convergence results for SA and TTSA typically rely on diminishing step-sizes. For two time-scale schemes, one usually assumes the step-sizes to have infinite ℓ1 -norm, finite ℓ2 -norm and that limk→∞ αβkk = 0 which ensures the time-scale separation of the associated limiting ODEs [4]. In many applications, however, constant step-sizes αk ≡ α and βk ≡ β are preferred, notably because they allow an exponential rate of forgetting initial conditions [30]. In this regime, one should not expect exact convergence to the equilibrium; instead, the iterates typically approach a neighborhood whose size depends primarily on the step-sizes and the noise. This motivates a finite-time analysis of the constant-step-size TTSA under Markovian noise. Previous works have quantified the bias induced by constant step-sizes for nonlinear single time-scale SA [1] and for linear TTSA [19]. The nonlinear TTSA setting remains more delicate, because the fast tracking error and the slow dynamics interact through the nonlinear manifold x = h(y). Our goal is to characterize this interaction quantitatively. We will measure the performance of the algorithm through two quantities: the mean-squared error (MSE), which is defined as h i h i 2 2 MSExk ≜ E ∥xk − h(yk )∥ and MSEyk ≜ E ∥yk − y ∗ ∥ , and the bias of the iterates around their limiting equilibrium, which we formally define as Biasxk ≜ E [xk ] − h(y ∗ )

and

Biasyk ≜ E [yk ] − y ∗ .

Our goal in this paper is to provide a finite-time analysis of nonlinear two time-scale stochastic approximation with constant step-sizes under a quite general Markovian noise. One original feature of our model is that we allow the kernel of the noise process ξk to depend on the iterates (xk , yk ). This is particularly relevant for applications to reinforcement learning where the exploration policy might depend on the current values of the parameters. Contributions. Our first main result, Theorem 1, proves that the mean-squared errors of both time-scale components are bounded, up to exponentially decaying initialization terms, by O α + β + β 2 /α2 , which reduces to O(α + β 2 /α2 ) when β ≤ α. The proof of Theorem 1 combines a contraction argument for both iterates xk and yk with a decomposition of the Markovian noise based on the problem-specific Poisson equation which allows us to handle the dependencies of the transition kernel on the model’s parameters. The bound we prove separates the exponentially vanishing contribution of the initialization from the persistent error induced by constant step-sizes. The resulting recursion makes explicit the different sources of error: the intrinsic fluctuations of the fast iterate, the slow motion of the quasi-stationary target h(yk ), and the bias induced by the noise. This extends the approaches developed for nonlinear single time-scale SA [1] and nonlinear TTSA with diminishing step-sizes [10] to the constant-step-size non-linear TTSA regime. When β ≤ α3/2 , our MSE bounds are of order O(α) for both the fast and the slow iterate. In Theorem 2, we show that this bound cannot be improved by exhibiting an example whose MSE is of the order Ω(α). This shows that the non-linear TTSA is significantly different from the linear TTSA for which the MSE of the slow iterate was shown to be O(β + α2 ) in [26]. Our second main result, Theorem 3, establishes matching-order bounds on the bias of the fast and slow iterates. The bias analysis is substantially more delicate than the MSE analysis because the analysis of the first-order dynamics of the two bias components remains coupled. We handle this difficulty by linearizing the nonlinear drift around the limiting equilibrium and then applying a sequence of coordinate transformations that separates the dominant fast and slow modes. In these transformed coordinates, the bias satisfies a recursion formula, which leads to upper bounds involving the Markovian-noise residuals and the MSE bounds fromTheorem 1. Transforming back to the original variables yields bias bounds of order O α + β 2 /α2 . When β ≤ α3/2 , this again reduces to O(α) which we prove to be tight in Theorem 4. 2

Outline. The remainder of the paper is organized as follows. In Section 2, we detail related work on single and two time-scale stochastic approximation, focusing on finite-time guarantees obtained with constant step-sizes under Markovian noise. Section 3 introduces the nonlinear TTSA model and states the assumptions used throughout the paper. Section 4 presents Theorem 1, a finite-time MSE bound, together with a proof overview and a discussion of tightness. Section 5 presents Theorem 3, a finite-time bias bound, and outlines the main ideas behind its proof. Finally, Section 6 illustrates the theoretical bounds by creating an example of nonlinear TTSA showing how the term β 2 /α2 persists in the bounds of our MSE.

Related Work

The study of stochastic approximation (SA) algorithms traces back to the seminal work of [33], which introduced a framework for finding a root x that satisfies M (x) = b for a given b ∈ R. Here, when referring to SA, we mean one time-scale stochastic approximation, where only the variable x is present (as opposed to to TTSA where two variables are present and evolve at different time-scales β ≪ α). In this classical setting, the monotonic function M : R → R is assumed to be unknown, requiring the algorithm to rely strictly on noisy observations. A pivotal development in the theoretical understanding of SA is through the ordinary differential equation (ODE) method [5]. This technique characterizes the asymptotic convergence properties of discrete SA updates by analyzing their continuous-time limit, represented by the analogous ODE ẋ = −f¯(x). The ODE method is the method of choice for analyzing SA dynamics, with its underlying theory and applications extensively detailed across numerous foundational texts [4, 25, 3]. Constant step-size SA. A line of work has been developed around the constant step-size assumption in SA due to practical advantages in algorithm design. Starting with Linear SA (LSA), [27] have shown that the MSE of the Polyak-Ruppert (PR) averages [31] decreases with order O( k1 ) under i.i.d noise. Previously, in [2] averaging was used on non-strongly convex SA to achieve the same rate. Under proved that with probability at least 1 − δ, |xk − x∗ | is bounded by p the same settings [15] ∗ O( α log(1/δ)), where x denotes the equilibrium of x. Later, the work of [14] generalized the high probability bounds of the latter result to PR averages under Markovian noise. In [19], it is shown that LSA under constant step size and Markovian noise is characterized by a bias of order O(α) that can be improved to higher orders by applying Richardson-Romberg extrapolation. Switching to non-linear SA, [7] showed an MSE of order O(α log(1/α)). This bound is improved in [1] to O(α), with a characterization of the bias vector by αV + O(α2 ), where V is a solution vector of a problem-specific Lyapunov equation. Linear TTSA. In the linear regime, the foundational work [24] established the asymptotic variance and normality of the iterates under i.i.d. noise and decreasing step-sizes. Subsequent research shifted toward finite-time analyses. In [21], the MSE of the fast iterate was shown to be of order O(αk ) and that of the slow iterate O(βk ) for step-sizes satisfying αk = O(k −v ) for some v < 1 and βk = O(k −1 ). A convergence rate of O(1/k 2/3 ) under Markovian noise was derived in [9], which was later tightened to O(1/k) in [17]. The constant step-sizes setting was studied in [26], which established an upper bound of order O(α) and O(β + α2 ) for the fast and slow iterates, respectively. They further proved upper bounds of order O(α + β) for the bias of both iterates. Our results indicate that these bounds are specific to this setting, since they degrade after removing the linearity assumption. Nonlinear TTSA. Parallel advancements have been made in non-linear TTSA. Extending asymptotic normality to this setting, [29] proved a Central Limit Theorem (CLT) for the iterates, a property recently expanded upon by [18] under decreasing step-sizes. Regarding finite-time guarantees, [10, 9] achieved an MSE rate of O(1/k 2/3 ) under both i.i.d. and Markovian noise. More recently, [6] analyzed contractive non-linear TTSA under Markovian noise, demonstrating that the MSE for both iterates is bounded by O(αk + βk2 /αk2 ). All of those works assume decreasing step-sizes. The remaining gap, which we address in this work, is to obtain finite-time bounds on both the MSE and the bias of nonlinear TTSA with constant step-sizes, under Markovian noise. 3

Assumptions

First, to guarantee algorithmic stability, we assume that the algorithm’s iterates remain within a compact set. This is often enforced in practice via explicit projections [25, Chapter 5], [28]. A. 1. There exist two compact sets X ⊂ Rdx and Y ⊂ Rdy such that for all k ≥ 0, xk ∈ X and yk ∈ Y a.s. Our second assumption concerns the stochastic process (ξk )k≥1 , which we assume to be a positive recurrent Markov chain. This allows us to account for temporally correlated data, in contrast with the classical i.i.d. or martingale-difference noise models, and is particularly relevant for applications such as reinforcement learning. The price of this added flexibility is that the analysis must handle the bias and dependence induced by the Markovian dynamics. A. 2. (ξk )k≥1 evolves as a Markov chain on a finite state space Ξ with a transition kernel K(xk ,yk ) that is allowed to depend on xk , yk , that is: for every ξ ′ , ξ ∈ Ξ and x ∈ X , y ∈ Y we have: P(ξk+1 = ξ ′ |ξk = ξ, Fk ) = P(ξk+1 = ξ ′ |ξk = ξ, xk = x, yk = y) ≜ K(x,y) (ξ, ξ ′ ), where {Fk }k≥1 is the natural filtration of the process (xk , yk , ξk )k≥1 . We further assume that the kernel K(x,y) corresponds to a unichain Markov chain for all (x, y) ∈ Rdx × Rdy . Also, we assume that for all ξ, ξ ′ ∈ Ξ the quantity K(.,.) (ξ, ξ ′ ) is LK -Lipschitz continuous: K(x1 ,y1 ) (ξ, ξ ′ ) − K(x2 ,y2 ) (ξ, ξ ′ ) ≤ LK (∥x1 − x2 ∥ + ∥y1 − y2 ∥). The assumption of having a Markovian noise is standard in stochastic approximation and covers a broad class of algorithms [26, 18, 10, 7]. The fact that we allow the kernel at time k to depend on (xk , yk ) is more original and is relevant in reinforcement learning algorithms where the exploration policy might depend on the current parameters [1]. Assumption A.2 guarantees the existence of a unique stationary distribution π(x,y) for each (x, y), thereby allowing us to define the averaged functions f¯ and ḡ by integrating f and g with respect to π: X f (x, y) ≜ E[f (x, y, ξ)] = π(x,y) (ξ)f (x, y, ξ) and g(x, y) ≜ E[g(x, y, ξ)] . (3) ξ∼π(x ,y )

ξ∼π(x ,y )

ξ∈Ξ

Furthermore, the assumption on Lipschitzness of K implies that the stationary distribution π(x,y) is also Lipschitz-continuous in (x, y), see e.g. [1]. Third, to ensure the existence and uniqueness of the limiting trajectories of the equivalent ODE, we need to assume Lipschitzness of f, g, which, combined with Lipschitzness of K and the definitions of f¯ and ḡ of (3) implies that the functions f¯ and ḡ are also Lipschitz continuous which is a fundamental prerequisite of the ODE method, ensuring the existence and uniqueness of solutions of the ODE. We start with our assumptions on f : A. 3. For all ξ ∈ Ξ, the function f (., ., ξ) is Lf -Lipschitz continuous: ∥f (x1 , y1 , ξ) − f (x2 , y2 , ξ)∥ ≤ Lf (∥x1 − x2 ∥ + ∥y1 − y2 ∥) In addition, there exists µf > 0 such that f is µf -strongly monotonic in its first argument, i.e., for every x1 , x2 ∈ X and y ∈ Y we have ⟨f (x1 , y) − f (x2 , y), x1 − x2 ⟩ ≥ µf ∥x1 − x2 ∥2 . This strong monotonicity assumption guarantees that for a fixed y, the ODE ẋ(t) = −f¯(x(t), y) possesses a unique, exponentially stable equilibrium, which we denote by h(y). The equilibrium map L h is µff -Lipschitz from Lipschitzness and strong monotonicity of f : 2 µf ∥h(y1 ) − h(y2 )∥ ≤ f¯(h(y1 ), y1 ) − f¯(h(y2 ), y1 ), h(y1 ) − h(y2 ) = f¯(h(y2 ), y2 ) − f¯(h(y2 ), y1 ), h(y1 ) − h(y2 )

≤ f¯(h(y2 ), y2 ) − f¯(h(y2 ), y1 ) ∥h(y1 ) − h(y2 )∥ ≤ Lf ∥y2 − y1 ∥ ∥h(y1 ) − h(y2 )∥ L

The result follows when h(y1 ) ̸= h(y2 ), we denote Lh ≜ µff . Moving to the slow iterate, we need Lipschitzness of g: 4

A. 4. For all ξ ∈ Ξ, the function g(., ., ξ) is Lg -Lipschitz continuous: ∥g(x1 , y1 , ξ) − g(x2 , y2 , ξ)∥ ≤ Lg (∥x1 − x2 ∥ + ∥y1 − y2 ∥) In addition, for µg > 0, g is 1-point µg -strongly monotonic with respect to a point y ∗ ∈ Y that satisfies ḡ(h(y ∗ ), y ∗ ) = 0: for all y ∈ Y: ⟨g(h(y), y), y − y ∗ ⟩ ≥ µg ∥y − y ∗ ∥2 . This assumption implies that the function ḡ(h(y), y) is also Lipschitz continuous and in particular implies that the ODE ẏ = −ḡ(h(y), y) has a unique solution. To guarantee that this ODE has a unique stationary point y ∗ , we added the 1-point strong monotonicity at y ∗ , for instance, as in [11, 8]. It can be viewed as a generalization of the non-linear case of having Hurwitz matrices in the linear setting [26, 19, 15], which are the conditions required for the ODE trajectories to exhibit global convergence. Pd Notation Let x ∈ Rd , ∥x∥ is the Euclidean norm ∥x∥ = ( i=1 x2i )1/2 . For y ∈ Rd we denote the inner product between x and y by ⟨x, y⟩. For a square matrix M ∈ Rn×n we define the operator norm ∥M ∥op = sup∥x∥=1 ∥M x∥. For a positive definite matrix L we define the condition number (L) of L as κ(L) = λλmax where λmax (L) and λmin (L) are the largest and smallest eigenvalues of min (L) L. For a multivariate function f¯ : Rdx × Rdy → Rdx , ∇x f (x, y) denotes the partial Jacobian with i (x,y) respect to the first argument (∇x f (x, y))i,j = ∂f∂x and ∇y f (x, y) denotes the partial Jacobian j with respect to the second argument.

The MSE Analysis of TTSA

Our first main result is a bound on the Mean Squared Error (MSE) for both the fast and the slow iterates. 4.1

Main Result

The following theorem provides non-asymptotic bounds on MSExk and MSEyk , explicitly separating the transient initial condition from the steady-state residual error caused by the Markovian noise and the constant step-sizes. Theorem 1. Assuming that A.1-3 hold, and that αµf < 1 and β ≤ α, then for any k ≥ 1 we obtain the following upper bound on the MSE of the fast iterate, B1 β2 2 2 k E[∥xk − h(yk )∥ ] ≤ (1 − αµf ) ∥x0 − h(y0 )∥ + α + B2 + β B3 + 2 B4 , (4) 1 − αµf α | {z } | {z } Decaying initial condition MSEx ∞ (α,β)

where the constants (Bk )k∈[4] are fully defined in appendix (Eq. (18)-(21)), and are independent of α, β and k: they only depend on the properties of the Markovian noise, the range, monotonicity and Lipschitzness parameters of the functions f, g, h. Further, assuming A.4, and βµg < 1, for any k ≥ 1 the MSE of the slow iterate is upper bounded by 2

E[∥yk − y ∗ ∥2 ] ≤ (1 − βµg )k ∥y0 − y ∗ ∥ + D1 β k(1 − (αµf ∧ βµg ))k ∥x0 − h(y0 )∥ | {z } Decaying initial conditions

+ D2 MSEx∞ (α, β) + β |

D3 + D4 1 − βµg {z

MSEy ∞ (α,β)

+ αD5 . }

where the constants (Dk )k∈[5] are also formally defined in appendix (Eq. (28)-(32)) and depend again on problem parameters but not on α, β, k. Proof Overview. To prove the MSE bounds for the fast iterate xk , we first use the strong monotonicity assumption A.3 to extract the following recursion on the tracking errors between steps k and k + 1, 2

∥xk+1 − h(yk+1 )∥ ≤ (1 − αµf ) ∥xk − h(yk )∥ + rkx , 5

(5)

where rkx is a residual term containing the deviation of the noisy dynamics from their stationary means in addition to higher order terms of α and β that we detail in the next paragraph. The goal of this recursion is to dissociate the impact of the initial conditions from the cumulative impact of the noise. Repeatedly applying the above inequality, we get 2

∥xk+1 − h(yk+1 )∥ ≤ (1 − αµf )k+1 ∥x0 − h(y0 )∥ +

k X (1 − µf α)k−j rjx . j=0

This first step is similar to the approach used for the decreasing step sizes in [10], except in the latter case the sum of the residuals under expectation vanishes as k → ∞ which is not the case here due to the assumption of fixed step-sizes. To bound this sum of residuals, we then use A.3 to show that each rjx is bounded by 2α xj − h(yj ), f (xj , yj , ξj ) − f (xj , yj ) , plus a term Mα,β that is of order O(α2 + β 2 + αβ + Pk β 2 /α). The sum j=0 (1 − µf α)k−j Mα,β leads to a O(α + β + β 2 /α2 ) error term. The main technical difficulty is therefore to control the sum induced by the first term. To do so, we follow the Poisson-equation technique used in [1, Lemma 8] to rewrite the noise as a Martingale difference term, a bounded term that scales linearly with the step-sizes, and a weighted telescoping term. Because of the geometric weights (1 − αµf )k−j , the telescoping term is not exact and leaves additional terms. Additionally, since the kernel is depending on the current iterates, the Martingale term creates a drift in the kernel that we control due to Lipschitzness of K from A.2. Lemma 2 shows that all these contributions are of the required order. To prove the MSE bound for the slow iterate, we apply a recipe similar to the one followed for the fast iterate. From assumption A.4 we obtain an equation for y similar to (5) with a contraction factor 1 − βµg , along with an additional term capturing the influence of the fast iterate: 2

∥yk+1 − y ∗ ∥ ≤ (1 − βµg ) ∥yk − y ∗ ∥ + β ∥xk − h(yk )∥ + rky . 2

We then use the results of the fast iterate to bound ∥xk − h(yk )∥ . The treatment of rky is then similar to the analysis done for the fast iterate. The full proof can be found in Appendix B. 4.2

Discussion and tightness

Under the constant step-size regime and the strong monotonicity assumptions A.3 and A.4, the initial conditions are forgotten at an exponential rate, respectively (1 − αµf )k for the fast iterate and (1 − βµg )k for the slow iterate. This illustrates one of the main practical advantages of constant step-sizes: the iterates rapidly enter a neighborhood of the equilibrium, although the size of this neighborhood is controlled by the persistent stochastic error. The two time-scale structure also creates a transient coupling between the initial fast error and the slow iterate, which appears through the term k 2 k 1 − (αµf ∧ βµg ) ∥x̂0 ∥ , which may initially increase because of the prefactor k, but still vanishes exponentially as k → ∞. Once the transient terms have vanished, the iterates remain, in expectation, within a neighborhood of the equilibrium whose squared radius scales as O α + β 2 /α2 . The presence of the term β 2 /α2 in both bounds is one of the main features of Theorem 1. It is consistent with finite-time results for nonlinear TTSA with decreasing step-sizes, where bounds of order O(αk + βk2 /αk2 ) are obtained for both iterates [6]. This term quantifies the cost of coupling the fast and slow recursions: without sufficient time-scale separation, it may dominate the steady-state error. However, when the slow step-size is sufficiently small relative to the fast one, the α-term becomes dominant. In particular, if β ≤ α3/2 , then β 2 /α2 ≤ α, and Theorem 1 immediately yields that, for any TTSA satisfying its assumptions, there exist constants C > 0 and α0 > 0 such that, for all α < α0 and all β ≤ α3/2 , lim sup E[∥xk − h(yk )∥2 ] ≤ Cα,

and

k→∞

lim sup E[∥yk − y ∗ ∥2 ] ≤ Cα.

(6)

k→∞

We emphasize that the steady-state MSE bound for the slow iterate is substantially larger than what is obtained in the linear TTSA setting, where the corresponding rate is O(β + α2 ) [26]. This raises a 6

natural question: is the O(α) term in the bound for the non-linear case unavoidable, or is it merely a proof artifact? The next result shows that it is not. It constructs a nonlinear TTSA instance satisfying the assumptions of Theorem 1 for which both the fast tracking error and the slow error have asymptotic MSE at least of order α. Thus, in the nonlinear setting, the rate in (6) cannot be improved in general: there is a qualitative separation from the linear TTSA setting, where sharper rates are possible for the slow iterate. Theorem 2. There exist constants C > 0, α0 > 0, and a TTSA instance satisfying all the assumptions of Theorem 1 such that, for all step-sizes β ≤ α < α0 , lim inf E ∥xk − h(yk )∥2 ≥ Cα, and lim inf E ∥yk − y ∗ ∥2 ≥ Cα. k→∞

k→∞

The proof is given in Appendix C and relies on a simple one-dimensional construction for both variables xk and yk .

The Bias Analysis of TTSA

Analyzing the bias requires a finer characterization of the dynamics. We therefore assume additional differentiability, which allows us to linearize the averaged update around the equilibrium. The first-order terms drive the bias recursion, while the second-order Taylor remainders are controlled by the MSE bounds from Theorem 1. We use the following notation: Jxx ≜ ∇x f (h(y ∗ ), y ∗ ), Jxy ≜ ∇y f (h(y ∗ ), y ∗ ), Jyx ≜ ∇x g(h(y ∗ ), y ∗ ) and Jyy ≜ ∇y g(h(y ∗ ), y ∗ ). B. 1. f and g are at least two times differentiable and their second derivatives are bounded operators. −1 In addition, the matrices −Jxx and −∆ ≜ −Jyy + Jyx Jxx Jxy are Hurwitz. The differentiability part of B.1 is standard in analyses that rely on a local expansion of the stochasticapproximation dynamics. It is used, for instance, to characterize the constant-step-size bias in nonlinear SA [1], and it is also a natural assumption in CLT analyses for nonlinear TTSA [18]. In our proof, this assumption allows us to linearize the averaged dynamics around the equilibrium: the first-order terms determine the leading bias recursion, while the second-order remainders are controlled through the MSE bounds. The Hurwitz condition on −Jxx and −∆ is essential when studying TTSA settings [26, 17, 18, 21, 24]. In particular, for the linear settings, it is equivalent (up to a linear transformation) to monotonicity in A.3 and A.4. It is needed to ensure that the trajectories under the approximate linear model converge exponentially fast to the equilibrium. 5.1

Main result

We now introduce the second main result of this paper, which provides an upper bound on the bias of TTSA. The result shows that, once the initial conditions vanish, the bias is of the same order as the bound on the MSE obtained in Theorem 1. Theorem 3. Assuming A.1-4 and B.1, there exists α∗ defined in (8), β ∗ defined in (9), and ω ∗ defined β in (10) such that for α < α∗ , β < β ∗ and α < ω ∗ we have the following bound on the bias of both iterates: max{∥E [xk − h(y ∗ )]∥ , ∥E [yk − y ∗ ]∥} ≤ T −1 S −1 op (∥sxk ∥ + ∥syk ∥) for some upper triangular matrix T defined in (35) and a lower triangular matrix S defined in (39), and where vectors sxk , syk are some transformed bias vectors defined in (43), (44) respectively. Furthermore, p p ˆ 2 2 ∥sxk ∥ + ∥syk ∥ ≤ e−kcα κ(L) ∥sx0 ∥ + e−kdβ κ(G) ∥sy0 ∥ + ae−αĉ(k−1) ∥x̂0 ∥ + be−β d(k−1) ∥ŷ0 ∥ | {z } Decaying initial conditions

∞ + F1 (α + β) + F2 MSE∞ x (α, β) + MSEy (α, β)

Asymptotic bias

+ Rα,β , }

(7)

symmetric positive definite matrices solutions of the Lyapunov equations (45), (49), c = (8 ∥L∥op )−1 and d = (8 ∥G∥op )−1 , constants a, b, ĉ, dˆ are related to decaying initial errors from the MSE and can be found in (53), and the MSE terms are defined in Theorem 1. We detail the proof of the theorem in Appendix D, and present its main arguments in the following proof overview. Proof Overview. Recall that Biasxk ≜ E [xk − h(y ∗ )] and biasyk ≜ E [yk − y ∗ ] are the bias vectors. As used in [29, 10], we define the concatenation of the two bias vectors as bk ≜ (Biasxk , biasyk )T . This manipulation allows us to analyze the coupled dynamics and to control cross impacts of both iterates on each other. Using assumption B.1 we use a Taylor expansion of the dynamics f¯ and ḡ around (h(y ∗ ), y ∗ ) to obtain bk+1 = Hbk + Rk , where H represents the linear dynamics around the equilibrium and the remainder term Rk contains the higher order terms of the Taylor expansion plus the expected impact of the Markovian noise. As we will see later, the remainder term Rk can be controlled somehow similarly to what we did for the MSE. The intricate part is to find a way to analyze the dynamics induced by H where the off diagonal term represents the cross influence between the biases of both iterates. To treat this term, we use a linear transformations of coordinates to decouple the fast and slow part, following techniques from control theory, see e.g. [22, Chapter 2]. First, we construct an upper triangular matrix T with identity on the diagonal and an unknown upper right block denoted by U , such that Π ≜ T HT −1 is lower triangular with a zero in the upper right block. This leads to solving a Riccati equation for U [22, Eq. 2.13]. Second, we search for a transform of the matrix Π that would make it block-diagonal, by removing the lower left block. To do that, we build a lower triangular matrix S that has an unknown block X, and consider Λ = SΠS −1 with the goal to make it diagonal matrix. We do that by solving a Sylvester equation in X [32, Lemma 8.3]. Using both transformations, we get sk+1 = Λsk + ST Rk , where sk ≜ ST bk . By having the matrix Λ being block diagonal, with two blocks Λ1 , Λ2 , we can unroll the recursion to give us access to individual updates of each component: sxk+1 = Λk1 sx0 +

k X

Λk−j Rjx , 1

syk+1 = Λk2 sy0 +

j=0

k X

Λk−j Rjy . 2

j=0

The last challenging part is to find the right norm so to ensure that the leading matrices behave as contractions for small enough step sizes. This is done by solving a set of Lyapunov equations to get symmetric positive definite matrices that contain information about the linear dynamics of the system and such that the leading matrices are contractions under the Lyapunov induced norms. We can then propagate back the bounds to the original bias vectors by applying the inverse operation bk = (ST )−1 sk . 5.2

Discussion

The asymptotic bias provided in Theorem 3 is of the order of O(α + β 2 /α2 ). This is not a direct consequence of the MSE bound of Theorem 1 which, by Jensen’s inequality, would only provide a p 2 O( α + β /α2 ) on the bias. To prove Theorem 3, we follow a different path from the one we used to prove the MSE: the approach we used was to treat the fast and slow iterates separately. We first obtained a bound on the MSE of x and then plugged this term in the analysis of y. For the bias, we believe that this approach does not work and that the two iterates must be treated at the same time. Hence, what we do in the proof of Theorem 3 is to couple the analysis of the fast and slow iterates by using a change of variables linked to linearization of f¯ and ḡ around (h(y ∗ ), y ∗ ). The first term in the bias bound (7) is an exponentially fast vanishing term. This term is similar to the one that we obtained for the MSE but is there because of a different reason. In the MSE analysis, this term came from the strong monotonicity assumption. Here, this term is linked to the Hurwitz stability of first order dynamics of the TTSA (Assumption B.1). The proof overview above explains how a 8

0.08 0.06

Fast iterate xk Slow iterate yk y* 0

2000 4000 6000 8000 10000 Iterate k

(a) Transient regime of (102)

1.0

Curve Curve = 1.5 MSE x (simulation) MSE y (simulation)

0.04 0.02 0.00

Curve ( / )2 = 0.01 ( = 0.95) = 0.5 ( = 0.95)

0.8 MSEx

0.10

0.7 0.6 0.5 0.4 0.3 0.2

0.6 0.4 0.2

0.02

0.04

0.06

0.08

(b) MSE·∞ of (102)

0.10

0.0

0.2

0.4 0.6 Ratio /

0.8

1.0

Figure 1: Illustration of the convergence results.

transformation of our bias vectors into new coordinates allowed us to get a contraction. Obtaining this contraction allows us to derive the two remaining terms that form the asymptotic bound on the bias. These terms emerge from two different sources: the first source is the Markovian noise and induces the term F1 (α + β). The other terms MSEx∞ (α, β) + MSEy∞ (α, β) are inherited errors from the MSE due to the high order terms from the Taylor approximation. Taken together, the bias of both iterates has a leading dependence of order O(α + β), which matches the results of linear TTSA under Markovian noise [26]. In addition to that, we get an extra term β2 α2 due to non linearity of the dynamics. This term is injected from the MSE results and can not be improved unless we improve the bound on the MSE. When assuming a clear separation of timescales such that the slow step-size β satisfies β ≤ α3/2 the bound on the bias reduces to O(α) for both iterates. Similarly to the MSE bound, we can show that this bound is tight: We prove in Theorem 4 that there exists a TTSA that has a lower bound on the bias of order Ω(α). The exact statement of the theorem along with its proof can be found in Appendix C. Note that, contrary to the MSE case where linear and non-linear TTSA behave very differently, the bias in non-linear TTSA is essentially of the same order as the bias in the linear setting of [26] where the order O(α + β) is proved to be tight.

Illustrations

In this section, we illustrate some implications of the main theorems of this paper. More numerical results are provided in Appendix F, which also contains the detailed definitions of the examples. In a first example, we consider the TTSA defined in Eq. (102) of Appendix F. In Figure 1(a), we plot a trajectory of the iterates (xk , yk ) as a function of k for α = 0.005 and β = α1.5 . In this example, h(y) = y and y ∗ ≈ 0.57. As observed on this figure, the fast iterate xk tracks h(yk ) = yk but with a relatively large error. When k goes to infinity, both iterates get close to y ∗ ≈ 0.57. This illustrates the vanishing term present in Theorems 1 and 3. In Figure 1(b), we show the MSE of x and y as a function of α, for the case β = α1.5 . We observe that the MSE of x is close to α and that the MSE of y is close to β = α1.5 . Hence, the MSE of y is smaller than what is suggested by Theorem 1. We believe that this is due to the fact that the dynamics of f¯ and ḡ are differentiable in y ∗ (and not just Lipschitz) and that for this particular case stronger guarantees for the MSE of y could be obtained. In our second example, we consider the TTSA where the fast variable xk tries to track yk by applying −xk ) xk+1 = xk + α(yk − xk ), while yk tries to escape from xk by applying yk+1 = yk + β |y(ykk−x γ for k| some γ < 1. As shown in Figure 1(c), for γ = 0.95, the MSE of x in this example is essentially equal to β 2 /α2 , independently of α (remark that the figure contains two curves for each of the values α = 0.01 and α = 0.5, but they are superposed). This shows that the dependence on β 2 /α2 cannot be avoided. To be more precise, this example does not satisfy A.4 (because g is not Lipschitz nor strong monotonic), but our method to obtain the bound on MSEx (in (4)) does not use A.4. This means that the term β 2 /α2 that is present in (4) cannot be improved without using Asssumption A.4. The bounds of all of our Theorems contain terms in β 2 /α2 that all come from MSEx . As we were unable to find an example that has this β 2 /α2 term and that satisfies A.4, it is possible that the dependence 9

on β 2 /α2 could be removed under A.4 but this would require a radically different analysis. We refer to Appendix F for a detailed discussion.

Conclusion

In this paper, we carry out a finite time analysis on the bias and MSE of nonlinear two time-scale stochastic approximation (TTSA) under constant step-sizes. We consider a general model with Markovian noise where the transition kernel can depend on the iterates (xk , yk ), this allows for broad applications, for instance in reinforcement learning. We show that the bias and MSE are upper bounded by vanishing initialization terms plus a residual term that scales as O(α + β 2 /α2 ) and we show that our analysis is tight. While our current framework establishes robust theoretical limits for TTSA, there are several directions for future work. The first is to improve the bounds by (slightly) strenghening the assumptions: for instance, could the MSE bounds be improved if we add B.1 to Theorem 1? Also, one may want to relax the finite state-space restriction on the Markovian noise. This will require careful technical adjustments, particularly in validating the Poisson equation solutions required to bound the noise impact in continuous or unbounded spaces. Furthermore, replacing the strict compactness assumption on the iterate trajectory with non-uniform bounds would capture a larger class of practical, unconstrained algorithms.

References [1] Sebastian Allmeier and Nicolas Gast. Computing the bias of constant-step stochastic approximation with markovian noise. Advances in Neural Information Processing Systems, 37:137873–137902, 2024. [2] Francis Bach and Eric Moulines. Non-strongly-convex smooth stochastic approximation with convergence rate o(1/n). In C.J. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger, editors, Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013. [3] Albert Benveniste, Michel Métivier, and Pierre Priouret. Adaptive algorithms and stochastic approximations. In Applied Mathematics, 1990. [4] Vivek S Borkar. Stochastic approximation: a dynamical systems viewpoint, volume 100. Springer, 2008. [5] Vivek S. Borkar and Sean P. Meyn. The o.d.e. method for convergence of stochastic approximation and reinforcement learning. SIAM J. Control. Optim., 38:447–469, 2000. [6] Siddharth Chandak, Shaan Ul Haque, and Nicholas Bambos. Finite-time bounds for two-timescale stochastic approximation with arbitrary norm contractions and markovian noise. In 2025 IEEE 64th Conference on Decision and Control (CDC), pages 6095–6101, 2025. [7] Zaiwei Chen, Sheng Zhang, Thinh T Doan, John-Paul Clarke, and Siva Theja Maguluri. Finitesample analysis of nonlinear stochastic approximation with applications in reinforcement learning. Automatica, 146:110623, 2022. [8] Zixi Chen, Yumin Xu, and Ruixun Zhang. Convergence rate in a nonlinear two-time-scale stochastic approximation with state (time)-dependence. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, pages 15993–16000, 2025. [9] Thinh T Doan. Finite-time analysis and restarting scheme for linear two-time-scale stochastic approximation. SIAM Journal on Control and Optimization, 59(4):2798–2819, 2021. [10] Thinh T. Doan. Finite-time convergence rates of nonlinear two-time-scale stochastic approximation under markovian noise, 2021. [11] Thinh T. Doan. Nonlinear two-time-scale stochastic approximation: Convergence and finite-time performance. IEEE Transactions on Automatic Control, 68(8):4695–4705, 2023. [12] Randal Douc, Eric Moulines, Pierre Priouret, and Philippe Soulier. Markov chains. Operation research and financial engineering. Springer, 2018. 10

[13] John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12(61):2121–2159, 2011. [14] Alain Durmus, Eric Moulines, Alexey Naumov, and Sergey Samsonov. Finite-time highprobability bounds for polyak–ruppert averaged iterates of linear stochastic approximation. Mathematics of Operations Research, 50(2):935–964, 2025. [15] Alain Durmus, Eric Moulines, Alexey Naumov, Sergey Samsonov, Kevin Scaman, and Hoi-To Wai. Tight high probability bounds for linear stochastic approximation with fixed stepsize. Advances in Neural Information Processing Systems, 34:30063–30074, 2021. [16] Ian J Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. Advances in neural information processing systems, 27, 2014. [17] Shaan Ul Haque, Sajad Khodadadian, and Siva Theja Maguluri. Tight finite time bounds of two-time-scale linear stochastic approximation with markovian noise, 2025. [18] Jie Hu, Vishwaraj Doshi, et al. Central limit theorem for two-timescale stochastic approximation with markovian noise: Theory and applications. In International Conference on Artificial Intelligence and Statistics, pages 1477–1485. PMLR, 2024. [19] Dongyan Huo, Yudong Chen, and Qiaomin Xie. Bias and extrapolation in markovian linear stochastic approximation with constant stepsizes. In Abstract Proceedings of the 2023 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems, pages 81–82, 2023. [20] Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In International conference on machine learning, pages 4882–4892. PMLR, 2021. [21] Maxim Kaledin, Eric Moulines, Alexey Naumov, Vladislav Tadic, and Hoi-To Wai. Finite time analysis of linear two-timescale stochastic approximation with markovian noise. In Conference on Learning Theory, pages 2144–2203. PMLR, 2020. [22] Petar Kokotović, Hassan K Khalil, and John O’reilly. Singular perturbation methods in control: analysis and design. SIAM, 1999. [23] Vijay Konda and John Tsitsiklis. Actor-critic algorithms. In S. Solla, T. Leen, and K. Müller, editors, Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999. [24] Vijay R. Konda and John N. Tsitsiklis. Convergence rate of linear two-time-scale stochastic approximation. The Annals of Applied Probability, 14(2), May 2004. [25] Harold J. Kushner and George Yin. Stochastic approximation algorithms and applications. In Applied Mathematics, 1997. [26] Jeongyeol Kwon, Luke Dotson, Yudong Chen, and Qiaomin Xie. Two-timescale linear stochastic approximation: Constant stepsizes go a long way. In The 28th International Conference on Artificial Intelligence and Statistics, 2025. [27] Chandrashekar Lakshminarayanan and Csaba Szepesvari. Linear stochastic approximation: How far does constant step-size and iterate averaging go? In Amos Storkey and Fernando Perez-Cruz, editors, Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1347–1355. PMLR, 09–11 Apr 2018. [28] Hamid Maei, Csaba Szepesvári, Shalabh Bhatnagar, Doina Precup, David Silver, and Richard S Sutton. Convergent temporal-difference learning with arbitrary smooth function approximation. In Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems, volume 22. Curran Associates, Inc., 2009. [29] Abdelkader Mokkadem and Mariane Pelletier. Convergence rate and averaging of nonlinear two-time-scale stochastic approximation algorithms. The Annals of Applied Probability, 16(3), August 2006. 11

[30] Deanna Needell, Nathan Srebro, and Rachel Ward. Stochastic gradient descent, weighted sampling, and the randomized kaczmarz algorithm. Advances in neural information processing systems, 27, 2014. [31] B. T. Polyak and A. B. Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4):838–855, 1992. [32] Alexander S. Poznyak. Chapter 8 - linear matrix equations. In Advanced Mathematical Tools for Automatic Control Engineers: Deterministic Techniques, pages 133–137. Elsevier, Oxford, 2008. [33] Herbert Robbins and Sutton Monro. A stochastic approximation method. Annals of mathematical statistics, 22:400–407, 1951. [34] Richard S Sutton, Hamid Maei, and Csaba Szepesvári. A convergent o(n) temporal-difference algorithm for off-policy learning with linear function approximation. In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008.

Main notations

In this section, we recall some notation and main variables that are used in all appendix. Pd Spaces and norms Let x ∈ Rd . The quantity ∥x∥ is the Euclidean norm ∥x∥ = ( i=1 x2i )1/2 . For y ∈ Rd we denote the inner product between x and y by ⟨x, y⟩. For a square matrix M ∈ Rn×n we define the operator norm ∥M ∥op = sup∥x∥=1 ∥M x∥. For a positive definite matrix L we √ define the norm induced by L as ∥x∥L = xT Lx. We define the operator norm induced by L by ∥X∥L = sup∥u∥L =1 ∥Xu∥L . The quantity σmax (L) denotes the maximal singular value of L and (L) λmax (L) is its maximal eigenvalue, and the condition number of L is κ(L) = λλmax . min (L)

Probabilities For a discrete random variable Pξ in a countable space Ξ, we have for a probability measure π on Ξ the notation Eξ∼π [f (x, ξ)] = ξ∈Ξ f (x, ξ)π(ξ) and for a kernel P : Ξ × Ξ → [0, 1] P we have for ϵ ∈ Ξ, Eξ∼P (.|ϵ) [f (x, ξ)] = ξ∈Ξ f (x, ξ)P (ξ|ϵ). For an event E the quantity 1{E} equals 1 if E is true and 0 otherwise. Derivatives For a multivariate function f¯ : Rdx × Rdy → Rdx , ∇x f (x, y) denotes the partial i (x,y) Jacobian with respect to the first argument (∇x f (x, y))i,j = ∂f∂x and ∇y f (x, y) denotes the j partial Jacobian with respect to the second argument.

Important variables.

In what follows, we will use the notation:

x̂k ≜ xk − h(yk )

h i 2 Tracking error of xk with respect to h(yk ). Note that MSExk = E ∥x̂k ∥

x̃k ≜ xk − h(y ∗ )

Tracking error of xk with respect to h(y ∗ ). Note that Biasxk = E [x̃k ] i h 2 Tracking error of yk . Note that MSEyk = E ∥ŷk ∥ and Biasyk = E [ŷk ]

ŷk ≜ yk − y ∗

Step sizes To obtain our result on the bias (Theorem 3), we require a few conditions on the step sizes from the contraction under the Lyapunov norms passing by the existence and boundedness of the solutions of the Riccati and Sylvester equations to the invertibility conditions on the bias dynamics. We compile all the conditions here so we can refer to them easily in the statement of our √1 Theorem 3. In what follows, we define Rric = where MJ is a positive constant 4∥G∥op

κ(G)MJ

satisfying ∥Jyx ∥op ≤ MJ . We start with the conditions on α∗ : 1 ∗ α ≜ min , µf 1 2

4 ∥L∥op (∥Jxx ∥L + ε∗4 2 κ(L)(Rric + ∥U0 ∥op )2 ∥Jyx ∥L ) 1 ∥Jxx ∥op + (Rric + ∥U0 ∥op ) ∥Jyx ∥op

(8)

where the first upper bound comes from the MSE Theorem, the second bound gives us the contraction of the fast bias dynamics and can be found in (81) where ε∗4 = min{(74), (75)} and the last condition is for the invertibility of the fast bias dynamics and can be found in (90). Second, we enumerate the conditions on β ∗ : 1 β ∗ ≜ min , µg 1 2

2 ) 4 ∥G∥op (2κ(G) ∥∆∥op MJ Rric + ∥∆∥G + κ(G)MJ2 Rric 1 ∥∆∥op + Rric ∥Jyx ∥op

(9)

where the first condition is required for the results of the MSE theorem to hold, the second condition is reqruired by the contraction of the slow bias dynamics and can be found in (84), the last condition is about the invertibility of the slow dynamics and can be found in (97). Finally, we define ω ∗ to be: ω ∗ ≜ min 1, Rric , −1 Jxx op (Rric + ∥U0 ∥op )(∥Jyy ∥op + (Rric + ∥U0 ∥op ) ∥Jyx ∥op ) 1

−1 Jxx op

, ∥Jyy ∥op + 2 ∥Jyx ∥op (Rric + ∥U0 ∥op ) 1

∥Jyy ∥op

−1 Jxx + 2(Rric + ∥U0 ∥op ) ∥Jyx ∥op op

−1 Jxx op

(10) ,

1 ∗ 2 ∥L∥op (Rric + ∥U0 ∥op )(2α̃ κ(L) ∥Jyx ∥op ∥Jxx ∥op + 2κ(L)(1/2) MJ )

where the first condition is to make the results of the MSE hold, the second and third conditions come from the existence and boundedness of the Riccati equation and can be found in (74), (75), the fourth condition, alongside the third, is for the existence and boundedness of the solution of the Sylvester equation and can be found in (78), the fifth condition is needed for the contraction of the fast bias dynamics and can be found in (82) with noting that α̃∗ in the denominator is independent of α and β and can be found in (81).

Proof of Theorem 1: MSE of TTSA

In this section we prove the upper bound of Theorem 1 on the MSE of TTSA. We restate the theorem, before detailing the proof. Theorem 1. Assuming that A.1-3 hold, and that αµf < 1 and β ≤ α, then for any k ≥ 1 we obtain the following upper bound on the MSE of the fast iterate, B1 β2 2 2 k E[∥xk − h(yk )∥ ] ≤ (1 − αµf ) ∥x0 − h(y0 )∥ + α + B2 + β B3 + 2 B4 , (4) 1 − αµf α | {z } | {z } Decaying initial condition MSEx ∞ (α,β)

E[∥yk − y ∗ ∥2 ] ≤ (1 − βµg )k ∥y0 − y ∗ ∥ + D1 β k(1 − (αµf ∧ βµg ))k ∥x0 − h(y0 )∥ {z } | Decaying initial conditions

D3 + D2 MSEx∞ (α, β) + β + D4 + αD5 . 1 − βµg | {z } MSEy ∞ (α,β)

where the constants (Dk )k∈[5] are also formally defined in appendix (Eq. (28)-(32)) and depend again on problem parameters but not on α, β, k. Proof. We divide the proof into two main parts, the first part is dedicated to the fast iterate MSE and the second part is for the slow iterate. Upper bounding the MSE of the fast iterate the shorthand

We recall that x̂k ≜ xk − h(yk ). We also introduce

ψk ≜ f (xk , yk , ξk+1 ) − f (xk , yk ), 14

∆hk ≜ h(yk+1 ) − h(yk ) .

Using the update xk+1 = xk − αf (xk , yk , ξk+1 ), we obtain x̂k+1 = xk+1 − h(yk+1 ) = x̂k − αf (xk , yk , ξk+1 ) − ∆hk , so expanding the squared norm gives 2

∥x̂k+1 ∥ = ∥xk+1 − h(yk+1 )∥

= x̂k − αf (xk , yk , ξk+1 ) − ∆hk 2

= ∥x̂k ∥ + α2 ∥f (xk , yk , ξk+1 )∥ − 2α x̂k , f (xk , yk ) + ∥∆hk ∥2 − 2 x̂k , ∆hk + 2α f (xk , yk , ξk+1 ), ∆hk − 2α ⟨x̂k , ψk ⟩ Rearranging the terms to facilitate the organization of the bound: 2

∥x̂k+1 ∥ = ∥x̂k ∥ − 2α x̂k , f (xk , yk ) − 2 x̂k , ∆hk

(11)

+ α ∥f (xk , yk , ξk+1 )∥

(12)

+ ∥∆hk ∥2 + 2α

(13)

f (xk , yk , ξk+1 ), ∆hk

− 2α ⟨x̂k , ψk ⟩ .

(14)

We first control the bounded second-order terms. By Lemma 1, there exist constants Mf , Mh > 0 such that ∥f (xk , yk , ξk+1 )∥ ≤ Mf , ∆hk ≤ βMh . Therefore, by the Cauchy–Schwarz inequality, 2

(12) + (13) ≤ α2 ∥f (xk , yk , ξk+1 )∥ + ∆hk

+ 2α ∥f (xk , yk , ξk+1 )∥ ∆hk

≤ α2 M2f + β 2 M2h + 2αβ Mf Mh . We now extract the contraction term. By property of h, for every y ∈ Y, f (h(y), y) = 0. Hence, recalling that x̂k = xk − h(yk ), we can write −2α x̂k , f (xk , yk ) = −2α xk − h(yk ), f (xk , yk ) − f (h(yk ), yk ) ≤ −2αµf ∥xk − h(yk )∥

= −2αµf ∥x̂k ∥ , where the inequality follows from the strong monotonicity assumption A.3 applied to the first argument of f , with y = yk . The third term in line (11) is controlled by Cauchy–Schwarz and Young’s inequality. Precisely: 2

−2 x̂k , ∆hk ≤ 2 ∥x̂k ∥ ∆hk ≤ αµf ∥x̂k ∥ +

1 β 2 M2h 2 2 ∆hk ≤ αµf ∥x̂k ∥ + , αµf αµf

where the last inequality follows from Lemma 1. Combining the previous bounds on (11)–(13), we obtain 2

∥x̂k+1 ∥ ≤ (1 − αµf ) ∥x̂k ∥ − 2α ⟨x̂k , ψk ⟩ + α2 M2f + β 2 M2h + 2αβMf Mh +

β 2 M2h . αµf

For brevity, set Mα,β ≜ α2 M2f + β 2 M2h + 2αβMf Mh +

β 2 M2h . αµf

Iterating the previous inequality gives 2

∥x̂k+1 ∥ ≤ (1 − αµf )k+1 ∥x̂0 ∥ − 2α

k k X X (1 − αµf )k−j ⟨x̂j , ψj ⟩ + (1 − αµf )k−j Mα,β . j=0

j=0

Taking expectations yields h i 2 2 E ∥x̂k+1 ∥ ≤(1 − αµf )k+1 ∥x̂0 ∥ + Mα,β

(15)

k X (1 − αµf )j

(16)

j=0

  k X − 2αE  (1 − αµf )k−j ⟨x̂j , ψj ⟩ .

(17)

j=0

The first term (15) is the contribution of the initial condition and decays geometrically. The two remaining terms correspond to the steady-state error induced by the constant step-size. We first bound the deterministic geometric term (16): Mα,β

k X

(1 − αµf )j = Mα,β

j=0

1 − (1 − αµf )k+1 Mα,β ≤ . αµf αµf

It remains to control the Markovian-noise term (17). To do so, we invoke Lemma 2, that we detail and prove in Appendix E.   k X Lv Mf (MS + Mh ) 3 − αµf Mf (1 + |Ξ|MS LK ) αE  (1 − αµf )k−j ⟨x̂j , ψj ⟩ ≤α Mξ + + f 1 − αµf µf µf j=0 Mξ Mh + Lvf Mg (MS + Mh ) + |Ξ|MS LK Mξ Mg +β µf Combining the above results provide h i 2 2 E ∥x̂k+1 ∥ ≤(1 − αµf )k+1 ∥x̂0 ∥ ! M2f + 2(Mξ (1 + |Ξ|MS LK ) + Lvf (MS + Mh ))Mf 4Mξ +α + 2Mξ + 1 − αµf µf 2(Mξ + Mf )Mh + 2Lvf Mg (MS + Mh ) + |Ξ|MS LK Mξ Mg β 2 M2 β 2 M2h +β + 2 2h . + µf α µf α µf 2

Using that βα ≤ β because β ≤ α, we are ready to formally define the constants B1 := 4Mξ , B2 := 2Mξ +

(18) M2f + 2(Mξ (1 + |Ξ|MS LK ) + Lvf (MS + Mh ))Mf µf

(19)

B3 :=

2(Mξ + Mf )Mh + 2Lvf Mg (MS + Mh ) + 2|Ξ|MS LK Mξ Mg + M2h , and µf

(20)

B4 :=

M2h . µ2f

(21)

This proves the desired MSE bound for the fast iterate, and allows us to define the non-vanishing term of the bound, B1 β2 + B2 + β B3 + 2 B4 . MSEx∞ (α, β) := α 1 − αµf α MSE of the slow iterate We now turn to the MSE of the slow iterate. The argument follows the same structure as for the fast iterate, with one additional difficulty: the update of yk is evaluated at (xk , yk ), whereas the contraction assumption is available for the reduced dynamics evaluated at (h(yk ), yk ). We therefore need to control the tracking error of the fast iterate, which enters the slow recursion as an additional residual term. 16

For convenience, introduce ϕk ≜ g(xk , yk , ξk+1 ) − g(xk , yk ). ∆gk ≜ g(xk , yk ) − g(h(yk ), yk ), g Thus, ∆k measures the error due to replacing xk by h(yk ) in the noisy averaged dynamics, while ϕk is the Markovian-noise fluctuations around the averaged slow dynamics. Using the update yk+1 = yk − βg(xk , yk , ξk+1 ), we have ŷk+1 = ŷk − βg(xk , yk , ξk+1 ). Expanding the squared norm and adding/subtracting the reduced dynamics gives 2

∥ŷk+1 ∥ = ∥ŷk − βg(xk , yk , ξk+1 )∥

= ∥ŷk ∥ − 2β ⟨ŷk , g(xk , yk , ξk+1 )⟩ + β 2 ∥g(xk , yk , ξk+1 )∥

= ∥ŷk ∥ − 2β ⟨ŷk , g(h(yk ), yk )⟩ − 2β ⟨ŷk , g(xk , yk ) − g(h(yk ), yk )⟩ 2

+ β 2 ∥g(xk , yk , ξk+1 )∥ − 2β ⟨ŷk , g(xk , yk , ξk+1 ) − g(xk , yk )⟩ 2

= ∥ŷk ∥ − 2β ⟨ŷk , g(h(yk ), yk )⟩ − 2β ⟨ŷk , ∆gk ⟩ 2

+ β ∥g(xk , yk , ξk+1 )∥ − 2β ⟨ŷk , ϕk ⟩ .

(22) (23) (24)

By the strong monotonicity assumption A.4, the first two terms in (22) yield a contraction. The additional term involving ∆gk measures the error induced by evaluating the slow update at xk instead of h(yk ). More precisely, 2

(22) = ∥ŷk ∥ − 2β ⟨ŷk , g(h(yk ), yk )⟩ − 2β ⟨ŷk , ∆gk ⟩ 2

≤ ∥ŷk ∥ − 2βµg ∥ŷk ∥ + 2β ∥ŷk ∥ ∥∆gk ∥ 2

≤ ∥ŷk ∥ − 2βµg ∥ŷk ∥ + βµg ∥ŷk ∥ + 2

≤ (1 − βµg ) ∥ŷk ∥ +

β 2 ∥∆gk ∥ µg

≤ (1 − βµg ) ∥ŷk ∥ + β

L2g 2 ∥x̂k ∥ . µg

In the first inequality, we used the strong monotonicity of the reduced slow dynamics and Cauchy– Schwarz. In the second inequality, we used Young’s inequality 2ab ≤ ϵa2 + b2 /ϵ with ϵ = µg . The last inequality follows from the Lipschitz property of g(·, ·): ∥∆gk ∥ = ∥g(xk , yk ) − g(h(yk ), yk )∥ ≤ Lg ∥xk − h(yk )∥ = Lg ∥x̂k ∥ We take β < µ−1 g so that the factor 1 − βµg is strictly contractive and bounded away from zero. The remaining deterministic term in (23) is bounded using Lemma 1: there exists Mg > 0 such that 2

β 2 ∥g(xk , yk , ξk+1 )∥ ≤ β 2 M2g . Combining the previous bounds on (22)–(24), we obtain 2

∥ŷk+1 ∥ ≤ (1 − βµg ) ∥ŷk ∥ +

L2g 2 β ∥x̂k ∥ + β 2 M2g − 2β ⟨ŷk , ϕk ⟩ . µg

Iterating this recursion from time 0 to time k, and then taking expectations, gives k h i X 2 2 E ∥ŷk+1 ∥ ≤(1 − βµg )k+1 ∥ŷ0 ∥ + M2g β 2 (1 − βµg )j

(25)

j=0 k h i L2g X 2 β (1 − βµg )k−j E ∥x̂j ∥ µg j=0   k X − 2βE  (1 − βµg )k−j ⟨ŷj , ϕj ⟩ .

j=0

(26)

(27)

The first term in (25) is the contribution of the initial condition and decays geometrically. The second term in (25) is a geometric sum, and therefore M2g β 2

k X 1 − (1 − βµg )k+1 (1 − βµg )j = M2g β 2 βµg j=0

≤

M2g β. µg

We now control the term (26), which is the geometrically weighted accumulation of the fast-iterate MSEs. Using the upper bound we obtained for the fast iterate, we get k k h i L2 X L2g X g 2 2 β (1 − βµg )k−j E ∥x̂j ∥ ≤ β (1 − βµg )k−j (1 − αµf )j+1 ∥x̂0 ∥ + MSEx∞ (α, β) µg j=0 µg j=0

≤

L2g MSEx∞ (α, β) µ2g k X L2g 2 + (1 − βµg )k−j (1 − αµf )j . β(1 − αµf ) ∥x̂0 ∥ µg j=0 | {z } Sk

For the term Sk , we simply remark that, for j ∈ [k], it holds that (1 − βµg )k−j (1 − αµf )j ≤ (1 − (αµf ∧ βµg ))k . Therefore, we can upper bound each term of the sum uniformly and obtain that Sk ≤ (k + 1)(1 − (αµf ∧ βµg ))k . Note that the bound is tight for αµf = βµg . This leads to k h i L2 L2g L2g X g 2 2 β (1 − βµg )k−j E ∥x̂j ∥ ≤ 2 MSEx∞ (α, β) + β ∥x̂0 ∥ (k + 1)(1 − (αµf ∧ βµg ))k+1 . µg j=0 µg µg

It remains to control the Markovian-noise term (27). As in the fast-iterate analysis, we use the Poisson-equation decomposition from Lemma 2. In particular,     k k X X −2βE  (1 − βµg )k−j ⟨ŷj , ϕj ⟩ ≤2β E  (1 − βµg )k−j ⟨ŷj , ϕj ⟩ j=0

j=0

(Lvg + Mξ LK |Ξ|)MS Mg ≤2β + Mξ µg (Lvg + Mξ LK |Ξ|)MS Mf + 2α µg

3 − βµg Mg + 1 − βµg µg

Combining this bound with the bounds on (25) and (26), we obtain h i M2g 2 2 E ∥ŷk+1 ∥ ≤(1 − βµg )k+1 ∥ŷ0 ∥ + β µg L2g 2 ∥x̂0 ∥ (k + 1)(1 − (αµf ∧ βµg ))k+1 µg L2g (Lvg + Mξ LK |Ξ|)MS Mg + Mξ Mg 2Mξ x + 2 MSE∞ (α, β) + 2β + + Mξ µg µg 1 − βµg

+β

+ 2α

(Lvg + Mξ LK |Ξ|)MS Mf . µg 18

It remains to identify the constants, D1 :=

L2g , µg

(28)

D2 :=

L2g , µ2g

(29)

D3 := 4Mξ , and (30) 2 Mg (Lvg + Mξ LK |Ξ|)MS Mg + Mξ Mg D4 := +2 + Mξ . (31) µg µg (Lvg + Mξ LK |Ξ|)MS Mf D5 := 2 . (32) µg This proves the desired MSE bound for the slow iterate and concludes the proof of the theorem.

Proof of Theorem 2: Lower bounds

Theorem 2 shows that there exists a TTSA whose MSE is larger than Ω(α) for both the fast and the slow iterates. This result is a consequence of Theorem 4 below where we exhibit an example that attains this lower bound, and where we also provide an example has an MSE of order Ω(α) for both x and y. Theorem 4. There exist a constant C and a constant α0 > 0 and: (i) a TTSA satisfying all Assumptions of Theorem 1 such that for all β ≤ α ≤ α0 : lim inf E ∥xk − h(yk )∥2 ≥ Cα and lim inf E ∥yk − y ∗ ∥2 ≥ Cα; k→∞

k→∞

(ii) a TTSA satisfying all Assumptions of Theorem 3 such that for all β ≤ α ≤ α0 : lim inf ∥E [xk ] − h(y ∗ )∥ ≥ C(α + β)

and

k→∞

lim inf ∥E [yk ] − y ∗ ∥ ≥ C(α + β). k→∞

Proof. The proof of (i) and (ii) uses the same fast iterate (xk ) that is defined as: xk+1 = xk − α(f (xk ) + ξk+1 ), x where f (x) = 1+x . The sequence of noises (ξk ) is a sequence of zero-mean i.i.d. scaled Rademacher random variables (i.e., such that P[ξk = 1/2] = P[ξk = −1/2] = 0.5). Note that the function f is Lipschitz continuous with a Lipschitz-continunous derivative. We assume that α ∈ (0, 1) and we set x0 = 0.

We first show by induction on k that the sequence xk is bounded almost surely, and more precisely that xk ∈ [−1/2, 1]. Assume that this holds for some k and recall that xk + ξk+1 . xk+1 = xk − α 1 + xk We distinguish two cases: xk 1. If 1+x + ξk ≥ 0, then k

xk+1 ≤ xk ≤ 1 xk x2k 1 1 xk+1 ≥ xk − + ξk = − ≥− , 1 + xk 1 + xk 2 2 where we used that α ≤ 1 in the second line and the fact that x2 /(1 + x) ≥ 0. xk 2. If 1+x + ξk ≤ 0, then k

1 xk+1 ≥ xk ≥ − 2 xk x2k 1 xk+1 ≤ xk − + ξk = + ≤ 1, 1 + xk 1 + xk 2 19

As the behavior of xk does not depend on the slow iterates (yk )k , the sequence (xk ) can be analyzed with standard tools from one-time-scale stochastic approximation. Here, we have x∗ = h(yk ) = 0. In particular, since near 0 this function has a second order approximation of x − x2 ,from [1], we know that the bias and the MSE of xk are of order α and asymptotically equal to: α lim E [xk ] = lim E (xk )2 = + O(α2 ). k→∞ k→∞ 2c for some positive constant c. Moreover, from [25, Chapter 10], when α is small and k is large, the variable xk is close to a Gaussian distribution of variance α/2c, which implies p (33) lim E [|xk |] = α/dπ + O(α) k→∞

for some positive constant d. For the variable yk , we consider two different cases to prove either (i) or (ii): (i) To prove (i), we set yk+1 = yk − β(yk − |xk |). Unrolling the recurrence equation yk+1 = yk (1 − β) + β|xk | shows that yk+1 = β

k X (1 − β)k−i |xi |. i=0

Combined with (33), this shows that lim E [yk ] = β

k→∞

∞ X

(1 − β)i lim E [|xk |] k→∞

i=0

= lim E [|xk |] k→∞ r α = + O(α). dπ By using Jensen’s inequality, this shows that: lim E (yk )2 ≥ lim (E [yk ])2 k→∞

k→∞

α + O(α3/2 ). dπ

(ii) To prove (ii), we define yk as: yk+1 = yk − β(yk − xk ) This shows that limk→∞ E [yk ] = limk→∞ E [xk ] = α/2c + O(α2 ).

The proof of the above theorem exhibits an example where the MSE of y can be of order Ω(α). To construct this example, we used a function g such that ḡ(x, y) = y − |x| that satisfies all assumptions of Theorem 1 but not the assumptions of Theorem 3 because it is not differentiable. In fact, we believe that having a non-differentiable function is needed to obtain the Ω(α) lower bound for MSE∞ y , and that under a stronger assumption of differentiability, one may obtain a tighter bound for the MSE. This is an interesting question for future work.

Proof of Theorem 3 : Bias of TTSA

We start by restating the theorem, before detailing its proof. 20

Theorem 3. Assuming A.1-4 and B.1, there exists α∗ defined in (8), β ∗ defined in (9), and ω ∗ defined β in (10) such that for α < α∗ , β < β ∗ and α < ω ∗ we have the following bound on the bias of both iterates: max{∥E [xk − h(y ∗ )]∥ , ∥E [yk − y ∗ ]∥} ≤ T −1 S −1 op (∥sxk ∥ + ∥syk ∥) for some upper triangular matrix T defined in (35) and a lower triangular matrix S defined in (39), and where vectors sxk , syk are some transformed bias vectors defined in (43), (44) respectively. Furthermore, p p ˆ 2 2 ∥sxk ∥ + ∥syk ∥ ≤ e−kcα κ(L) ∥sx0 ∥ + e−kdβ κ(G) ∥sy0 ∥ + ae−αĉ(k−1) ∥x̂0 ∥ + be−β d(k−1) ∥ŷ0 ∥ | {z } Decaying initial conditions

∞ + F1 (α + β) + F2 MSE∞ x (α, β) + MSEy (α, β)

Asymptotic bias

+ Rα,β , }

(7)

where F1 , F2 are positive constants independent of α, β and k, whose values may depend on the problem parameters and on ω ∗ , Rα,β contains terms of higher order (α + β)2 . Additionally, L, G are symmetric positive definite matrices solutions of the Lyapunov equations (45), (49), c = (8 ∥L∥op )−1 and d = (8 ∥G∥op )−1 , constants a, b, ĉ, dˆ are related to decaying initial errors from the MSE and can be found in (53), and the MSE terms are defined in Theorem 1. Proof. We start by detailing the general outline of the proof. We begin by concatenating the TTSA updates to jointly analyze the bias of the fast and slow variables. Exploiting the differentiability of the underlying dynamics (Assumption B. 1), we apply a first-order Taylor expansion to linearize the updates, which yields a highly coupled linear dynamical system. To mathematically decouple these dynamics, we introduce a two-stage coordinate transformation. First, we block-triangularize the system by solving an algebraic Riccati equation. Then, we fully diagonalize the system by solving a Sylvester equation. This transformation projects the bias into a new coordinate space governed by a diagonalized recursion. While this recursion features strong contracting terms, it is corrupted by higher-order perturbation residuals arising from the diagonalization process and the multiplication by the transformation matrices. However, the order of these perturbations will not affect the nature of the contractions due to our assumptions on the step sizes. The diagonal structure then allows us to unroll the recursion and bound each transformed component independently. Finally, because our transformation matrices possess finite norms, and by a simple triangular inequality we get the final bounds on the bias. Concatenating the iterates. We start by deriving a recursion for the joint bias vector. In this proof, we work with deviations from the equilibrium and write x̃k ≜ xk − h(y ∗ ),

and we recall that ŷk = yk − y ∗ .

We also introduce the centered Markovian-noise terms ψk ≜ f (xk , yk , ξk+1 ) − f (xk , yk ),

φk ≜ g(xk , yk , ξk+1 ) − g(xk , yk ).

For the fast iterate, we obtain E [x̃k+1 ] = E [x̃k ] − αE f (xk , yk ) − αE [ψk ] = E [x̃k ] − αE f h(y ∗ ) + x̃k , y ∗ + ŷk − αE [ψk ] = E [x̃k ] − αE f (h(y ∗ ), y ∗ ) + Jxx x̃k + Jxy ŷk + Rf (x̃k , ŷk ) − αE [ψk ] = E [x̃k ] − αJxx E [x̃k ] − αJxy E [ŷk ] − αE [Rf (x̃k , ŷk )] − αE [ψk ] , where we used f (h(y ∗ ), y ∗ ) = 0 in the last line. Similarly, for the slow iterate, E [ŷk+1 ] = E [ŷk ] − βE [g(xk , yk )] − βE [φk ] = E [ŷk ] − βE g h(y ∗ ) + x̃k , y ∗ + ŷk − βE [φk ] = E [ŷk ] − βE [g(h(y ∗ ), y ∗ ) + Jyx x̃k + Jyy ŷk + Rg (x̃k , ŷk )] − βE [φk ] = E [ŷk ] − βJyx E [x̃k ] − βJyy E [ŷk ] − βE [Rg (x̃k , ŷk )] − βE [φk ] , 21

where we used g(h(y ∗ ), y ∗ ) = 0. In both expansions, we applied a first-order Taylor expansion of the averaged dynamics around the equilibrium (h(y ∗ ), y ∗ ). Since the second derivatives of f¯ and ḡ are bounded by Assumption B.1, there exists constant LR > 0 such that, for all k ≥ 0, 2 2 ∥Rf (x̃k , ŷk )∥ ≤ LR ∥x̃k ∥ + ∥ŷk ∥ , ∥Rg (x̃k , ŷk )∥ ≤ LR ∥x̃k ∥ + ∥ŷk ∥ . Concatenating the two bias vectors yields the joint recursion E [x̃k+1 ] I − αJxx −αJxy E [x̃k ] −αE [Rf (x̃k , ŷk )] − αE [ψk ] = + . (34) E [ŷk+1 ] −βJyx I − βJyy E [ŷk ] −βE [Rg (x̃k , ŷk )] − βE [φk ] {z } {z } | | Rk

Hα,β

The matrix Hα,β is the linearized discrete-time dynamics of the joint bias. Its block structure shows that the two components are still coupled: the fast bias depends on the slow bias through the block Jxy , while the slow bias depends on the fast bias through the block Jyx . This coupling is the main reason why one cannot simply analyze the two bias coordinates separately. The next step is therefore to find a change of coordinates adapted to this two-time-scale linear system. More precisely, following ideas from singular perturbation theory [22, Chapter 2], we construct invertible transformations that turn Hα,β into a block-diagonal matrix. In the transformed coordinates, the fast and slow bias components can then be controlled separately, with their respective contraction rates. In the following, we drop the subscript α, β from the matrix H for conciseness of notation. Triangularization through a Riccati equation. The first step is to remove the upper-right coupling in the linearized dynamics. We look for an invertible change of coordinates T such that Π ≜ T HT −1 is block lower triangular, or equivalently H = T −1 ΠT . We choose T to be upper triangular with identity blocks on the diagonal: I U I −U −1 T = , T = , 0 I 0 I

(35)

where U ∈ Rdx ×dy is to be determined. To simplify the notation, write the block decomposition of H as A B H= . C D Then a direct computation gives T HT

−1

I 0

A + UC C

= =

U I

A C

B D

R(U ) D − CU

I 0

−U I

where R(U ) ≜ −(A + U C)U + B + U D. Thus, the upper-right block is zero if and only if U solves R(U ) = −(A + U C)U + B + U D = 0.

(36)

This is an algebraic Riccati equation for the unknown block U [22, Eq. 2.8]. Any solution of (36) yields the block lower-triangular form A + UC 0 −1 Π = T HT = . C D − CU Writing the Riccati equation in terms of the blocks of the TTSA dynamics gives β β U Jyx U + Jxx U − U Jyy = Jxy . α α 22

(37)

This equation is the singular-perturbation Riccati equation associated with the linearized two-timescale dynamics. Under Assumption B.1, it is shown in [22, Chapter 2, Eq. 2.13] that when β/α small enough, this Riccati equation admits a solution of the form ∞ ℓ X β U= Ũℓ , α ℓ=0

where the matrices (Ũℓ )ℓ≥0 are independent of α and β. The matrix Ũ0 can be found by setting β/α = 0 in (37) which gives Jxx Ũ0 = Jxy . −1 Since −Jxx is Hurwitz by Assumption B.1, the matrix Jxx is invertible and Ũ0 = Jxx Jxy . This shows that ∞ ℓ X β U = Ũ0 + Ũℓ α ℓ=1

−1 = Jxx Jxy +

where Ûα,β ≜

P∞ β ℓ−1 ℓ=1

β Ûα,β , α

Ũℓ .

β From Assumption B.1 and since we are assuming that α < min{ε∗0 , ε∗1 } where ε∗0 , ε∗1 are defined √1 in (74), (75) and are dependent on the predefined radius Rric = we can apply 4∥G∥op

κ(G)MJ

Lemma 3 to show that the Riccati equation accepts a unique solution that lives inside the closed ball β U0 BR = {U : ∥U − Ũ0 ∥ ≤ Rric }. This implies that for the selected range of the ratio α we have : ric ∥U ∥op ≤ Rric + ∥U0 ∥op We now apply the transformation to the linear recursion. Recall that x bk E [x̃k ] = . y E [ŷk ] bk Using H = T −1 ΠT , the recursion (34) becomes x bk+1 I −U I − αJxx − βU Jyx = byk+1 0 I −βJyx Define the transformed bias

rkx rky

≜T

bxk byk

0 I − β(Jyy − Jyx U )

I 0

U I

bxk byk

I 0

U I

bxk byk

+ Rk

Multiplying the bias recursion by T gives x x α,β rk+1 I − αJxx 0 rk = + T Rk , y α,β rk+1 rky −βJyx I − βJyy where β α,β Jxx ≜ Jxx + U Jyx α β α,β −1 Jyy ≜ Jyy − Jyx Jxx Jxy − Jyx Ûα,β . α In particular, the leading-order matrix for the slow variable’s bias is

(38)

−1 ∆ ≜ Jyy − Jyx Jxx Jxy . This matrix is the effective first-order drift of the slow dynamics. The assumption that −∆ is Hurwitz is standard in two-time-scale stochastic approximation and appears, for instance, in the linear and CLT analyses of TTSA [26, 17, 18, 21].

For lighter notation, we write (38) as where rk ≜

rkx rky

rk+1 = Πrk + T Rk , α,β I − αJxx 0 , and Π ≜ . α,β −βJyx I − βJyy 23

Block diagonalization through a Sylvester equation. After the Riccati transformation, the linear part of the recursion is block lower triangular. The remaining coupling is the bottom left block −βJyx , which still lets the fast transformed bias influence the slow transformed bias. We now remove this coupling by a second change of coordinates. Let I 0 I 0 S≜ , S −1 = , (39) X I −X I where X ∈ Rdy ×dx is to be determined. A direct computation gives α,β I − αJxx SΠS −1 = α,β α,β X(I − αJxx ) − (I − βJyy )X − βJyx

0 α,β I − βJyy

Hence, the bottom left block vanishes if and only if α,β α,β X(I − αJxx ) − (I − βJyy )X − βJyx = 0.

Equivalently, after simplifying and dividing by α, α,β XJxx −

β α,β β Jyy X = − Jyx . α α

(40)

This is a Sylvester equation for the unknown block X [22, Chapter 2, Eq 4.3], [32, Lemma 8.3]. From β < min{ε∗0 , ε∗1 , ε∗2 , ε∗3 } that are dependent on Rric [32, Lemma 8.3]. Since we are assuming B.1 and α and are defined in (74), (75), (78), (79) then the above Sylvester equation accepts a unique solution −1 X that lives inside a closed ball BεRsyl = {X : ∥X∥op ≤ εRsyl } where Rsyl = 2 ∥Jyx ∥op Jxx , op this directly implies that: ∥X∥op ≤

β Rsyl α

(41)

Furthermore, the solution admits the expansion X=

∞ ℓ X β ℓ=1

X̃ℓ ,

β We can solve for X̃1 by replacing the sum into (40) then dividing by α then setting it to zero to get: −1 X̃1 = −Jyx Jxx

With this choice of X, the transformed matrix Λ ≜ SΠS −1 is block diagonal: α,β 0 I − αJxx Λ= . α,β 0 I − βJyy We now apply the second change of coordinates to the transformed bias vector rk . Define x x x sk rk rk I 0 . ≜S = y y X I rk rky sk By the choice of X, the linear part of the recursion is block diagonal. Hence, using I U ST = , X XU + I we obtain

sxk+1 syk+1

=Λ

sxk syk

−

I X

U XU + I

α E [Rf (x̃k , ŷk )] + E [ψk ] . β E [Rg (x̃k , ŷk )] + E [φk ]

Hence, the concatenated transformed bias sk satisfies the compact recursion sk+1 = Λsk + ST Rk . 24

(42)

Recursion on the transformed system. We now roll out the transformed recursion (42). Since the matrix Λ is block diagonal, we obtain sk+1 = Λk+1 s0 +

k X

Λk−j ST Rj .

j=0

This representation is the main payoff of the two transformations: the linear part of the recursion is now decoupled into a fast block and a slow block. Writing the two components separately gives α,β k+1 x sxk+1 = (I − αJxx ) s0 − α

k X α,β k−j (I − αJxx ) (E [Rf (x̃j , ŷj )] + E [ψj ]) j=0

−β

k X

α,β k−j (I − αJxx ) U (E [Rg (x̃j , ŷj )] + E [φj ]) ,

(43)

j=0 α,β k+1 y syk+1 = (I − βJyy ) s0 − α

k X α,β k−j (I − βJyy ) X (E [Rf (x̃j , ŷj )] + E [ψj ]) j=0

−β

k X

α,β k−j (I − βJyy ) (XU + I) (E [Rg (x̃j , ŷj )] + E [φj ]) .

(44)

j=0

The original bias recursion was fully coupled: the fast and slow biases influenced each other already at the level of the linear dynamics. After the changes of coordinates defined by T and S, this coupling has been pushed into the residual terms, while the leading linear dynamics are block diagonal. This is precisely what allows us to study the two components separately. The next step is to show that the two diagonal blocks are contractions in suitable Lyapunov norms. More precisely, for our choice on the step sizes, the matrices describing the perturbed dynamics: α,β I − αJxx

and

α,β I − βJyy

inherit the stability of I − αJxx and I − β∆, respectively. Their powers therefore produce the exponentially decaying initial-condition terms in the final bias bound, with decay rates of order α for the fast component and β for the slow component. Bound of the fast transformed component. We start from the recursion (43). For readability, write α,β Ax ≜ I − αJxx . Applying the triangle inequality, submultiplicativity of the operator norm, and Jensen’s inequality gives sxk+1 ≤

Ak+1 sx0 + α x

k X

Ak−j E [∥Rf (x̃j , ŷj )∥] x op

j=0

+β

k X

Ak−j ∥U ∥op E [∥Rg (x̃j , ŷj )∥] x op

j=0

+α

k X

Ak−j x E [ψj ] + β

j=0

k X

Ak−j x U E [φj ] .

j=0

In the two terms involving Rf and Rg , we used ∥E [Z]∥ ≤ E [∥Z∥]. By contrast, we keep the Markovian-noise terms as weighted sums, since these will be controlled later through the Poisson equation argument. α,β The Euclidean operator norm does not directly reveal the contraction of Ax = I − αJxx . We therefore use a Lyapunov norm adapted to the unperturbed fast linear dynamics. Since −Jxx is Hurwitz by Assumption B.1, there exists a unique symmetric positive definite matrix L solving ⊤ LJxx + Jxx L = I.

(45)

We denote the associated norm by ∥z∥L ≜

√

z ⊤ Lz.

β Since we are assuming B.1, α < α (Rric ) as defined in (81) and the ratio satisfies α < ∗ ∗ ∗ ∗ ∗ ∗ min{ε0 , ε1 , ε5 } where ε0 , ε1 , ε5 depend on Rric and are fully defined in (74), (75), (82) then after applying submultiplicativity of the norm we can apply Lemma 5 to get: we have: ∗

(I − αJxx − βU Jyx )k L ≤ ∥I − αJxx − βU Jyx ∥L ≤ e−kcα where c = (8 ∥L∥op )−1 , see Lemma 5 and the proof for full derivation. Equivalently, Akx L ≤ e−kcα which yields for every ℓ ≥ 0, Aℓx op ≤

p p κ(L) Aℓx L ≤ κ(L)e−cαℓ .

Applying this bound to (43) gives sxk+1 ≤

κ(L)e−(k+1)cα ∥sx0 ∥ + α

k X p κ(L) e−cα(k−j) E [∥Rf (x̃j , ŷj )∥]

(46)

j=0

+β

κ(L) ∥U ∥op

k X

e−cα(k−j) E [∥Rg (x̃j , ŷj )∥]

(47)

j=0

+α

k X

Ak−j x E [ψj ] + β

j=0

k X

Ak−j x U E [φj ] .

(48)

j=0

We now bound the terms on the right hand side separately. The first term in (46) is the exponentially decaying contribution of the initial condition. The second term involves the Taylor remainder. Since Rf is a second order remainder, there exists a constant LR > 0 such that: h 2 i E [∥Rf (x̃j , ŷj )∥] ≤ LR E ∥x̃j ∥ + ∥ŷj ∥ h i 2 2 ≤ 2LR E ∥x̃j ∥ + ∥ŷj ∥ h i 2 2 = 2LR E ∥xj − h(yj ) + h(yj ) − h(y ∗ )∥ + ∥ŷj ∥ h i 2 2 2 ≤ 2LR E 2 ∥xj − h(yj )∥ + 2 ∥h(yj ) − h(y ∗ )∥ + ∥ŷj ∥ h i 2 2 ≤ 2LR E 2 ∥xj − h(yj )∥ + (2L2h + 1) ∥ŷj ∥ . Setting: d1 ≜ max{2, 2L2h + 1}, h i 2 2 E [∥Rf (x̃j , ŷj )∥] ≤ 2LR d1 E ∥xj − h(yj )∥ + ∥ŷj ∥ .

we obtain:

Plugging this into the second term of (46) yields: α

κ(L)

k X

e−cα(k−j) E [∥Rf (x̃j , ŷj )∥] ≤ 2α

k h i X p 2 2 κ(L)LR d1 e−cα(k−j) E ∥xj − h(yj )∥ + ∥ŷj ∥ .

j=0

By Lemma 8, and using the notation ∞ Mα,β ≜ MSE∞ x (α, β) + MSEy (α, β),

we obtain k X p p c,α α κ(L) e−cα(k−j) E [∥Rf (x̃j , ŷj )∥] ≤ α 2 κ(L)LR d1 (Mα,β + uc,α k + wk ) | {z } j=0

≜Cx,1

+2 |

p κ(L)LR d1 c−1 Mα,β . {z } ≜Cx,2

The term involving Rg is controlled in the same way. We denote MU ≜ Rric + ∥U0 ∥op . Thus, (47) ≤ β 2 |

β p c,α κ(L)MU LR d1 (Mα,β + uc,α 2 κ(L)MU LR d1 c−1 Mα,β . k + wk ) + {z } {z } α| Cx,4

≜Cx,3

It remains to control the two weighted Markovian-noise terms in (48). By Lemma 7, β2 B4 . α where the constants are defined in detail in the proof of the lemma. Merging all the bounds from above we obtain p c,α sxk+1 ≤ e−(k+1)cα κ(L) ∥sx0 ∥ + (α + β)Fx(1) + Fx(2) Mα,β + (α + β)Fx(3) (uc,α k + wk ) (48) ≤ αA1 + β(A2 + B1 ) + α2 A3 + β 2 B2 + αβ(A4 + B3 ) +

x + Rα,β (1)

(2)

(3)

x where the quantities Fx , Fx , Fx , Rα,β are defined by:

Fx(1) := A1 ∨ (A2 + B1 ) ∨ ω ∗ B4 , Fx(2) := Cx,2 + ω ∗ Cx,4 + Fx(3) (α∗ + β ∗ ) , Fx(3) := Cx,3 ∨ Cx,1 , x := α2 A3 + β 2 B2 + αβ(A4 + B3 ) . Rα,β

Bound of the slow transformed component.

We now turn to the slow component. Let I ′ ≜ I + XU.

α,β Ay ≜ I − βJyy ,

Applying the triangle inequality, submultiplicativity, and Jensen’s inequality to (44) gives syk+1 ≤ Ak+1 sy0 + α y

k X

Ak−j ∥X∥op E [∥Rf (x̃j , ŷj )∥] y op

j=0

+β

k X

Ak−j ∥I ′ ∥op E [∥Rg (x̃j , ŷj )∥] y op

j=0

+α

k X

Ak−j XE [ψj ] + β y

j=0

k X

Ak−j I ′ E [φj ] . y

j=0

We use a Lyapunov norm adapted to the effective slow matrix −1 ∆ ≜ Jyy − Jyx Jxx Jxy .

Since −∆ is Hurwitz by Assumption B.1, there exists a unique symmetric positive definite matrix G solving the following Lyapunov equation: G∆ + ∆⊤ G = I.

(49)

β Since we are assuming B.1, and we have β ≤ β ∗ as defined in (84) that also satisfies α < min{ε∗0 , ε∗1 } then we can apply Lemma 6 to have:

∥I − β(∆ − Jyx (U − U0 ))∥G ≤ e−dβ where d = (8 ∥G∥op )−1 , see Lemma 6 and its proof for full derivation. Equivalently, Aℓy G ≤ e−dβℓ ,

∀ℓ ≥ 0.

we obtain, for every ℓ ≥ 0, Aℓy op ≤

p p κ(G) Aℓy G ≤ κ(G)e−dβℓ . 27

Therefore, syk+1 ≤ e−dβ(k+1)

+ α ∥X∥op

κ(G) ∥sy0 ∥

κ(G)

k X

e−dβ(k−j) E [∥Rf (x̃j , ŷj )∥]

(50)

j=0

+ β ∥I ′ ∥op

κ(G)

k X

e−dβ(k−j) E [∥Rg (x̃j , ŷj )∥]

(51)

j=0

+α

k X

Ak−j XE [ψj ] + β y

k X

Ak−j I ′ E [φj ] . y

(52)

j=0

We now bound the Taylor remainder terms. As in the fast component analysis, there exists a constant LR > 0 such that h i 2 2 E [∥Rf (x̃j , ŷj )∥] ≤ 2LR d1 E ∥xj − h(yj )∥ + ∥ŷj ∥ , d1 ≜ max{2, 2L2h + 1}. The same bound holds for Rg , up to increasing LR if necessary. By Lemma 8, k X

−dβ(k−j)

E [∥Rf (x̃j , ŷj )∥] ≤ 2LR d1

j=0

1 1+ dβ

d,β Mα,β + ud,β k + wk

At first sight, the factor 1/β is problematic, however, the term (50) still has a multiplication by α ∥X∥op . Since we are assuming B.1, we showed previously in (41) that α ∥X∥op ≤ βRsyl −1 with Rsyl = 2 ∥Jyx ∥op Jxx . Consequently, op p d,β (50) ≤ β 2Rsyl κ(G)LR d1 Mα,β + ud,β + w k k | {z } ≜Cy,1

+ 2Rsyl |

p κ(G)LR d1 d−1 Mα,β . {z } ≜Cy,2

β < min{ε∗0 , ε∗1 , ε∗2 , ε∗3 } we have The term (51) is controlled similarly. For α

∥I ′ ∥op ≤ ∥I∥op + ∥X∥op ∥U ∥op ≤ 1 + Rsyl (Rric + ∥U0 ∥op ) ≜ MI ′ Lemma 8 gives (51) ≤ β 2MI ′ |

p d,β κ(G)LR d1 Mα,β + ud,β + w k k {z } ≜Cy,3

+ 2MI ′ |

p κ(G)LR d1 d−1 Mα,β . {z } ≜Cy,4

It remains to bound the weighted Markovian-noise terms in (52). By Lemma 7, (52) ≤ α(C2 + D4 ) + β(C1 + D1 ) + αβ(C4 + D3 ) + β 2 (C3 + D2 ). We define : Fy(1) := (C2 + D4 ) ∨ (C1 + D1 ) , Fy(2) := Cy,2 + Cy,4 + β ∗ (Cy,1 + Cy,3 ) , Fy(3) := Cy,3 + Cy,1 , y := αβ(C4 + D3 ) + β 2 (C3 + D2 ) . Rα,β

We wrap all the bounds in the final one p d,β syk+1 ≤ e−dβ(k+1) κ(G) ∥sy0 ∥ + (α + β)Fy(1) + Fy(2) Mα,β + βFy(3) (ud,β k + wk ) y + Rα,β

we merge the bounds of ∥sxk ∥ and ∥syk ∥ to get p p sxk+1 + syk+1 ≤e−(k+1)cα κ(L) ∥sx0 ∥ + e−dβ(k+1) κ(G) ∥sy0 ∥ c,α + Fx(1) + Fy(1) (α + β) + (Fx(2) + Fy(2) ) Mα,β + (α + β)Fx(3) (uc,α k + wk ) | {z } {z } | ≜F1

≜F2

d,β y x + βFy(3) (ud,β k + wk ) + Rα,β + Rα,β {z } | ≜Rα,β

We define the new sequence ŜMSE to be: k d,β c,α (3) c,α ŜMSE = βFy(3) (ud,β k k + wk ) + (α + β)Fx (uk + wk )

where we have showed that by Lemma 8 that the limit exist and it is zero lim ŜMSE =0 k

k→∞

which yields our bound: sxk+1 + syk+1 ≤e−(k+1)cα

p p κ(L) ∥sx0 ∥ + e−dβ(k+1) κ(G) ∥sy0 ∥

+ F1 (α + β) + F2 Mα,β + ŜMSE + Rα,β k c,α where Rα,β contains terms of higher order O((α + β)2 ). Moreover, from the definition of uc,α k , wk in Lemma 8 there exists a, b, ĉ, dˆ such that: 2

ŜMSE ≤ ae−αĉk ∥x̂0 ∥ + be−β dk ∥ŷ0 ∥ k

(53)

So our final bound will be sxk+1 + syk+1 ≤e−(k+1)cα

p p ˆ 2 2 κ(L) ∥sx0 ∥ + e−dβ(k+1) κ(G) ∥sy0 ∥ + ae−αĉk ∥x̂0 ∥ + be−β dk ∥ŷ0 ∥

+ F1 (α + β) + F2 Mα,β + Rα,β Unwrapping the results.

It remains to return to the original bias coordinates. By construction, x x sk bk sk = = ST . syk byk

Since the matrices S and T are invertible, we can write equivalently, x bk = T −1 S −1 sk . byk Therefore, by applying the norm on both sides, x bk ≤ T −1 S −1 sk byk ≤ T −1 S −1 op ∥sk ∥ ≤ T −1 S −1 op (∥sxk ∥ + ∥syk ∥) . where in the last inequality we have used the triangular inequality since x x sk 0 sk = + ≤ ∥sxk ∥ + ∥syk ∥ 0 syk syk Plugging the bounds obtained above for the fast and slow transformed components gives the desired bias bound: max {∥bxk ∥ , ∥byk ∥} ≤ T −1 S −1 op (∥sxk ∥ + ∥syk ∥) . which concludes the proof of our theorem. 29

Supporting Lemmas

The first lemma is merely a collection of various Lipschitzness/boundedness properties induced by the assumptions we make in the paper, introducing several constants that we use in the proofs of the main theorems. Lemma 1. Under assumptions A.1-4, there exist finite non-negative constants Mf , Mg , Mh , Mξ , MS dependent on the problem settings satisfying : 1. ∥f (xk , yk )∥ ≤ 2Lf MS + ∥f (0, 0)∥ ≤ Mf 2. ∥g(xk , yk )∥ ≤ 2Lg MS + ∥g(0, 0)∥ ≤ Mg 3. ∥f (xk , yk , ξk+1 )∥ ≤ 2Lf MS + MS ≤ Mf 4. ∥g(xk , yk , ξk+1 )∥ ≤ 2Lg MS + MS ≤ Mg 5. ∥∆hk ∥ ≤ βLh Mg ≤ βMh 6. ∥ψk ∥ ≤ Mξ 7. ∥h(yk )∥ ≤ Lh MS + ∥h(y0 )∥ ≤ Mh 8. Given assumption A.2-4, for every x1 , x2 ∈ Rdx , y1 , y2 ∈ Rdy , vf and vg (as defined in (54) and (65)) are Lvf and Lvg -Lipschitz continuous: ∥vf (x1 , y1 , ξ) − vf (x2 , y2 , ξ)∥ ≤ Lvf (∥x1 − x2 ∥ + ∥y1 − y2 ∥) ∥vg (x1 , y1 , ξ) − vg (x2 , y2 , ξ)∥ ≤ Lvg (∥x1 − x2 ∥ + ∥y1 − y2 ∥) 9. ∥vf (x1 , y1 , ξ)∥ ≤ 2Lf MS + f0 ≤ Mξ 10. ⟨x̂k , vf (xk , yk , ξk+1 )⟩ ≤ Mξ Proof. We provide one line proof for each point: 1. ∥f (xk , yk )∥ ≤ ∥f (xk , yk ) − f (0, 0)∥ + ∥f (0, 0)∥ ≤ Lf (∥xk ∥ + ∥yk ∥) + ∥f (0, 0)∥ ≤ 2Lf MS + ∥f (0, 0)∥ ≤ Mf 2. ∥g(xk , yk )∥ ≤ 2Lg MS + ∥g(0, 0)∥ ≤ Mg 3. ∥f (xk , yk , ξk+1 )∥ ≤ 2Lf MS + MS ≤ Mf 4. ∥g(xk , yk , ξk+1 )∥ ≤ 2Lg MS + MS ≤ Mg 5. ∥∆hk ∥ = ∥h(yk+1 ) − h(yk ))∥ ≤ Lh ∥βg(xk , yk , ξk+1 )∥ ≤ βMh 6. ∥ψk ∥ = ∥f (xk , yk , ξk+1 ) − f (xk , yk )∥ ≤ ∥f (xk , yk , ξk+1 )∥ + ∥f (xk , yk )∥ ≤ 2Mf ≤ Mξ 7. ∥h(yk )∥ ≤ Lh MS + ∥h(y0 )∥ ≤ Mh 8. See [1, Lemma 8] and [17, Lemma D.4]. 9. ∥vf (x1 , y1 , ξ)∥ ≤ 2Lf MS + f0 ≤ Mξ 10. ⟨x̂k , vf (xk , yk , ξk+1 )⟩ ≤ ∥x̂k ∥ ∥vf (xk , yk , ξk+1 )∥ ≤ Mξ

The next lemma proves bounds on the terms containing the weighted sum of deviations of the dynamics estimators from their stationary averages f (xk , yk , ξk+1 )−f (xk , yk ) and g(xk , yk , ξk+1 )− g(xk , yk ) in Theorem 1. 30

Lemma 2. Taking the sequence (xk , yk )k≥0 defined by (1) and under assumptions A.1-3, for every αµf < 1 and βµg < 1:   k X Lvf Mf (MS + Mh ) 3 − αµf Mf (1 + |Ξ|MS LK ) k−j   αE (1 − αµf ) ⟨x̂j , ψj ⟩ ≤α Mξ + + 1 − αµf µf µf j=0 Mξ Mh + Lvf Mg (MS + Mh ) + |Ξ|MS LK Mξ Mg +β µf   k X (Lvg + Mξ LK |Ξ|)MS Mg 3 − βµg Mg k−j   βE (1 − βµg ) ⟨ŷj , ϕj ⟩ ≤β + Mξ + µg 1 − βµg µg j=0 +α

(Lvg + Mξ LK |Ξ|)MS Mf µg

where ψk = f (xk , yk , ξk+1 ) − f (xk , yk ) and ϕk = g(xk , yk , ξk+1 ) − g(xk , yk ). Proof. To prove this lemma, we use the same technique in [1, Lemma 8] and [21] by invoking the solution of a Poisson equation that allow us to decompose the deviation terms f (xj , yj , ξj+1 ) − f (xj , yj ) inside the sum into a Martingale Difference Sequence (MDS), a bounded term and a telescoping term. Since our terms in the sum are weighted by (1 − αµf )k−j we don’t have exactly the telescoping terms that we want, instead, we have additional residuals that we needed to deal with. Moreover, since our transition kernel K is depending on the iterates, we will have a drift between the Poisson solution time step and the required term to build the Martingale Difference. In this proof we first start by introducing the core idea of the Poisson equation, then we derive each component individually to get their respective bounds. Moreover, throught this lemma we bound images of the functions g and ḡ by continuity and compactness and not using A. 4 and we use the same upper bound constant Mg . Under assumption A.2, by [12, Definition 21.2.1] and [1, Lemma 8] vf is a solution of the following Poisson equation: vf (xk , yk , ξk+1 ) −

K(xk ,yk ) (ξk+1 , ξ ′ )vf (xk , yk , ξ ′ ) = f (xk , yk , ξk+1 ) − f (xk , yk )

(54)

ξ ′ ∈Ξ

In the following, we use the notation, with Kk ≜ K(xk ,yk ) X Eξ′ ∼Kk (ξk+1 ,·) [vf (xk , yk , ξ ′ )] ≜ K(xk ,yk ) (ξk+1 , ξ ′ )vf (xk , yk , ξ ′ ). ξ ′ ∈Ξ

It is shown in [1, Lemma 8] that vf is Lipschitz-continuous. Our goal from introducing this solution is to be able to write this decomposition: k X

(1 − αµf )k−j ⟨x̂j , ψj ⟩ =

j=0

k X

MDS + Telescoping Term+ Bounded Term + Kernel Drift

j=0

First, we start by substituting (54) in the sum: k X j=0

(1 − αµf )k−j ⟨x̂j , ψj ⟩ =

k X

(1 − αµf )k−j x̂j , vf (xj , yj , ξj+1 ) − Eξ′ ∼Kj (ξj+1 ,·) [vf (xj , yj , ξ ′ )]]

j=0

We notice that the term Eξ′ ∼Kj (ξj+1 ,·) [vf (xj , yj , ξ ′ )] does not make a Martingale Difference with the term vf (xj , yj , ξj+1 ) since ξ ′ is sampled at the step j + 2, we add and subtract vf (xj , yj , ξj+2 ): vf (xj , yj , ξj+1 ) − Eξ′ ∼Kj (ξj+1 ,·) [vf (xj , yj , ξ ′ )] =vf (xj , yj , ξj+2 ) − Eξ′ ∼Kj (ξj+1 ,·) [vf (xj , yj , ξ ′ )] + vf (xj , yj , ξj+1 ) − vf (xj , yj , ξj+2 ) However, this is not enough since ξj+2 ∼ K(xj+1 ,yj+1 ) (ξj+1 , ·) and not K(xj ,yj ) (ξj+1 , ·). To obtain a Martingale Difference we add and subtract Eξ′ ∼Kj+1 (ξj+1 ,·) [vf (xj , yj , ξ ′ )]. Finally, to use 31

Lipschitzness of vf and get a telescoping term we add and subtract vf (xj+1 , yj+1 , ξj+2 ). This will give us the following decomposition of the Markovian noise deviation D ⟨x̂j , ψj ⟩ = x̂j ,vf (xj , yj , ξj+2 ) − Eξ′ ∼Kj+1 (ξj+1 ,·) [vf (xj , yj , ξ ′ )] (55) + vf (xj+1 , yj+1 , ξj+2 ) − vf (xj , yj , ξj+2 ) + vf (xj , yj , ξj+1 ) − vf (xj+1 , yj+1 , ξj+2 )

(56) (57)

E + Eξ′ ∼Kj+1 (ξj+1 ,·) [vf (xj , yj , ξ ′ )] − Eξ′ ∼Kj (ξj+1 ,·) [vf (xj , yj , ξ ′ )]

(58)

the first term (55) is an MDS, so it reduces to zero under expectation:   k X E  (1 − αµf )k−j x̂j , vf (xj , yj , ξj+2 ) − Eξ′ ∼Kj+1 (ξj+1 ,·) [vf (xj , yj , ξ ′ )]  = 0 j=0

The second term (56) scales linearly with the sum of step-sizes, since f and K are Lipschitz, so vf will be Lvf -Lipschitz too for some positive constant Lvf and this is due to Lemma 1. So we can bound (56) by: |(56)| ≤ ∥x̂j ∥ ∥vf (xj+1 , yj+1 , ξj+2 ) − vf (xj , yj , ξj+2 )∥ = ∥xj − h(yj )∥ ∥vf (xj+1 , yj+1 , ξj+2 ) − vf (xj , yj , ξj+2 )∥ ≤Lvf ∥xj − h(yj )∥ (∥xj+1 − xj ∥ + ∥yj+1 − yj ∥) ≤Lvf (∥xj ∥ + ∥h(yj )∥)(∥xj+1 − xj ∥ + ∥yj+1 − yj ∥) ≤Lvf (MS + Mh )(∥αf (xj , yj , ξj+1 )∥ + ∥βg(xj , yj , ξj+1 )∥) ≤Lvf (MS + Mh )(αMf + βMg ) For brevity we set zj = ⟨x̂j , vf (xj+1 , yj+1 , ξj+2 ) − vf (xj , yj , ξj+2 )⟩. Now we apply our bound on the sum to have α

k X

(1 − αµf )k−j zj ≤α

j=0

k X (1 − αµf )k−j |zj | j=0

≤α

k X (1 − αµf )k−j Lvf (MS + Mh )(αMf + βMg ) j=0

=αLvf (MS + Mh )(αMf + βMg )

1 − (1 − αµf )k+1 µf α

Lvf (MS + Mh ) (αMf + βMg ) µf To see why the term (58) is bounded we use our assumptions A 2 and A 3 about Lipschitzness of K and f ≤

(58) =

≤

(Kj+1 (ξj+1 , ξ ′ ) − Kj (ξj+1 , ξ ′ ))vf (xj , yj , ξ ′ )

ξ ′ ∈Ξ

∥Kj+1 (ξj+1 , ξ ′ ) − Kj (ξj+1 , ξ ′ )∥ ∥vf (xj , yj , ξ ′ )∥

ξ ′ ∈Ξ

≤

LK (∥xj+1 − xj ∥ + ∥yj+1 − yj ∥) Mξ

ξ ′ ∈Ξ

≤ LK Mξ (αMf + βMg )|Ξ| which gives the bound α

k X

(1 − µf α)k−j x̂j , Eξ′ ∼Kj+1 (ξj+1 ,·) [vf (xj , yj , ξ ′ )] − Eξ′ ∼Kj (ξj+1 ,·) [vf (xj , yj , ξ ′ )]

j=0

≤

|Ξ|MS LK Mξ (αMf + βMg ) µf

Finally, we handle the third term (57) which is supposed to be a telescoping term inside the sum: k X α (1 − αµf )k−j ⟨x̂j , vf (xj , yj , ξj+1 ) − vf (xj+1 , yj+1 , ξj+2 )⟩

(59)

j=0

This is not perfectly telescoping. This hurdle emerges from two quantities: first, in the inner product itself we have x̂j which depends on j and second, the weights (1 − αµf )k−j also depend on j. To get over this, we will proceed in two steps. First, we fix the telescoping term inside the inner product by using the update rule in (1). Due to this step, we will have a remainder term dependent on the dynamics of xj and controlled by α. The second step, is to handle the weights to have our telescoping terms. These two steps, will introduce remainders that we will need to control. Now we detail every step, so we start with the inner product and we set vj = vf (xj , yj , ξj+1 ) so we can easily write the computations: ⟨x̂j , vj+1 − vj ⟩ = ⟨xj − h(yj ), vj+1 ⟩ − ⟨x̂j , vj ⟩ = ⟨xj+1 + αf (xj , yj , ξj+1 ) − h(yj ), vj+1 ⟩ − ⟨x̂j , vj ⟩ = x̂j+1 + αf (xj , yj , ξj+1 ) + ∆hj , vj+1 − ⟨x̂j , vj ⟩ = ⟨x̂j+1 , vj+1 ⟩ − ⟨x̂j , vj ⟩ + αf (xj , yj , ξj+1 ) + ∆hj , vj+1 where ∆hj = h(yj+1 ) − h(yj ). Now, we have a telescoping term plus a remainder. Substituting the last line into (59): (59) = α

k X

(1 − αµf )k−j (⟨x̂j+1 , vj+1 ⟩ − ⟨x̂j , vj ⟩)

(60)

j=0

+ α2

k X (1 − αµf )k−j ⟨f (xj , yj , ξj+1 ), vj+1 ⟩ j=0

+α

k X

(1 − αµf )k−j ∆hj , vj+1

j=0

Starting with the first term (60). Because of the weights in the sum (1 − αµf )k−j we can’t directly have a telescoping sum. An extra step has to be done for this term to be able to make it telescope with the cost of introducing new remainders. Hopefully, these remainders can be controlled. We proceed by: (60) =α

k k X X (1 − αµf )k−j ⟨x̂j+1 , vj+1 ⟩ − α (1 − αµf )k−j ⟨x̂j , vj ⟩ j=0

j=0

=α(1 − αµf )

k X

(1 − αµf )k−(j+1) ⟨x̂j+1 , vj+1 ⟩ − α

j=0

k X (1 − αµf )k−j ⟨x̂j , vj ⟩ j=0

k k X X =α (1 − αµf )k−(j+1) ⟨x̂j+1 , vj+1 ⟩ − α (1 − αµf )k−j ⟨x̂j , vj ⟩ j=0

− α2 µf

j=0 k X

(1 − αµf )k−(j+1) ⟨x̂j+1 , vj+1 ⟩

j=0

=α(1 − αµf )−1 ⟨x̂k+1 , vk+1 ⟩ − α(1 − αµf )k ⟨x̂0 , v0 ⟩ − α2 µf

k X

(1 − αµf )k−(j+1) ⟨x̂j+1 , vj+1 ⟩

j=0

where in the second equality we factored out 1 − αµf to make the inside weight match the weight in the other sum. In the last line, we removed the telescoping terms to end up with the final three terms. 33

Now, into (59) we apply absolute value to get: |(59)| ≤ α (1 − αµf )−1 ⟨x̂k+1 , vk+1 ⟩ − (1 − αµf )k ⟨x̂0 , v0 ⟩ +α

k X

(1 − αµf )k−j

(61)

∆hj , vj+1

(62)

(1 − αµf )k−(j+1) |⟨x̂j+1 , vj+1 ⟩|

(63)

j=0

+ α2 µf

k X j=0

+ α2

k X (1 − αµf )k−j |⟨f (xj , yj , ξj+1 ), vj+1 ⟩|

(64)

j=0

We bound each term independently. Starting with the first term (61): 2 − αµf Mξ + αMξ = αMξ (61) ≤ α 1 − αµf 1 − αµf the second term (62) (62) ≤ α

k X

(1 − αµf )k−j ∥∆hj ∥∥vj+1 ∥ ≤ αβMξ Mh

j=0

j=0 k+1

= αβMξ Mh

k X (1 − αµf )j

1 − (1 − αµf ) αµf

≤β

Mξ Mh µf

the third term (63) (63) ≤ α2 µf Mξ

k k X X (1 − αµf )k−(j+1) = α2 µf Mξ (1 − αµf )−1 (1 − αµf )k−j j=0

j=0

Mξ 1 − αµf and the last term (64) ≤α

(64) ≤ α2 Mξ Mf

k X Mξ Mf (1 − αµf )j ≤ α µf j=0

After bounding each term individually, using bounds from Lemma 1, we apply these bounds on the main sum (59) 3 − αµf Mf Mξ Mh |(59)| ≤ αMξ + +β 1 − αµf µf µf Since expectation preserves order, we apply it on both sides to get the final result:   k X Lvf Mf (MS + Mh ) Mf (1 + |Ξ|MS LK ) 3 − αµf k−j   αE (1 − αµf ) ⟨x̂j , ψj ⟩ ≤α Mξ + + 1 − αµf µf µf j=0 Mξ Mh + Lvf Mg (MS + Mh ) + |Ξ|MS LK Mξ Mg +β µf which concludes the proof of the first part. As above, let vg be the solution of the Poisson equation: g(xj , yj , ξj+1 ) − g(xj , yj ) = vg (xj , yj , ξj+1 ) − Eξ′ ∼Kj (ξj+1 ,·) [vg (xj , yj , ξ ′ )] (65) we rewrite: ⟨ŷj , ϕj ⟩ = ⟨ŷj , vg (xj , yj , ξj+2 ) − Eξ′ ∼Kj+1 (ξj+1 ,·) [vg (xj , yj , ξ ′ )] (66) + vg (xj+1 , yj+1 , ξj+2 ) − vg (xj , yj , ξj+2 ) + vg (xj , yj , ξj+1 ) − vg (xj+1 , yj+1 , ξj+2 )

(67) (68)

+ Eξ′ ∼Kj+1 (ξj+1 ,·) [vg (xj , yj , ξ ′ )] − Eξ′ ∼Kj (ξj+1 ,·) [vg (xj , yj , ξ ′ )]⟩

(69)

the first term (66) is 0 under expectation because it is a Martingale Difference. The last term (69) is treated in the same way as in (58) and it has the following bound β

k X

(1 − µg β)k−j ŷj , Eξ′ ∼Kj+1 (ξj+1 ,·) [vg (xj , yj , ξ ′ )] − Eξ′ ∼Kj (ξj+1 ,·) [vg (xj , yj , ξ ′ )]

j=0

≤

LK Mξ MS |Ξ|(αMf + βMg ) µg

The second term (67) is of order O(α + β): |(67)| ≤ ∥ŷj ∥ ∥vg (xj+1 , yj+1 , ξj+2 ) − vg (xj , yj , ξj+2 )∥ ≤Lvg MS (∥yj+1 − yj ∥ + ∥xj+1 − xj ∥) ≤Lvg MS (∥αf (xj , yj , ξj+1 )∥ + ∥βg(xj , yj , ξj+1 )∥) ≤Lvg MS (αMf + βMg ) where in the third line we have used the definition of the updates and in the last line we have used the bounds from Lemma 1. So we have: β

k k X X (1 − βµg )k−j ⟨ŷj , vg (xj+1 , yj+1 , ξj+2 ) − vg (xj , yj , ξj+2 )⟩ ≤β (1 − βµg )k−j |(67)| j=0

j=0

≤

Lvg MS (αMf + βMg ) µg

where in the first inequality we have used the triangular inequality. Moving to the third term (68). Let wj = vg (xj , yj , ξj+1 ), the term in the sum becomes: −(68) = ⟨ŷj , wj+1 − wj ⟩ = ⟨ŷj , wj+1 ⟩ − ⟨ŷj , wj ⟩ = ⟨ŷj+1 + βg(xj , yj , ξj+1 ), wj+1 ⟩ − ⟨ŷj , wj ⟩ = ⟨ŷj+1 , wj+1 ⟩ − ⟨ŷj , wj ⟩ + ⟨βg(xj , yj , ξj+1 ), wj+1 ⟩ We substitute in the sum: β

k k X X (1 − βµg )k−j ⟨ŷj , wj+1 − wj ⟩ = β (1 − βµg )k−j (⟨ŷj+1 , wj+1 ⟩ − ⟨ŷj , wj ⟩) j=0

(70)

j=0

+ β2

k X (1 − βµg )k−j ⟨g(xj , yj , ξj+1 ), wj+1 ⟩ j=0

treating the first term (70): (70) =β(1 − βµg )−1 ⟨ŷk+1 , wk+1 ⟩ − β(1 − βµg )k ⟨ŷ0 , w0 ⟩ 2

− β µg

k X

(1 − βµg )k−(j+1) ⟨ŷj+1 , wj+1 ⟩

j=0

substituting: β

k X

(1 − βµg )k−j ⟨ŷj , wj+1 − wj ⟩ =β (1 − βµg )−1 ⟨ŷk+1 , wk+1 ⟩ − (1 − βµg )k ⟨ŷ0 , w0 ⟩

j=0

(71) + β 2 µg

k X (1 − βµg )k−(j+1) ⟨ŷj+1 , wj+1 ⟩

(72)

j=0

+ β2

k X (1 − βµg )k−j ⟨g(xj , yj , ξj+1 ), wj+1 ⟩ j=0

(73)

we bound the term (71): (71) ≤ βMξ

2 − βµg 1 − βµg

the second term (72): (72) ≤ β

Mξ 1 − βµg

and the last term (73): (73) ≤ β

Mξ Mg µg

where we have used the bounds from Lemma 1 as we have done in detail for the previous section regarding the fast iterate. The last bound can be written as:   k X (Lvg + Mξ LK |Ξ|)MS Mg 3 − βµg Mg k−j   βE (1 − βµg ) ⟨ŷj , ϕj ⟩ ≤β + Mξ + µg 1 − βµg µg j=0 +α

(Lvg + Mξ LK |Ξ|)MS Mf µg

which concludes the proof. The following two lemmas are describing a uniform upper bounds of the solutions of the Riccati β equation (37) and Sylvester equation (40) for a range of values of α . Lemma 3. Let ε > 0 and let Jyx ∈ Rm×n , Jxx ∈ Rn×n , Jyy ∈ Rm×m , Jxy ∈ Rn×m with Jxx being invertible, for every Rric > 0 we define ε∗0 , ε∗1 as follows : ε∗0 = ε∗1 =

Rric , −1 Jxx op (Rric + ∥U0 ∥op )(∥Jyy ∥op + (Rric + ∥U0 ∥op ) ∥Jyx ∥op ) 1 −1 Jxx op

(74) (75)

∥Jyy ∥op + 2 ∥Jyx ∥op (Rric + ∥U0 ∥op )

−1 with U0 = Jxx Jxy . Then for every ε < min{ε∗0 , ε∗1 } there exists a matrix U ∈ Rn×m that belongs U0 to the closed ball BR = {U : ∥U − U0 ∥op ≤ Rric } that is the unique solution of the following ric Riccati equation:

εU Jyx U + Jxx U − εU Jyy = Jxy

(76)

−1 Proof. We multiply the Riccati equation (76) by Jxx ,since we assumed it is invertible: −1 −1 −1 εJxx U Jyx U + U − εJxx U Jyy = Jxx Jxy | {z } ≜U0

which gives the fixed point equation: −1 −1 U = U0 + ε Jxx U Jyy − Jxx U Jyx U

(77)

We define T to be the operator of the fixed point equation. We define the closed ball: U0 BR = {U : ∥U − U0 ∥op ≤ Rric } ric U0 U0 U0 ) ⊆ BR . Let M ∈ BR : First let’s find ε∗0 such that T (BR ric ric ric −1 −1 T (M ) = U0 + ε Jxx M Jyy − Jxx M Jyx M

calculating the norm −1 −1 ∥T (M ) − U0 ∥op = ε Jxx M Jyy − Jxx M Jyx M op −1 −1 ≤ ε( Jxx M Jyy op + Jxx M Jyx M op ) −1 −1 (M − U0 + U0 )Jyx (M − U0 + U0 ) op ) (M − U0 + U0 )Jyy op + Jxx = ε( Jxx −1 ≤ ε( Jxx (∥M − U0 ∥op + ∥U0 ∥op ) ∥Jyy ∥op op −1 (∥M − U0 ∥op + ∥U0 ∥op )2 ∥Jyx ∥op ) + Jxx op −1 −1 (Rric + ∥U0 ∥op )2 ∥Jyx ∥op ) (R + ∥U0 ∥op ) ∥Jyy ∥op + Jxx ≤ ε( Jxx op op U0 now we condition on ε to make T (M ) ∈ BR , we need ε ≤ ε∗0 such that : ric

ε∗0 =

Rric −1 Jxx op (Rric + ∥U0 ∥op )(∥Jyy ∥op + (Rric + ∥U0 ∥op ) ∥Jyx ∥op )

U0 Now we prove T is a contraction. Let U1 , U2 ∈ BR , we calculate the norm of their difference: ric −1 −1 −1 −1 U1 Jyy − Jxx U1 Jyx U1 − ε Jxx U2 Jyy − Jxx U2 Jyx U2 op ∥T (U1 ) − T (U2 )∥op = ε Jxx

We have: U1 AU1 − U2 AU2 = U1 AU1 − U1 AU2 + U1 AU2 − U2 AU2 = U1 A(U1 − U2 ) + (U1 − U2 )AU2 We replace: −1 ∥T (U1 ) − T (U2 )∥op = ε Jxx ((U1 − U2 )Jyy − ((U1 − U2 )Jyx U2 + U1 Jyx (U1 − U2 ))) op −1 ∥T (U1 ) − T (U2 )∥op ≤ ε Jxx ∥U1 − U2 ∥op ∥Jyy ∥op + ∥Jyx ∥op (∥U1 ∥op + ∥U2 ∥op ) op

Bounding the norms: ∥U1 ∥ ≤ ∥U1 − U0 ∥ + ∥U0 ∥ ≤ Rric + ∥U0 ∥op ∥U2 ∥ ≤ ∥U2 − U0 ∥ + ∥U0 ∥ ≤ Rric + ∥U0 ∥op −1 ∥T (U1 ) − T (U2 )∥ ≤ ε Jxx ∥J ∥ + 2 ∥J ∥ (R + ∥U ∥ ) ∥U1 − U2 ∥op yy yx ric 0 op op op op | {z } r(ε)

To make r(ε) < 1 we need: ε<

1 = ε∗1 −1 Jxx op ∥Jyy ∥op + 2 ∥Jyx ∥op (Rric + ∥U0 ∥op )

U0 For ε < min{ε∗0 , ε∗1 , }, the solution U of the Riccati equation is unique and is inside the ball BR . ric It can be obtained by successive application of T on U0 . This shows that the uniform bound on ∥U ∥ ≤ Rric + ∥U0 ∥op is proved.

Lemma 4. Let ε > 0 and let Jyx ∈ Rm×n , Jxx ∈ Rn×n , Jyy ∈ Rm×m , Jxy ∈ Rn×m with Jxx −1 being invertible, define U0 = Jxx Jxy , and for all Rric > 0 there exists ε∗ such that for all ε < ε∗ n×m −1 the matrix Uε ∈ R is bounded by ∥Uε ∥op ≤ Rric + ∥U0 ∥op , we set Rsyl = 2 ∥Jyx ∥op Jxx op ∗ ∗ and we define ε2 , ε3 : ε∗2 = ε∗3 =

1 ∥Jyy ∥op

−1 Jxx + 2(Rric + ∥U0 ∥op ) ∥Jyx ∥op op

−1 Jxx op

1 −1 Jxx op

∥Jyy ∥op + 2 ∥Jyx ∥op (Rric + ∥U0 ∥op )

(78) (79)

then for all ε < min{ε∗ , ε∗2 , ε∗3 } there exists a matrix X ∈ Rm×n that belongs to the closed ball BεRsyl = {X : ∥X∥op ≤ εRsyl } which is the unique solution of the following Sylvester equation: X(Jxx + εUε Jyx ) − ε(Jyy − Jyx Uε )X = −εJyx 37

Proof. To make the notation lighter we drop ε from Uε such that U = Uε . We multiply the Sylvester −1 equation by Jxx on the right: −1 −1 −1 −1 X + εXU Jyx Jxx − εJyy XJxx + εJyx U XJxx = −εJyx Jxx

which gives the fixed point equation: −1 −1 −1 −1 X = ε Jyy XJxx − XU Jyx Jxx − Jyx U XJxx − Jyx Jxx | {z } T (X)

−1 We define the ball BεRsyl = {X : ∥X∥op ≤ εRsyl } where Rsyl = 2 ∥Jyx ∥op Jxx (This choice op U0 is explained in the next lines). Let M ∈ BεRsyl and U ∈ BR : ric −1 −1 −1 −1 ∥T (M )∥op = ε Jyy M Jxx − M U Jyx Jxx − Jyx U M Jxx − Jyx Jxx op −1 −1 −1 ≤ ε(∥Jyy ∥op ∥M ∥op Jxx + 2 ∥M ∥op ∥U ∥op ∥Jyx ∥op Jxx + ∥Jyx ∥op Jxx ) op op op −1 −1 −1 = ε(∥M ∥op (∥Jyy ∥op Jxx + 2 ∥U ∥op ∥Jyx ∥op Jxx ) + ∥Jyx ∥op Jxx ) op op op

≤ εRsyl

−1 −1 −1 εRsyl (∥Jyy ∥op Jxx + 2(Rric + ∥U0 ∥op ) ∥Jyx ∥op Jxx ) + ∥Jyx ∥op Jxx op op op

Rsyl

a first condition to set is that −1 −1 −1 ) + ∥Jyx ∥op Jxx + 2(Rric + ∥U0 ∥op ) ∥Jyx ∥op Jxx εRsyl (∥Jyy ∥op Jxx op op op

Rsyl

≤1

to make the operator T maps into BεRsyl which is equivalent to ε<

−1 Rsyl − ∥Jyx ∥op Jxx op −1 −1 Rsyl (∥Jyy ∥op Jxx + 2(Rric + ∥U0 ∥op ) ∥Jyx ∥op Jxx ) op op

−1 where we extract two conditions , the first is that Rsyl > ∥Jyx ∥op Jxx since ε is strictly positive, op −1 here we justify the choice of the radius Rsyl = 2 ∥Jyx ∥op Jxx , once we have done that we can op ∗ have a condition on ε ≤ ε2 to ensure that image of T stays in the ball BεRsyl such that:

ε∗2 =

1 ∥Jyy ∥op

−1 Jxx + 2(Rric + ∥U0 ∥op ) ∥Jyx ∥op op

−1 Jxx op

Now we move to the contraction part of T . Let X1 , X2 ∈ BεRsyl : −1 ∥T (X1 ) − T (X2 )∥op ≤ ε Jxx ∥Jyy ∥op + ∥U Jyx ∥op + ∥Jyx U ∥op ∥X1 − X2 ∥op op −1 ≤ ε Jxx ∥J ∥ + 2 ∥J ∥ (R + ∥U ∥ ) ∥X1 − X2 ∥op yy yx ric 0 op op op op The condition for the contraction is the same as the Riccati equation ε∗3 = ε∗1 then by Banach fixed point theorem, the map T has a fixed point solution X that is unique inside the ball BεRsyl for every ε < min{ε∗2 , ε∗3 }. In the next lemmas, we extract contractions of the bias dynamics using Lyapunov norms suitable for a given range of step sizes. The choice of the Lyapunov equations (80) and (83) appears in [22, Chapter 7] to study the stability of the equivalent singular perturbation system. The bound on the L-norm of I − γA is common and appears in [17, 15, 26, 21]. However, in our case we deal with nonlinear systems so we have extra perturbations emerging from the coupling between the two variables, so we extend the classical way to prove the contraction property under the suitable Lyapunov norms to account for the new added perturbations which are dependent on the solutions of the Riccati and the Sylvester equations that we showed previously to exist and to be bounded under certain assumptions on the step sizes. 38

Lemma 5. Let −Jxx be a Hurwitz matrix in Rd×d , Jyx ∈ Rn×d such that ∥Jyx ∥op ≤ MJ for some positive constant MJ , L ∈ Rd×d a symmetric positive definite matrix a solution to the following Lyapunov equation: T Jxx L + LJxx = I.

(80)

Let Rric > 0, there exists α∗ = α∗ (Rric ) where α∗ (Rric ) is defined by: α∗ (Rric ) =

(81)

2 2 4 ∥L∥op (∥Jxx ∥L + ε∗4 2 κ(L)(Rric + ∥U0 ∥op )2 ∥Jyx ∥L )

β such that for every 0 < α ≤ α∗ and for every β such that α < min{ε∗0 , ε∗1 , ε∗5 } where ε∗0 , ε∗1 , ε∗5 depend on Rric and are fully defined in (74), (75), (82), we define the matrix U β ∈ Rd×n to be α bounded by ∥U β ∥ ≤ Rric + ∥U0 ∥op , then we have: α

I − αJxx − βU β Jyx α

≤ e−cα L

where c = (8 ∥L∥op )−1 . Proof. Under the Lyapunov norm: 2

∥I − αJxx − βU Jyx ∥L = sup uT (I − αJxx − βU Jyx )T L(I − αJxx − βU Jyx )u ∥u∥L =1

T T = sup {uT (L − αJxx L − β(U Jyx )T L − αLJxx + α2 Jxx LJxx ∥u∥L =1

T T T + αβJyx U T LJxx − βLU Jyx + αβJxx LU Jyx + β 2 Jyx U T LU Jyx )u} T T = sup {uT Lu − αuT (Jxx L + LJxx )u + α2 uT Jxx LJxx u ∥u∥L =1

T T T − βuT (Jyx U T L + LU Jyx )u + αβuT (Jyx U T LJxx + Jxx LU Jyx )u {z } {z } | | A1

T + β 2 uT (Jyx U T LU Jyx )u} T = 1 + sup {−αuT u + α2 uT Jxx LJxx u − βuT A1 u + αβuT A2 u ∥u∥L =1

T + β 2 uT (Jyx U T LU Jyx )u} 2

≤ 1 + sup {−αuT u} + α2 ∥Jxx ∥L + β 2 ∥U Jyx ∥L + αβ L−1/2 A2 L−1/2 ∥u∥L =1

+ β sup −uT A1 u ∥u∥L =1 2

≤ 1 − α(2 ∥L∥op )−1 + α2 ∥Jxx ∥L + β 2 ∥U Jyx ∥L − α(2 ∥L∥op )−1 + αβ L−1/2 A2 L−1/2

+ β sup −uT A1 u ∥u∥L =1

−1

We have used the fact that inf ∥u∥L =1 {uT u} = ∥L∥op = (2 ∥L∥op )−1 + (2 ∥L∥op )−1 .We will use the symmetry of A1 to bound the last term such that β sup −uT A1 u ≤ β sup |uT A1 u| = β sup |z T L−1/2 A1 L−1/2 z| = β L−1/2 A1 L−1/2 ∥u∥L =1

∥u∥L =1

∥z∥=1

we have 2

∥I − αJxx − βU Jyx ∥L ≤ 1 − α(2 ∥L∥op )−1 + α2 (∥Jxx ∥L +

β2 2 ∥U Jyx ∥L ) α2

− α(2 ∥L∥op )−1 + βα L−1/2 A2 L−1/2

+ β L−1/2 A1 L−1/2

Since U is a solution of the Riccati equation (37) then by Lemma 3 for all Rric > 0 there exists β ε∗0 and ε∗1 defined in (74) and (75), respectively, such that for α < ε∗4 ≜ min{ε∗0 , ε∗1 } we have ∥U ∥ ≤ Rric + ∥U0 ∥op . We apply the result of lemma: 2

∥I − αJxx − βU Jyx ∥L ≤ 1 − α(2 ∥L∥op )−1 + α2 (∥Jxx ∥L + ε∗4 2 κ(L)(Rric + ∥U0 ∥op )2 ∥Jyx ∥L ) − α(2 ∥L∥op )−1 + βα L−1/2 A2 L−1/2 −1

≤ 1 − α(2 ∥L∥op )

+ β L−1/2 A1 L−1/2

op 2 2 ∗2 2 + α (∥Jxx ∥L + ε4 κ(L)(Rric + ∥U0 ∥op ) ∥Jyx ∥L ) 2

− α(2 ∥L∥op )−1 + 2βακ(L)(∥Jyx ∥op ∥Jxx ∥op (Rric + ∥U0 ∥op )) + 2βκ(L)1/2 (MJ (Rric + ∥U0 ∥op )) First, we look for values of α such that the first line is less than one to get the contraction property. This is true for the choice of α ≤ α∗ (Rric ) defined in the body of the lemma so we have: 2

∥I − αJxx − βU Jyx ∥L ≤ 1 − c′ α − α(2 ∥L∥op )−1 + 2βα∗ κ(L)(∥Jyx ∥op ∥Jxx ∥op (Rric + ∥U0 ∥op )) + 2βκ(L)1/2 (MJ (Rric + ∥U0 ∥op )) with c′ = (4 ∥L∥op )−1 . Second, we need the term in the last two lines to be negative so this will imply a condition of the type: β 1 ≤ ε∗5 ≜ α 2 ∥L∥op (Rric + ∥U0 ∥op )(2α∗ κ(L) ∥Jyx ∥op ∥Jxx ∥op + 2κ(L)1/2 MJ )

(82)

Once this is established we have: 2

∥I − αJxx − βU Jyx ∥L ≤ 1 − c′ α ≤ e−cα this implies that: c′

∥I − αJxx − βU Jyx ∥L ≤ e− 2 α which establishes the proof. Lemma 6. Let −∆ be a Hurwitz matrix in Rn×n , Jyx ∈ Rn×m such that ∥Jyx ∥op ≤ MJ for some positive constant MJ , G ∈ Rn×n a symmetric positive definite matrix a solution to the following Lyapunov equation: ∆T G + G∆ = I. Let Rric =

4∥G∥op

√1

κ(G)MJ

β∗ =

(83)

, we define β ∗ = β ∗ (Rric ) where β ∗ (Rric ) is defined by: 1 2

2 ) 4 ∥G∥op (2κ(G) ∥∆∥op MJ Rric + ∥∆∥G + κ(G)MJ2 Rric

(84)

β for all α > 0 and β ≤ β ∗ with α < min{ε∗0 , ε∗1 } where ε∗0 , ε∗1 are dependent on Rric and defined U0 in (74) and (75), the matrix U β ∈ BR is a solution of the Riccati equation defined in (76) with ric α

β −1 ε= α and U0 = Jxx Jxy , and we have:

I − β(∆ − Jyx (U β − U0 )) α

where d = (8 ∥G∥op )−1 . 40

≤ e−dβ G

Proof. Taking the Lyapunov norm 2

∥I − β(∆ − Jyx (U − U0 ))∥G = sup {uT (I − β∆ − βJyx (U − U0 ))T G(I − β∆ − βJyx (U − U0 ))u} ∥u∥G =1

= sup {uT Gu − βuT G∆u − βuT GJyx (U − U0 )u − βuT ∆T Gu ∥u∥G =1

+ β 2 uT ∆T G∆u + β 2 uT ∆T GJyx (U − U0 )u T T − βuT (U − U0 )T Jyx Gu + β 2 uT (U − U0 )T Jyx G∆u T + β 2 uT (U − U0 )T Jyx GJyx (U − U0 )u}

= 1 + sup {−βuT u − 2βuT GJyx (U − U0 )u + 2β 2 uT ∆T GJyx (U − U0 )u ∥u∥G =1 T + β 2 uT ∆T G∆u + β 2 uT (U − U0 )T Jyx GJyx (U − U0 )u} −1

≤ 1 − β ∥G∥op + 2β sup {−uT GJyx (U − U0 )u} ∥u∥G =1

+ 2β

sup {u ∆T GJyx (U − U0 )u} + β 2 ∥∆∥G + β 2 ∥Jyx (U − U0 )∥G

∥u∥G =1 −1

≤ 1 − β ∥G∥op + 2β G1/2 Jyx (U − U0 )G−1/2 + 2β

−1/2

∆ GJyx (U − U0 )G

op 2

+ β 2 ∥∆∥G + β 2 ∥Jyx (U − U0 )∥G

−1

≤ 1 − β ∥G∥op + 2βκ(G)1/2 MJ ∥U − U0 ∥op + 2β 2 G−1/2 ∆T GJyx (U − U0 )G−1/2

+ β 2 ∥∆∥G + β 2 ∥Jyx (U − U0 )∥G

Here we set the radius of convergence of the Riccati equation to calibrate for the negative term in the first line of the last inequality above. We want: 2κ(G)1/2 MJ ∥U − U0 ∥op ≤ (2 ∥G∥op )−1 which is equivalent to: ∥U − U0 ∥op ≤

1 = Rric 4 ∥G∥op κ(G)1/2 MJ

So under this condition we have 2

2 ∥I − β(∆ − Jyx (U − U0 ))∥G ≤ 1 − β(2 ∥G∥op )−1 + β 2 (2κ(G) ∥∆∥op MJ Rric + ∥∆∥G + κ(G)MJ2 Rric )

So the contraction is established for β ≤ β ∗ where β ∗ is defined in the body of the lemma and for d′ = (4 ∥G∥op )−1 we have ′

∥I − β(∆ − Jyx (U − U0 ))∥G ≤ 1 − dβ ≤ e−d β this implies that: d′

∥I − β(∆ − Jyx (U − U0 ))∥G ≤ e− 2 β

In the next lemma we handle the norms of the sum of deviations appearing in the proof of Theorem 3. Lemma 7. Taking the sequence {xk , yk }k≥0 defined by (1) and under assumptions A.1-2 and B.1. β Assume also that the step sizes α and β satisfy α < α∗ , β < β ∗ and α < ω ∗ such that α∗ , β ∗ , ω ∗ are defined in (8),(9),(10) then, the term (48) is bounded by: (48) ≤αA1 + β(A2 + B1 ) + α2 A3 + β 2 B2 + αβ(A4 + B3 ) +

β2 B4 α

in addition the term (52) is bounded by: (52) ≤α(C2 + D4 ) + β(C1 + D1 ) + αβ(C4 + D3 ) + β 2 (C3 + D2 ) where the constants are independent of α and β and are detailed in the proof. 41

Proof. We start with bounding the first term: α

k X

α,β k−j (I − αJxx ) E [ψj ]

(85)

j=0

In Lemma 2 we have seen that there exists a function vf the solution of the Poisson equation (54), which allowed us to write: E [ψs ] = E vf (xs , ys , ξs+1 ) − Eξ′ ∼Ks (ξs+1 ,·) [vf (xs , ys , ξ ′ )] we could write the term inside the expectation as: ψs =vf (xs , ys , ξs+2 ) − Eξ′ ∼Ks+1 (ξs+1 ,·) [vf (xs , ys , ξ ′ )]

(86)

+ vf (xs+1 , ys+1 , ξs+2 ) − vf (xs , ys , ξs+2 ) + vf (xs , ys , ξs+1 ) − vf (xs+1 , ys+1 , ξs+2 )

(87) (88)

+ Eξ′ ∼Ks+1 (ξs+1 ,·) [vf (xs , ys , ξ ′ )] − Eξ′ ∼Ks (ξs+1 ,·) [vf (xs , ys , ξ ′ )]

(89)

where the first term (86) is a Martingale Difference Sequence and we denote Ks ≜ K(xs ,ys ) . By taking the expectation, we get: E[vf (xs , ys , ξs+2 ) − Eξ′ ∼Ks+1 (ξs+1 ,·) [vf (xs , ys , ξ ′ )]] =E Eξ′ ∼Ks+1 (ξs+1 ,·) [vf (xs , ys , ξ ′ )] − Eξ′ ∼Ks+1 (ξs+1 ,·) [vf (xs , ys , ξ ′ )] =0 moving to the second term (87) and for legibility, we set vf (xs , ys , ξs+2 ) = vs and vf (xs+1 , ys+1 , ξs+2 ) = vs+1 we have: α

k X

α,β k−s (I − αJxx ) (vs − vs+1 ) ≤α

s=0

k X

α,β I − αJxx

k−s

∥vs − vs+1 ∥

s=0

≤α

κ(L)

k X

e−cα(k−s) ∥vs − vs+1 ∥

s=0

κ(L)Lvf

k X

e−cα(k−s) (∥xs − xs+1 ∥ + ∥ys − ys+1 ∥)

s=0

κ(L)Lvf (αMf + βMg )

k X

e−cαs

s=0

Z ∞ κ(L)Lvf (αMf + βMg )(1 + e−cαt dt) 0 p p =αc−1 κ(L)Lvf Mf + βc−1 κ(L)Lvf Mg p p + α2 κ(L)Lvf Mf + αβ κ(L)Lvf Mg

≤α

Next, we handle the pseudo-telescoping term (88), but before that we set vf (xs , ys , ξs+1 ) = vs , β α,β α,β vf (xs+1 , ys+1 , ξs+2 ) = vs+1 , Jxx = Jxx + α U Jyx and Aα,β = I − αJxx . From the Neumann series convergence, Aα,β is invertible if we have a condition of the type α,β α Jxx <1 op

since we are assuming α<

1 ∥Jxx ∥op + (Rric + ∥U0 ∥op ) ∥Jyx ∥op

(90)

then the condition is satisfied and Aα,β is invertible and the inverse can be given by the Neumann P∞ α,β k series A−1 . However, we will not need the explicit form of the inverse, so we k=0 αJxx α,β = 42

use it directly: k X

Ak−s α,β (vs − vs+1 ) =

s=0

k X

k−(s+1)

Ak−s α,β vs − Aα,β Aα,β

vs+1

s=0

k X

k−(s+1)

Ak−s α,β vs − Aα,β

α,β −1 vs+1 + αJxx Aα,β

s=0

k X

Ak−s α,β vs+1

s=0

α,β −1 =Akα,β v0 − A−1 α,β vk+1 + αJxx Aα,β

k X

Ak−s α,β vs+1

s=0

taking the norm: α

k X

2 α,β A−1 Jxx Ak−s α,β α,β (vs − vs+1 ) ≤α op

s=0

+α ≤α2

A−1 α,β

k X op

Ak−s α,β

s=0

∥vs+1 ∥ + α Akα,β op ∥v0 ∥

∥vk+1 ∥

α,β κ(L) Jxx A−1 α,β op

Mξ + α A−1 α,β op p α,β ≤ κ(L) Jxx A−1 α,β op

Mξ

k X

e−cαs + αe−cαk

κ(L)Mξ

s=0

p α Mξ α2 + + αe−cαk κ(L)Mξ c op

+ α A−1 Mξ α,β op p p α,β A−1 Mξ + αe−cαk κ(L)Mξ =α2 κ(L) Jxx α,β op op p α,β −1 κ(L) Jxx op A−1 Mξ + α A−1 Mξ + αc α,β α,β op op p where in the second inequality we used ∥D∥op ≤ κ(L) ∥D∥L . The last term (89) is capturing the drift in the kernel and is bounded in the same way as (58) by Eξ′ ∼Ks+1 (ξs+1 ,·) [vf (xs , ys , ξ ′ )] − Eξ′ ∼Ks (ξs+1 ,·) [vf (xs , ys , ξ ′ )] ≤ Mξ MS LK |Ξ|(αMf + βMg ) which yields the bound α

k X p α,β k−s (I − αJxx ) (89) ≤ κ(L)MS Mξ LK |Ξ|c−1 (αMf + βMg ) s=0

We pack all the bounds together, under our assumptions on the step sizes there is CA , CJ that satisfy α,β A−1 ≤ CA and Jxx ≤ CJ so we have: α,β op op p p −1 (85) ≤α κ(L) c−1 Mf (Lvf + MS Mξ LK |Ξ|) + Mξ 1 + CA c−1 CJ + κ(L) {z } | ≜A1

−1

+β c |

p κ(L)Mg (Lvf + MS Mξ LK |Ξ|) +α2 κ(L) Lvf Mf + CJ CA Mξ | {z } {z } ≜A3

≜A2

p + αβ κ(L)Lvf Mg | {z } ≜A4

The second part of (48) concerns bounding the term: β

k X α,β k−j (I − αJxx ) U E [φj ] j=0

(91)

that is similar to the first one, except that now we have a multiplication by β and the deviations concern the dynamics of the slow variable. There exists a function vg solution of the Poisson equation: vg (xj , yj , ξj+1 ) − Eξ′ ∼Kj (ξj+1 ,·) [vg (xj , yj , ξ ′ )] = φj so we can write: φj =vg (xj , yj , ξj+2 ) − Eξ′ ∼Kj+1 (ξj+1 ,·) [vg (xj , yj , ξ ′ )]

(92)

+ vg (xj+1 , yj+1 , ξj+2 ) − vg (xj , yj , ξj+2 ) + vg (xj , yj , ξj+1 ) − vg (xj+1 , yj+1 , ξj+2 )

(93) (94)

+ Eξ′ ∼Kj+1 (ξj+1 ,·) [vg (xj , yj , ξ ′ )] − Eξ′ ∼Kj (ξj+1 ,·) [vg (xj , yj , ξ ′ )]

(95)

where term (92) is an MDS, under expectation it will equal zero. The second term (93), similarly to the previous part, we set vg (xj , yj , ξj+2 ) = vj and vg (xj+1 , yj+1 , ξj+2 ) = vj+1 we have

k X

α,β k−j (I − αJxx ) U (vj − vj+1 ) ≤β 2

κ(L) ∥U ∥op Lvg Mg + αβ

κ(L) ∥U ∥op Lvg Mf

j=0

p β2 p κ(L) ∥U ∥op Lvg Mg + βc−1 κ(L) ∥U ∥op Lvg Mf cα for the third term (94) it will follow the same bound as the one used in (88), we just need to multiply it by β, we set vg (xj , yj , ξj+1 ) = vj , vg (xj+1 , yj+1 , ξj+2 ) = vj+1 : +

k X

Ak−j α,β U (vj − vj+1 ) ≤αβ

p α,β κ(L) Jxx ∥U ∥op A−1 α,β op

j=0

+ βc−1

Mξ + βe−cαk

p α,β κ(L) Jxx ∥U ∥op A−1 α,β op

κ(L) ∥U ∥op Mξ

Mξ + β A−1 α,β

∥U ∥op Mξ

Finally, we bound the term β

k X p β α,β k−s (I − αJxx ) U (95) ≤ ∥U ∥op κ(L)MS Mξ LK |Ξ|c−1 (αMf + βMg ) α s=0

Putting all the bounds together, we get: p p −1 (91) ≤β ∥U ∥op κ(L) c−1 (Lvg + MS Mξ LK |Ξ|)Mf + Mξ CA c−1 CJ + κ(L) +1 {z } | ≜B1

+β

κ(L) ∥U ∥op Lvg Mg +αβ κ(L) CJ ∥U ∥op CA Mξ + ∥U ∥op Lvg Mf {z } | {z } |

p 2

≜B2

≜B3

β2 p + κ(L) ∥U ∥op Mg c−1 (Lvg + MS Mξ LK |Ξ|) α | {z } ≜B4

We merge the major two bounds together to get: β2 B4 α ending the first part of the proof. Now we start with the second part to bound (52). We start by the first term: (85) + (91) ≤αA1 + β(A2 + B1 ) + α2 A3 + β 2 B2 + αβ(A4 + B3 ) +

k X

α,β k−j (I − βJyy ) XE [ψj ]

(96)

j=0

as we have discussed above about rewriting the inside terms with the solutions of the Poisson equation, we set vf (xs , ys , ξs+2 ) = vs and vf (xs+1 , ys+1 , ξs+2 ) = vs+1 we have: 44

k k X X α,β k−s α,β k−s ∥X∥ ∥vs − vs+1 ∥ α (I − βJyy ) X(vs − vs+1 ) ≤α I − βJyy s=0

s=0

≤α

κ(G) ∥X∥

k X

e−dβ(k−s) ∥vs − vs+1 ∥

s=0

1 ) dβ p 1 ) ≤βRsyl κ(G)Lvf (αMf + βMg )(1 + dβ p p =αβRsyl κ(G)Lvf Mf + β 2 Rsyl κ(G)Lvf Mg p p + αd−1 Rsyl κ(G)Lvf Mf + βd−1 Rsyl κ(G)Lvf Mg

≤α

κ(G)Lvf ∥X∥ (αMf + βMg )(1 +

where in the fourth inequality we have used the bound (41), now we target the almost telescoping α,β term. We set vf (xs , ys , ξs+1 ) = vs , vf (xs+1 , ys+1 , ξs+2 ) = vs+1 , Jyy = ∆ − Jyx (U − U0 ) and α,β Bα,β = I − βJyy . Since we are assuming

β<

1 ∥∆∥op + Rric ∥Jyx ∥op

(97)

then Bα,β is invertible using the Neumann series, we use the inverse to have: k X

k−s Bα,β X(vs − vs+1 ) =

s=0

k X

k−(s+1)

k−s Bα,β Xvs − Bα,β Bα,β

Xvs+1

s=0

k X

k−(s+1)

k−s Bα,β Xvs − Bα,β

α,β −1 Xvs+1 + βJyy Bα,β

k X

k−s Bα,β Xvs+1

s=0

s=0 −1 k α,β −1 =Bα,β Xv0 − Bα,β Xvk+1 + βJyy Bα,β

k X

k−s Bα,β Xvs+1

s=0

applying the norm:

k X

k−s −1 k Bα,β X(vs − vs+1 ) ≤α Bα,β ∥X∥op ∥v0 ∥ + α Bα,β op

s=0 −1 α,β + αβ Jyy Bα,β op

≤β

κ(G)Rsyl e

+ β2

−dβk

∥X∥op

k X

∥X∥op ∥vk+1 ∥ k−s

∥Bα,β ∥op ∥vs+1 ∥

s=0

p −1 Mξ + β κ(G)Rsyl Bα,β Mξ

α,β κ(G)Rsyl Jyy

−1 Bα,β Mξ

k X

e−dβ(k−s)

s=0

p −1 α,β ≤βRsyl κ(G)Mξ e−dβk + Bα,β + Jyy p −1 α,β β 2 Rsyl κ(G)Mξ Jyy Bα,β

−1 Bα,β d−1

and we get our bound for the first part, under our assumptions on the step sizes there is CA , CJ that −1 α,β ≤ CJ so we have: ≤ CB and Jyy satisfy Bα,β op p α,β (96) ≤β Rsyl κ(G) d−1 (Lvf + MS Mξ LK |Ξ|)Mg + Mξ 1 + CB + Jyy CB d−1 | {z } ≜C1

+ α Rsyl |

κ(G)(Lvf + MS Mξ LK |Ξ|)Mf d−1 +β 2 Rsyl {z } |

≜C2

+ αβ Rsyl |

κ(G) Mξ CJ CB + Lvf Mg {z } ≜C3

κ(G)Lvf Mf {z } ≜C4

our last bound is for the term : β

k X

α,β k−j ′ (I − βJyy ) I E [φj ]

(98)

j=0

from the Poisson decomposition the first term will equal zero under expectation. For the second term (93), similarly to the previous part, we set vg (xj , yj , ξj+2 ) = vj and vg (xj+1 , yj+1 , ξj+2 ) = vj+1 we have: β

k X

α,β k−j ′ (I − βJyy ) I (vj − vj+1 ) ≤αβ ∥I ′ ∥op

κ(G)Lvg Mf + βd−1 ∥I ′ ∥op

κ(G)Lvg Mg

j=0

+ αd−1 ∥I ′ ∥op

κ(G)Lvg Mf + β 2 ∥I ′ ∥op

κ(G)Lvg Mg

and the last term, we set vj = vg (xj , yj , ξj+1 ), vj+1 = vg (xj+1 , yj+1 , ξj+2 ): k X p p −1 −1 k−s ′ α,β Bα,β κ(G)e−dβk + Bα,β + κ(G) Jyy β Bα,β I (vs − vs+1 ) ≤β ∥I ′ ∥op Mξ op op

s=0

+ β 2 ∥I ′ ∥op Mξ

−1 α,β κ(G) Jyy Bα,β op

which implies the bound on the term: p (98) ≤β ∥I ′ ∥op Mξ CB + κ(G) 1 + d−1 CJ CB + ∥I ′ ∥op (Lvg + MS Mξ LK |Ξ|)Mg {z } | ≜D1

+β

′

∥I ∥op Mξ

| + αβ ∥I ′ ∥op |

′

κ(G)CJ CB + ∥I ∥op {z

κ(G)Lvg Mg

}

≜D2

κ(G)Lvg Mf +α d−1 ∥I ′ ∥op {z } |

≜D3

κ(G)(Lvg + MS Mξ LK |Ξ|)Mf {z } ≜D4

which implies the following bound: (96) + (98) ≤ α(C2 + D4 ) + β(C1 + D1 ) + αβ(C4 + D3 ) + β 2 (C3 + D2 ) which concludes the proof. Lastly, this lemma is needed to bound the weighted sum of the MSEs of the fast and slow iterates. We need this result in the proof of Theorem 3. Lemma 8. Taking the sequence {xk , yk }k≥0 defined by (1) and under assumptions A.1-4, let q, γ be positive real numbers, we have: k h i X 1 2 2 q,γ −qγ(k−s) e E ∥ŷs ∥ + ∥x̂s ∥ ≤ 1 + (MSEx∞ (α, β) + MSEy∞ (α, β)) + uq,γ k + wk qγ s=0 46

d−1 op

where MSEx∞ (α, β) and MSEy∞ (α, β) are defined in Theorem 1 and 2

−kqγ uq,γ k = ∥x̂0 ∥ e

k X

(eqγ (1 − αµf ))s

s=0

wkq,γ =

k X

e−qγ(k−s) (1 − βµg ) ∥ŷ0 ∥ +

k X

s=0

e−qγ(k−s) D1 βs(1 − αµf ∧ βµg )s ∥x̂0 ∥

s=0

are vanishing sequences as k → ∞. Proof. First, we start by analyzing the fast iterate, we invoke the result of Theorem 1 to bound its MSE, we start by the sum of MSEs of the fast iterate: k X

k h i X 2 2 e−qγ(k−s) E ∥x̂s ∥ ≤ e−qγ(k−s) (1 − αµf )s ∥x̂0 ∥ + MSEx∞ (α, β)

s=0

≤

k X

−qγ(k−s)

(1 − αµf ) ∥x̂0 ∥

+ MSEx∞ (α, β)

Z ∞ −qγt 1+ e dt 0

s=0 k X

1 2 = e−qγ(k−s) (1 − αµf )s ∥x̂0 ∥ + MSEx∞ (α, β) 1 + qγ s=0 2

≤ ∥x̂0 ∥ e−kqγ

k X

(eqγ (1 − αµf ))

(99)

s=0

1 + 1+ qγ

MSEx∞ (α, β)

Now from our assumptions we have 1 − αµf ∈ (0, 1]. If eqγ (1 − αµf ) = 1 then the sum equals Pk e−kqγ s=0 (eqγ (1 − αµf ))s = e−kqγ (k + 1) which will go to 0 as k → ∞. For the other case eqγ (1 − αµf ) ̸= 1 so the sum : −kqγ

k X e−kqγ − eqγ (1 − αµf )k+1 (eqγ (1 − αµf ))s = 1 − eqγ (1 − αµf ) s=0

which will also disappear as k goes to infinity. So we set the sequence uq,γ k to be defined by: 2

−kqγ uq,γ k = ∥x̂0 ∥ e

k X

(eqγ (1 − αµf ))

s=0

such that

lim uq,γ k =0

k→∞

Same reasoning goes with the MSE of the slow iterate, from Theorem 1 following the same steps in the fast MSE above we get: k X

k h i X 2 s 2 e−qγ(k−s) E ∥ŷs ∥ ≤ e−qγ(k−s) (1 − βµg ) ∥ŷ0 ∥

s=0

k X

e−qγ(k−s) D1 βs(1 − αµf ∧ βµg )s ∥x̂0 ∥

(100)

e−qγ(k−s) MSEy∞ (α, β)

(101)

s=0

k X s=0

The first term follows the same reasoning as the term (99), we focus on the second term (100): 2

(100) ≤ ∥x̂0 ∥ D1 βke−qγk

k X s=0

(eqγ (1 − αµf ∧ βµg ))

if eqγ (1 − αµf ∧ βµg ) = 1: ke−qγk

k X

(eqγ (1 − αµf ∧ βµg )) = e−qγk k(k + 1)

s=0

which vanishes as k → ∞. otherwise, if eqγ (1 − αµf ∧ βµg ) ̸= 1: ke−qγk

k X

ke−qγk − keqγ (1 − βµg ∧ αµf )k+1 1 − eqγ (1 − βµg ∧ αµf )

(eqγ (1 − αµf ∧ βµg )) =

s=0

which also vanishes as k → ∞. We define the new sequence wkq,γ as: wkq,γ =

k X

e−qγ(k−s) (1 − βµg ) ∥ŷ0 ∥ +

s=0

k X

e−qγ(k−s) D1 βs(1 − αµf ∧ βµg )s ∥x̂0 ∥

s=0

such that

lim wkq,γ = 0

k→∞

and the last term (101) is bounded by (101) ≤

1 1+ qγ

MSEy∞ (α, β)

now we can write the final bound: k h i X 1 2 2 q,γ −qγ(k−s) e E ∥ŷs ∥ + ∥x̂s ∥ ≤ 1 + (MSEx∞ (α, β) + MSEy∞ (α, β)) + uq,γ k + wk qγ s=0

Illustrations

F.1

First example

We investigate the empirical performance of a TTSA where both variables are 1-dimensional: xk ∈ R and yk ∈ R and satisfy the following recursion x xk+1 = xk − α(xk − yk + ξk+1 ) (102) y −yk yk+1 = yk − β(2yk − e − xk + ξk+1 ), where the noise terms ξkx , ξky ∈ {−1, 1} form the Markovian noise with transition probability kernel P[ξk+1 = s′ | ξk = s] = 0.3 for s ̸= s′ . Initially, we set x0 = y0 = 0. For this example, h(y) = y and y ∗ is the solution of y = exp(−y) which is equal to W0 (1) ≈ 0.567, where W0 is the standard Lambert function. In Figure 2, we display a trajectory of the TTSA for three values of α ∈ {0.1, 0.01, 0.002} and a fixed β = 0.001. We observe two phenomena that are shown in Theorem 1: 1.5

0.8

1.0

0.6

0.5

0.4

0.0

Fast iterate xk Slow iterate yk y*

0.5 0

2000 4000 6000 8000 10000 Iterate k

(a) α = 0.1, β = 0.001

Fast iterate xk Slow iterate yk y*

0.2 0.0 0

2000 4000 6000 8000 10000 Iterate k

(b) α = 0.01, β = 0.001

0.6 0.5 0.4 0.3 0.2 0.1 0.0

Fast iterate xk Slow iterate yk y* 0

2000 4000 6000 8000 10000 Iterate k

Figure 2: Sample trajectory of the TTSA defined in (102). 48

1. First, the vanishing term vanishes at a slower rate as α gets smaller. For instance, when α = 0.002 (Figure 2(c)), the fast variable is lagging behind h(yk ) = yk where when α = 0.01 (Figure 2(b)), the variable xk has no lag yk . 2. Second, the larger is α, the larger are the oscillations of xk and yk : in Figure 2(a), the variable xk is not lagging behind h(yk ) = yk but it is oscillating a lot. In Figure 3, we focus on the asymptotic behavior of the MSE and of the bias as a function of the step size α and for the case where there is a time scale separation β = α1.5 . For each value, we simulate a trajectory of length k = 1.2 · 105 steps and compute the average values of xk for the last half of the trajectory. We perform 5000 independent trajectories to obtain the confidence intervals (they are shown on both panels but are too small to be visible on the MSE plot in Figure 3(a)). This figure shows that for this example, MSEx seems to be of the order of α whereas MSEy and the bias terms are of the order of β = α1.5 .

0.10 0.08 0.06

Curve Curve = 1.5 MSE x (simulation) MSE y (simulation)

0.04 0.02 0.00

0.02

0.04

0.06

Curve 0.2 1.5 Biasx (simulation) Biasy (simulation)

0.006 0.005 0.004 0.003 0.002 0.001 0.000

0.08

0.10

0.02

(a) MSE as a function of α (for β = α3/2 )

0.04

0.06

0.08

0.10

(b) Bias as a function of α (for β = α3/2 )

Figure 3: MSE and Bias of the TTSA defined in (102) for β = α3/2 F.2

Example showing a lower bound in Ω(β 2 /α2 )

The example used in the paper (Figure 1(c)) to illustrate the Ω(β 2 /α2 ) term in the MSE is a TTSA whose recurrence equation is given by xk+1 = xk − α(xk − yk ) (xk − yk ) yk+1 = yk − β , |xk − yk |γ

(103)

with x0 = 0 and y0 = 1. In this example, xk chases the variable yk (h(yk ) = yk ), while yk tries to escape xk . This example can be easily analyzed: Rewriting δk ≜ yk − xk , one obtains δk+1 = (1 − α)δk + β

(δk ) |δk |γ

As k goes to infinity, this shows that limk→∞ δk = (β/α)1/γ . This shows that: 2/γ β x lim MSE = . k→∞ α When γ is close to 1, this is essentially equal to (β/α)2 . This is what is illustrated in Figure 4: for this example, the quantity MSEx does not depend on α and is equal to (β/α)2/γ for all α, β and γ. Discussion on the assumptions The TTSA (103) satisfies the assumptions about the kernel and β2 the fast iterate, since here we care about the tracking error of the fast iterate that is of order α 2 and 49

0.6 0.4 0.2 0.0

1.0

Curve ( / )2.50 = 0.01 = 0.5

0.8 MSEx

MSEx

1.0

Curve ( / )4.00 = 0.01 = 0.5

0.8

0.6 0.4 0.2

0.2

0.4 0.6 Ratio /

0.8

1.0

0.0

Curve ( / )2.11 = 0.01 = 0.5

0.8 MSEx

1.0

0.6 0.4 0.2

0.2

(a) γ = 0.5

0.4 0.6 Ratio /

0.8

1.0

0.0

0.2

(b) γ = 0.8

0.4 0.6 Ratio /

0.8

1.0

Figure 4: MSE of the TTSA given by (103) as a function of the ratio β/α. We observe that MSEx does not depend on α and is equal to (β/α)2/γ for all α, β and γ. x

we got it from the MSE analysis of x without using information about Lipschitzness nor the strong monotonicity of g so A.4 was not needed in that analysis. In addition to that we relax the assumption on the boundedness of the iterates as we could not find an example that guarantees boundedness. The fact that xk and yk do not live in a compact set is more problematic but can be solved by slightly modifying the above example. Note that this new example cannot be easily solved in closed form which is why we decided to keep the TTSA (103) as the main example. To construct this new example, we consider 2-dimensional variables xk , yk ∈ R2 with x0 = [1, 0], y0 = [1, .1] and xk+1 = xk − α(xk − yk + 0.1ξk+1 ) (xk − yk ) 2 + y tanh(∥y ∥ − 1) , yk+1 = yk − β k k |xk − yk |γ

(104)

where ξk+1 is a sequence of Rademacher random variables. This example is essentially the same as the TTSA in (103) with two differences. First, we add a small white noise in the dynamics of xk (this does not change anything, but shows that the counter-example does not depend on the dynamics being deterministic). Second, and more importantly, the evolution 2 of yk contains the term yk tanh(∥yk ∥ − 1). The role of this term is to push yk towards the unit circle. As an example, we display sample trajectories of (xk , yk ) in Figure 5 for various values of α and β = 0.1α. We observe that both iterates xk and yk oscillate on the circle (with more noise when α is larger). 1.0 0.5 x y

0.0

1.0

0.5

0.5 x y

0.0

0.5

1.0

0.5

1.0 1.0

0.5

0.0

0.5

(a) α = 0.1

1.0

x y

0.0

1.0 1.0

0.5

0.0

0.5

(b) α = 0.01,

1.0

0.5

0.0

0.5

1.0

Figure 5: Sample trajectories of xk and yk for the TTSA (104) for various values of α and β = α/10 For this dynamics, we again have h(yk ) = yk . In Figure 6(a), we plot a sample trajectory of the tracking error xk − h(yk ) = xk − yk for various values of α while keeping the ratio β/α = 0.1 fixed. We observe that, as α goes to 0, the error does not vanish but seems to converge to a circular dynamics. Similarly to what was shown in Figure 4 for the TTSA (103), we observe in Figure 6(b) that the MSE of x grows like (β/α)2/γ , independently of α. Note that this model does not satisfy the strong monotonicity assumption for ḡ (A.4). It is an open question if one can construct an example that satisfies A.4, and whose asymptotic MSE for x has a dependence on (β/α). 50

0.15

0.005

0.00

0.003

0.05

0.002

0.10

0.001

0.15

= 0.01 = 0.001 Curve ( / )2/

0.004

0.05 MSEx

xk, 1 yk, 1

0.10

0.006

alpha=0.1 alpha=0.01 alpha=0.001

0.000

0.1

0.0 xk, 0 yk, 0

0.1

(a) Error xk − h(yk ) for β = α/10

0.02

0.04 0.06 Ratio /

0.08

(b) MSE as a function of β/α

Figure 6: Error and MSE of the TTSA (104).

0.10

Record · ID 978447 · SHA-256 13820a0b46c6e08a
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.