Stochastic Approximation (SA) is an iterative procedure to find the root of unknown operators from their noisy samples [39]. There has been a long line of work understanding the convergence behavior of SA both asymptotically [4,21] and in a finite-time [40], with a wide range of applications including stochastic optimization [21,37] and reinforcement learning [42,32,40].
Two-Timescale Stochastic Approximation (TTSA) is a variant of the SA algorithm, designed to find the root of two intertwined operators [3]. In particular, given two operators F and G, we aim to find the solution (x * , y * ) satisfying the fixed-point equations F (x * , y * ) = 0, G(x * , y * ) = 0. This work considers linear TTSA with constant stepsizes driven by Markovian data as the following:
x t+1 = x t -α t (F (x t , y t ) + w x (x t , y t ; ξ t )), y t+1 = y t -β t (G(x t , y t ) + w y (x t , y t ; ξ t )), t ≥ 0,
where α t ≡ α, β t ≡ β > 0 are constant stepsizes for slower and faster iterates respectively, F and G are linear operators, and w x and w y are linear Markovian noises driven by exogenous Markovian states ξ t (see Section 2 for precise formulation). The updates in (1) arise in many applications: examples include popular reinforcement learning algorithms such as actor-critic [29,18] and gradient temporal-difference (GTD) methods [35,43], and iterative algorithms for stochastic Bilevel optimization [7,16,22,31]. The core idea of TTSA is the use of different stepsizes for two iterations, which establishes a hierarchy between the two fixed-point equations. For example, in actor-critic algorithms [18], the y-variable minimizes the temporal-difference (TD) error, while the x-variable represents policy parameters to maximize long-term Existing Results for TTSA. TTSA arises as a popular iterative solution in various domain; from the classical iterateaveraging schemes [38] and off-policy reinforcement learning algorithms [42] to gradient descent-ascent algorithms for saddle-point problems [27] and single-loop algorithms for Bilevel optimization [22]. Asymptotic convergence and central limit theorems for TTSA with diminishing step sizes were initially established for linear cases with i.i.d. noise [30], followed by extensions to non-linear and Markovian noise settings [36,23].
More recent work has shifted focus to non-asymptotic results, deriving finite-time convergence rates for both linear [9,8] and nonlinear cases [28,19,20]. However, these studies primarily address MSE bounds with diminishing stepsizes. In contrast, we investigate distributional convergence under constant stepsizes, providing explicit decoupling of biases and variances. Additionally, we establish new results for tail-averaging and extrapolation in TTSA schemes.
Let F : R dx × R dy → R dx and G : R dx × R dy → R dy be linear mean-field operators in the following form:
where J 11 , . . . , J 22 (resp., b 1 , b 2 ) are fixed matrices (resp., vectors), and linear Markovian noises defined as the following: w x (x, y; ξ) = W 11 (ξ)x + W 12 (ξ)y + u 1 (ξ), w y (x, y; ξ) = W 21 (ξ)x + W 22 (ξ)y + u 2 (ξ).
Let J max := max i,j∈{1,2} ∥J ij ∥ op be the smoothness parameter of the system. The first assumption is on the mean-field operators being Hurwitz: Assumption 1 The matrices -J 22 and -∆ := -J 11 + J 12 J -1 22 J 21 are Hurwitz, that is, all real parts of the eigenvalues of J 22 and ∆ are strictly positive.
Therefore, a fixed point in the slower timescale is uniquely defined y * (x) = -J -1 22 (J 21 x + b 2 ) for every x, and the target joint fixed point (x * , y * ) is given as:
Furthermore, for all ξ ∈ Ξ, the following holds:
for some constants W max , u max ≥ 0. For simplicity, we further assume that W max ≤ J max .
The above two assumptions are common in the analysis of SA schemes with Markovian noises [8,24]. We introduce the notion of noise variances in our setting:
which reflect the mean-squared fluctuation of the stochastic update around the fixed point.
We study the convergence of TTSA iterations (1) via L 2 -Wasserstein distance [44]. Let P 2 (R d ) denote the space of square-integrable distributions on R d where d := d x + d y . Note that L 2 -Wasserstein distance between two distributions µ and ν in P 2 (R d ) is defined as the following:
, where Π(µ, ν) is a set of all possible couplings with marginal distributions µ and ν. To study the distribution convergence of the joint iterate-data sequence (x t , y t , ξ t ) t≥0 , we slightly extend the definition above to add hamming distance in Ξ.
Let P 2 (R d ×Ξ) be the set of distributions μ on R d ×Ξ with the property that the marginal of μ on R d is square-integrable.
Definition 1 For any two probability measures µ, ν in P 2 (R dx+dy × Ξ) over (x, y, ξ), we define the distance between µ and ν as
where Π(µ, ν) is a set of all possible couplings with marginal distributions µ, ν.
To establish the finite-time convergence of TTSA iterations (1), we define a few error metrics. Let Q x , Q y ≻ 0 be the unique solutions of the Lyapunov equations
The solutions Q x , Q y , which are guaranteed to exist since -∆, -J 22 are Hurwitz under Assumption 1 [5], are used for constructing the drift of potentials in our analysis. For the slower and faster iterates, we use ∥ • ∥ Qx and ∥ • ∥ Qy norms respectively, and define µ x := ∥Q x ∥ -1 op and µ y := ∥Q y ∥ -1 op . Note that σ min (∆) ≥ µ x /2 and σ min (Q x ) ≥ ∥∆∥ -1 op /2, and similarly for Q y and µ y . Consequently, we let the condition number of two iterations as κ x := κyJmax µx and κ y := Jmax µy .
Notation. For a positive definite matrix Q ≻ 0 let ∥x∥ Q := x ⊤ Qx for a vector x. With a general real-valued matrix A, we define
For two real-valued matrices A, B, we denote ⟨A, B⟩ = Tr(A ⊤ B). We define 1-Schatten norm ∥A∥ 1 := i |σ i (A)| as the absolute sum of singular values (sometimes we call it S 1 -norm), and ∞-Schatten norm ∥A∥ ∞ := max i |σ i (A)| be the maximum singular value, which is equivalent to matrix operator norm ∥A∥ op . For a positive semi-definite matrix Q ⪰ 0,
ii is the sum of diagonal elements. For a random vector x, we denote the covariance
We often use shorthands w x t := w x (x t , y t ; ξ t ) and w y t := w y (x t , y t ; ξ t ). We denote the fixed point of the faster iterates given x as y * (x), such that G(x, y * (x)) = 0. If we just write y * , then it means y * (x * ). For two probability distributions p, q, we denote ∥p -q∥ 1 as the total-variation distance between p and q. We use the notation O(•) to hide absolute constants, and O P (•) to omit up to polynomial factors in instance-dependent constants (smoothness, minimum eigenvalues, and noise variances) and up to logarithmic factors in stepsizes.
We start with two conditions for stepsizes to ensure the stability of TTSA iterations: Assumption 4 We assume that the stepsizes (α, β) satisfy the following:
where τ α := log(αµx/cρ) log ρ with some sufficiently small absolute constants c 1 , c 2 > 0.
The first condition in (4) ensures β less than the inverse smoothness of operators, and the second condition bounds the ratio between two-timescale iterations. We mention that the dependence on the condition numbers is not fully optimized. In the sequel, we start with a fine-grained convergence in MSE in Section 3.1. We then show the convergence in distribution and characterize the biases and variances of the limit distribution in 3.2, which is followed by our final result on tail-averaging and extrapolation in Section 3.3.
We analyze the MSE convergence of linear TTSA in terms of the centered iterates xt := x t -x * , ȳt := y t -y * (x t ).
To this end, we first rewrite the stochastic recursion as the following:
Then equation (1) can be rewritten as:
Note that the slower iterates view the error in faster iterates as an additional noise. We are now ready to state our first main convergence theorem with constant step-sizes.
Theorem 3.2 Suppose Assumptions 1-3 hold, and the step sizes α, β satisfy Assumption 4. Then, for all t ≥ 0 following the TTSA recursion (1), we have
where we define potential functions as
The theorem states that after sufficiently large iterations t ≫ α -1 , the convergence of TTSA in MSE can be characterized as the following:
To our best knowledge, this is the first result that explicitly characterizes the fine-grained scaling of MSE w.r.t. the stepsizes and noise variances of each iteration. The work in [9,40] only obtained an O(β 2 /α) asymptotic bound for the slower iterate. More recent work in [28,20] obtained an O(α) asymptotic bound but required β 2 ≤ α, hence not strong enough to reveal the dependence on β. Our result shows that noises from slower iterates only change x t by O(α), while noises from faster iterates influence x t by O(α + β 2 ), without requiring β 2 ≤ α.
Now we state the distributional convergence of the process (x t , y t , ξ t ) in Wasserstein distance as defined in Definition 1.
We require a mild assumption on the fourth-order moments of initial distributions:
Assumption 5 We assume that the fourth-order moments of the initial distribution are bounded, i.e., E[∥x 0 ∥ 4 2 +∥ȳ 0 ∥ 4 2 ] < ∞.
Our main theorem establishes the linear convergence of the Markovian process (x t , y t , ξ t ) t≥0 in W2 -distance to a unique stationary distribution: Theorem 3.3 Suppose Assumptions 1-3 hold, and step sizes α, β satisfy Assumption 4. If we start from an arbitrary initial distribution (x 0 , y 0 , ξ 0 ) ∼ µ 0 satisfying Assumption 5, then there exists a unique stationary distribution µ such that the process (x t , y t , ξ t ) ∼ µ t linearly converges in W2 -distance:
Furthermore, there exists vectors bx i , by
and variances of x ∞ and y ∞ are bounded by
A few remarks follow below. First, the theorem states that any sequence following TTSA (1) converges to some unique stationary distribution depending on problem instances and step sizes. Given the existence of the unique stationary distribution µ, henceforth, we can define random variables from the limit distribution (x ∞ , y ∞ , ξ ∞ ) ∼ µ.
Second, the limit distribution has a bias, whose dominating term grows linearly with the stepsizes. The β-wise growth in the bias of faster iterates y t is not surprising in light of known results for Markovian single-timescale SA [33,24]. More interesting is the bias of the slower iterates x t , which also grows linearly with β, even though the size of the update is only O(α) in each slow iteration. This is a unique phenomenon of two-timescale SA: the slower iterate effectively views the error from faster iterates, y t -y * (x t ), as additional "biased" noise.
Finally, the theorem shows that the limit distribution of slower iterates has an interesting property: the bias in x (slower iterates) is dominated by the faster step-size β, while its variance only scales with the slower step-size α. This is another key property of two-timescale SA that has been overlooked in prior work. In particular, we can deduce that the asymptotic MSE of slower iterates is resulted from two factors:
Focusing separately on the two iterates, we have the following more fine-grained convergence results:
Then for all t ≥ 0, we have the bounds
This corollary explicitly states how the optimization error decays from arbitrary initial points, and will be used in showing the convergence of tail-averaging next.
Using the explicit characterization of bias and variance in Theorem 3.3, we derive improved convergence rates for tail-averaging and extrapolation.
We first consider the tail-averaging variant of Polyak-Ruppert averaging [26]:
where t 0 ≳ α -1 is the length of the warm-up period. With the result from Theorem 3.2, we can analyze the MSE of tail-averaged sequence:
Theorem 3.5 Suppose Assumptions 1-5 hold and t 0 > C(αµ x ) -1 for some sufficiently large absolute constant C > 0.
Then for all t > t 0
In the above result, we omitted an additional optimization error exp(-cαµ x t 0 ) since it is dominated by other terms with t 0 ≫ 1/(αµ x ). As we can observe, O(β 2 ) is attribted to the squared-bias, and O(1/t) convergence is the variance decaying effect of tail-averaging. We also observe that the faster iterates has extra O( 1 t β 2 /α)-term. In part, this is because we measure the MSE of ỹt from y * = y * (x * ), not from y * (x t ). However, we are not fully aware whether this is an artifact of an analysis, or can be removed, and we leave the question as an open problem. Note that when β 2 ≤ α, both iterates enjoy the same O(1/t)-decaying rate of variances as if the two iterates are decoupled.
When tail-averaging can reduce the variance, extrapolation can reduce the biases of each iterate. As one example, using the fact that biases of iterates grow linearly with step sizes, we can extrapolate two sequences, (x α,β t , y α,β t ) and (a) Absolute error in the slower timescale. (b) Absolute error in faster timescale. (x 2α,2β t , y 2α,2β t
) with pairing stepsizes (α, β) and (2α, 2β). The extrapolated iterates are computed as
As a corollary of our main theorems, we have the following result characterizing the MSE of the extrapolated sequences. Extrapolation achieves reduced biases by canceling out the leading α and β terms in the asymptotic biases (6), improving the MSE bounds of both iterates from β 2 to β 4 . Corollary 3.6 Suppose Assumptions 1-5 hold and t 0 > C(αµ x ) -1 for some sufficiently large absolute constant C > 0. Then for all t > t 0 ,
Remark 1 If one uses pairing stepsizes (α, β) and (α, 2β), then only the leading β terms in the asymptotic biases (6) are cancelled.
Remark 2 It is possible to further reduce the order of bias via higher-order extrapolation using more than two sets of stepsizes as in [24,25], though it comes at the price of potentially slower convergence and higher variance due to using additional stepsizes [15,40].
We consider the TTSA iteration (1) in dimension d x = d y = 2 driven by a 10-state, irreducible, aperiodic Markov chain. We construct the transition matrix randomly and choose J 11 , J 12 , J 21 , J 22 such that Assumption 1 hold.
We tested the dependence of the bias and variance of both iterates with respect to α and β by varying each individually. After the tail-averaged iterates converged, we calculated the bias as the average distance between the averaged iterate and the true solution, and calculated the variance Tr V(•) as the average square distance from the iterate to the sample mean of the iterates. For the dependence on β, we held α constant and varied β between 0.03 and 0.07. For the dependence on α, we held β constant and varied α between 0.0001 and 0.0005.
For the slower iterate x t , Figure 1 shows that the bias scales with both β and α, while the variance is dependent mostly on α only. For the faster iterate y t , Figure 2 shows that the bias depends both on β and α, and the variance depends on β. Both results are consistent with our theory.
We also tested the effects of tail-averaging (TA) and Richardson-Romberg (RR) extrapolation with a similar setup. We fixed α = 0.0003 and let β = {0.01, 0.02, 0.04, 0.08}. In Figure 3, for each β, we plotted the absolute errors achieved by tail-averaging at stepsize β (labeled as "TA β = stepsize"), as well as the errors achieved by RR extrapolation with stepsizes β and 2β (labeled as "RR β = stepsize, 2 * stepsize"), which aims to cancel the β term in the bias. Compared to the TA iterates (solid line), the corresponding RR extrapolated iterate (the dashed line of the same color) achieved lower errors, corresponding to reduced asymptotic biases.
In addition, we examined the effectiveness of applying RR extrapolation to cancel both the α and β bias terms. Letting α = 0.0003, β = 0.02, we compared RR extrapolating on only β (using stepsizes β, α and 2β, α) with RR extrapolating on both β and α (using stepsizes β, α and 2β, 2α). In Figure 4, we see that while the former (red curves) already reduced a large amount of the bias, the latter (black curves) reduced it even further, as predicted by our theoretical results.
We outline the proofs of our main theorems. We focus on slower iterates; similar ideas apply to faster iterates.
The first step is to analyze the descent formula for each iterate separately. For the slower iterate, we have
The term T 1 would have been 0 if the noise sequence were martingale, and can be effectively handled with Markovian noises in a standard way by exploiting Assumption 2. More pressing issue is handling T 2 : with naively applying Young's inequality to bound (ii), i.e., with ⟨x t , J 12 ȳt ⟩ ≤ (c∥x t ∥ 2 + Jmax 4c ∥ȳ t ∥ 2 ), the asymptotic error easily end up being O(β 2 /α) as in [9,17], and such an approach can be improved up to at best O(β) [11].
Recent results in [28,20] directly analyzed the descent behavior of ∥T 2 ∥ op , achieving O(α) asymptotic error for the slower iterate. However, using operator norm often results in extra dependence on dimensions d x , d y , despite the smoothness condition J max = O(1) in operator norm.
Our tweak for this issue is simple: to track the convergence of cross-correlation norm, we employ the Schatten S 1 -measure for ∥Q
term is incorporated to ensure decreasing Lyapunov potential with asymmetric operators. The S 1 -norm is the best suited for exploiting the smoothness condition without incurring dimension dependence, thanks to the Holder's inequality for matrix Schattern norm:
Leveraging this property, we can construct the potential function as the sum of three terms (omitting constants):
With similar techniques for analyzing the faster iterates and cross-correlation norms, we can obtain a clean O(α) asymptotic error without additional dimension dependence. The full proof is given in Appendix B.
Once we have the MSE convergence result, extending the strategies in prior work for the single-timescale SA [10,24], we first consider two coupling sequences via sharing the common noise sequence (x
1 t , y 1 t , ξ t ) and (x 2 t , y 2 t , ξ t ). The idea is to show that the coupled sequences δ x t := x1 t -x2 t , δ y y := ȳ1 t -ȳ2 t converge linearly (Lemma B.1),
Then we can design two sequences coupled in such a way that (x 2 t , y 2 t , ξ t )
Combining the two results, the sequence (x 1 t , y 1 t , ξ t ) converges in L 2 -Wasserstein distribution to a unique stationary point. The remaining details can be found in Appendix B.2.
Turning to the stationary distributions of the iterates, we observe that x ∞ satisfies
Conditioned on the event ξ ∞+1 = ξ, we have ξ ∞ ∼ P † (•|ξ), where P † is the adjoint of the transition kernel P. Using this relation, we can construct a stationary equation for E[x ∞ |ξ ∞ = ξ], and find the explicit expression for biases by integrating the conditional expectation over a stationary distribution π, i.e.,
The variance of ȳ∞ is relatively simple to bound:
However, showing the variance upper bound O(α) can not be derived in the same fashion since the MSE bound for x∞ is O(α + β 2 ). To derive this, we also construct a stationary equation for the covariance:
and show that S 1 -norm of the above is O(α). Using the inequality Tr(A) ≤ ∥A∥ 1 completes the proof.
Lemma A.1 For any two real matrices A, B, we have
Lemma A.2 For a positive definite matrix Q ≻ 0 and any real matrix A, the following holds:
Lemma A. 3 For a positive definite matrix Q ≻ 0, and for any vectors x, y and a matrix M ,
is the condition number of Q.
Let -A be a Hurwitz matrix and Q be the solution to
Then for all ϵ ∈ 0,
, for any matrix B, we have
Lemma A. 5 For any two positive definite matrices Q 1 , Q 2 and a vector x, we have
We list some useful facts and lemmas here.
Lemma A. 6 For any t ≥ τ , for all i, j ∈ {1, 2}, we have
where v t , u t are any vectors that can be constructed at the t th iteration. Lemma A. 7 Let two intermediate variables:
Then, w x t , w y t can be rewritten as
Lemma A.8 For any t ≥ τ ≥ τ α , we have
The following corollary is convenient:
Lemma A.10 For any t ≥ τ ≥ c log( κy αµy ) with an absolute constant c > 0, we have
Similarly, we can derive the same upper bound for ∥E[w y t ȳ⊤ t ]∥ 1 .
Lemma A.11 For any t ≥ τ ≥ c log( κx αµx ) with an absolute constant c > 0, we have
Similarly, we can derive the same upper bound for ∥E[w
We recall the definition of σ x , σ y in ( 2)
Recall that we assume β/α ≫ κ y in Assumption 4, and ∥∆∥ op ≤ J max κ y .
The proof first investigates the convergence of three terms
Then, by constructing the potential function as the following:
and show that they decay in exponential rates.
We first study the descent behavior of ȳt :
We bound each term:
1. Using Lemma A.4, we have
2. Using the formula in Lemma A.7,
where we also used Lemma A.3 to have
Qy , and κ(Q y ) = O(κ y ). 3. Using Cauchy-Schwarz inequality, we have
4. We separate the cross-product term across ȳt and xt :
|.
For (i), we can derive that
For (ii), we can simply apply Cauchy-Schwarz inequality with
For (iii), we bound the term as
Combining (i)-(iii), we get
5. For the cross-product with noise, we get
6. For the last term, we simply apply Cauchy-Schwartz inequality and use inequalities used before:
Hence to summarize, with α ≪ β/κ 3 y and β ≪ 1/(J max κ 2 y ), we get
Then we can invoke Lemma A.10 with τ = O(τ α ), and noting that βJ 2 max κ y ≪ µ y to conclude that
We start with taking squared-∥ • ∥ Qx norm for the slower iterates:
Following the similar steps for the analysis of ȳt , we show the followings:
1. The main drift term satisfies
2. For the squared terms,
3. For the cross-product term,
4. For the product term with noise, we have
Writing down the intermediate result, with α ≪ 1/(J max κ x κ 3 y ), we have
Invoke Lemma A.11 with τ = O(τ α ), and we can conclude that
We start with unfolding the equation:
βJ 22 )ȳ t x⊤ t (I -α∆) -α(I -βJ 22 )ȳ t (J 12 ȳt + w x t ) ⊤ -αJ -1 22 J 21 (J 12 ȳt + ∆x t + w x t )x ⊤ t -βw y t x⊤ t + α 2 J -1 22 J 21 (J 12 ȳt + ∆x t + w x t )(∆x t + J 12 ȳt + w x t ) ⊤ + αβ • w y t (∆x t + J 12 ȳt + w x t ) ⊤ . The target norm is ∥ • ∥ 1 bound on the expectation of the cross-product term. The trick is to multiply Q 1/2 y from left on both sides, and use identity
We observe the following:
∥ op = ∥I -βJ 22 ∥ Qy ≤ 1 -µ y β, and therefore
). We omit some algebraic details, and state the desired bounds:
Applying Lemma A.11 and A.10, and using α ≪ β/κ 2 y in Assumption 4, we can conclude that
Recall the potential function V t and U t in (10):
We note that β ≪ 1/(κ 3 y κ x J max ) and β/α ≪ 1/(κ 3 y κ x ) in Assumption 4, and
Putting altogether, we have
Solving this recursion,
x , for all t, hence for all sufficiently large t ≫ α -1 , we have bounds for
Next, we consider the potential for faster iterates U t . We see that
E[∥y t ∥ 2 Qy ] ≤ U t ≤ exp(-βµ y t/2)U 0 + βκ 4 y τ α exp(-αµ x t/4)V 0 + O κ 2 y τ µ 2 y α β + τ 2 J max κ 4 y µ 2 x β ασ 2 x + O τ µ 2 y β σ 2 y , assuming βµ y ≫ αµ x . Thus for sufficiently large t ≫ α -1 log(1/β), we have E[∥ȳ t ∥ 2 Qy ] = O P (α 2 /β + αβ)σ 2 x + O P (β)σ 2 y . This concludes the final error rates as t → ∞.
Showing the distributional convergence consists of two steps. First, we setup two sequences {(x 1 t , y 1 t , ξ t )} t≥0 , {(x 2 t , y 2 t , ξ t )} t≥0 coupled with the same sequence of Markovian states {ξ t } t≥0 . We show that these two sequences will converge in the squared-L 2 expectation sense:
t , y 1 t , ξ t ) and (x 2 t , y 2 t , ξ t ), the following holds:
We first prove the above lemma and use it to conclude that the distribution of iteration variables converges in Wasserstein distance to a unique stationary distribution.
B.1 Let us define δ x t = x1 t -x2 t , δ y t = ȳ1 t -ȳ2 t . Then the stochastic recursion (1) becomes δ x t+1 = (I -α∆)δ x t -αJ 12 δ y t -αδ wx t , and for y, we have δ y t+1 = (I -βJ 22 )δ y t -αJ -1 22 J 21 (J 12 δ y t + ∆δ x t ) -(αJ -1 22 J 21 δ wx t + βδ wy t ),
where the noise differences are given by:
where we used the expression in Lemma A.7. This can be considered as the same TTSA recursion with σ x = σ y = 0. Therefore, the remaining steps are equivalent to the pilot result with σ x = σ y = 0, and it leads to
where we define
This shows that (x 1 t , ȳ1 t ), (x 2 t , ȳ2 t ) converges exponentially fast with the noise coupling, which in turn means (x 1 t , y 1 t ) and (x 2 t , y 2 t ) couples exponentially fast since
The steps here mostly follows the proof steps in [24], Appendix A.2.2. We first consider a sequence (ξ 1 t , x 1 t , y 1 t ) t≥0 that starts at (x 1 0 , y 1 0 , ξ 1 0 ) ∼ µ 0 sampled from some initial distribution µ 0 where ξ 1 0 ∼ π and (x 1 0 , y 1 0 ) are statistically independent. Then, we similarly define the initial point distribution of the second sequence (x 2 -1 , y 2 -1 ) as the same as (x 1 0 , y 1 0 ) and set (x 2 0 , y 2 0 ) be the result of one-step stochastic recursion (1), where ξ 2 -1 ∼ P † (•|ξ 1 0 ). Then we couple the Markovian states ξ 1 t = ξ 2 t for all t ≥ 0. Now that we have
1 follws a stationary distribution π, and ξ 1 t = ξ 2 t is coupled. Then by definition of Wasserstein distance (with the optimal coupling), using Lemma B.1, we get
and therefore (omitting superscript)
with α > 0. The probability space over Ξ × R dx × R dy equipped with W2 -norm is known to be a Polish space where every Cauchy sequence converges ( [44], Theorem 6.18). Furthermore, convergence in Wasserstein distance implies weak convergence ( [44], Theorem 6.9), hence weak convergence to some distribution µ ∈ P 2 (Ξ × R dx × R dy ).
We take conditional expectation on ξ ∞+1 = ξ, we have the backward conditional probability ξ ∞ ∼ P † (•|ξ ∞+1 = ξ).
Let T : Ξ × R * → Ξ × R * an unnormalized Markov operator over ξ:
Using the above notation, we can rewrite the recursion as
Let Π = 1 π, and note that Π{W ij } = 0 for all i, j ∈ {1, 2}. We eventually want to characterize zx := E[x ∞ -x * ] = π{z x } = ξ z x (ξ)dπ(ξ) and zy := π{z y }. Since πP † = π by the time-reversing property of the geometrically mixing chain, this implies ∆z x + J 12 zy + π{w x } = 0,
To further proceed, let δ x (ξ) = z x (ξ) -π{z x }, δ y (ξ) = z y (ξ) -π{z y }, and since (P † -Π){z x } = (P † -Π){δ x }, we can observe that
Then we note that
Plugging this back into (17) yields
In turn, we have
Rearranging ( 18) yields
The operation (I -P † + Π) is invertible (see Corollary B.6), and thus we can invert the operator (I -P † + Π). Putting all relationships together leads us to the recursion:
where
are independent of the choice of α, β. Next, we bound the norm of δ x , δ y , d x , d y , and thus the norm of zx , zy .
Before we proceed, we define the notion of norms that we use in the proof. For vector-valued quantities, let us define ∥v∥ L 2 (π) as
and for the Markov kernel T ,
For matrices, we use the conjugate norm-pair ∥ • ∥ 1 and ∥ • ∥ ∞ = ∥ • ∥ op . Specifically, for matrix-valued quantities, we define ∥A∥ S 1 (π) as
and
The following holder's inequality is crucial to obtain dimension-free bounds on variances: 4 For Markov kernel T and conditional matrix A(ξ), We have
where
Using the results from Markov chain literature, we have the following lemma:
Lemma B.5 (Proposition 22.3.5 in [13]) Let P be a Markov Kernel on a Boral state-space Ξ with invariant probability π. Under Assumption 2, we have
The following is the corollary:
We show that ∥δ x ∥ L2(π) = O P (α), ∥δ y ∥ L2(π) = O P (β). First, we note that
where
The last inequality is because
by Theorem 3.2, and similarly, we can also show that ∥z x ∥ L 2 (π) = O P (α + β 2 ). Furthermore,
for the same problem-dependent constants C 1 , C 2 , C 3 as defined above. This concludes that for α ≪ β ≪ 1/ max(C 1 , C 2 ), we have
Similarly, we can show that
Similarly, we have ∥ by 1 ∥ 2 = O P (1) and bx i = O P (1) for i = 1, 2. We can plug this result back to (19) to conclude the bias part of Theorem 3.3.
We note that the variance of x ∞ is measured by
where the expectation is taken over the stationary distribution (x ∞ , y ∞ , ξ ∞ ) ∼ µ, and thus we aim bound Tr(V(x ∞ )). For y, it is sufficient to bound y ∞ by O(β). To see this, note that
To proceed, we also get the expression for Σ xy :
where
and C ′ is appropriately defined. The steady-state equation is
and the system equation is
Combining these results, we can conclude that
Now plugging this back to (24), we have
Then using these results, from (23), we can derive that
Therefore, we can conclude that ∥ Σx ∥ 1 = O P (α). Since ∥ Σx ∥ 1 = Tr( Σx ) = Tr(V(x ∞ )), we obtain the last part of the theorem.
This is in fact a corollary of Lemma B.3. To see this, apply Lemma B.3 with µ 2 0 = µ, and then note that under optimal coupling between µ 0 and µ,
Similarly,
Then from Theorem 3.3, applying Tr(V(x ∞ )) = O P (α) and Tr(V(ȳ ∞ )) = O P (β), we have the lemma.
First, let us define V k , U k as:
For all k ≥ t 0 ≫ (αµ x ) -1 log(1/(αµ x )), with Theorem 3.4, we ensure that under an optimal coupling,
and similarly,
Slower Iterate: We want to analyze
where x ∞ ∼ µ(x). To show that this quantity is O(α), it suffices to bound
Qx ] under the optimal coupling. Rewriting this term,
We first note that by Theorem 3.4, under the optimal coupling between x k and x ∞ , we get
To proceed, we note that
To bound the second term, we first note that for any k ′ > 0, we use an optimal coupling between x k+k ′ |F k and x ∞ , and again apply Theorem 3.4:
Using Cauchy-Schwarz inequality, we have
where we used that V k = O P (α) for all k ≥ t 0 . Plugging this, we can conclude that
Faster Iterate: In this case, we first note that
where ỹk := 1 t-t0 t t ′ =t0 ȳt . The second term is squared-bias in order O P (β 2 ), and the third term inherits the error analysis from slower iterates. Thus, we focus on bounding the first term.
Following the same process for slower iterates, we first note that
For the first term, we invoke Corollary 3.4, under optimal coupling, we have
For the second term, with Corollary 3.4, we have
and again using Cauchy-Schwarz inequality, we can show that
On the other hand, we can apply Cauchy-Schwarz inequality in different ways:
Summarizing the results, we can conclude that
This concludes Theorem 3.5.
We note that
Note that from Theorem 3.3,
The first and second terms in ( 25) can be bounded by O P (1)/(t -t 0 ), following exactly same steps in the proof of Theorem 3.5. The result for faster iterates can also be derived similarly.
The stochastic approximation equation becomes
Using H(x * ) = 0, G(x, y * (x)) = 0, we can rewrite the recursion as (5).
We can start with a coarse bound on ∥ȳ t+1 ∥
2 Qy : ∥ȳ t+1 ∥ 2 Qy ≤ ∥(I -βJ 22 )ȳ t ∥ 2 Qy + α 2 ∥J -1 22 J 21 (J 12 ȳt + ∆x t + w x t )∥ 2 Qy + β 2 ∥w y t ∥ 2 Qy + 2α ⟨(I -βJ 22 )ȳ t , -J -1 22 J 21 (J 12 ȳt + ∆x t + w x t )⟩ Qy + 2β ⟨(I -βJ 22 )ȳ t , -w y t ⟩ Qy + 2αβ ⟨J -1 22 J 21 (J 12 ȳt + ∆x t -w x t ), w y t ⟩ Qy ≤ (1 -βµ y /2)∥ȳ t ∥ 2 Qy + O P (β 2 )∥x t ∥ 2 2 + O P (α)σ 2 x + O P (β)σ 2 y .
Thus, taking the square on both sides, we get
Given the above results, let c 2 > 0 be another sufficiently large absolute constant. Now if t ≤ c 2 log(β/α)/(αµ x ), then we set τ = c 1 τ α /4 such that t > 4τ . With TV(µ 1) ≤ exp(-αµ x t/4), we have
On the other hand, if t > c 2 log(β/α)/(αµ x ), then we take τ = t/8. In this case, we instead invoke the MSE result in Theorem 3.2, which gives
Together with TV(µ 1 τ (ξ 1 τ ), µ 2 τ (ξ 2 τ )) < ρ t/8 ≤ exp(-αµ x t/4), we get the same conclusion that
The inequality for y in ( 16) can also be similarly proven.
This result follows immediately from the fact that Tr(AB) ≤ ∥AB∥ 1 , and Hölder's inequality applied to matrix p-Schattern norm.
By definition of ∥A∥ Q , we start from
By definition for ⟨
Next, we observe that
Then since Qxx ⊤ is a rank-1 matrix, we have
Then, note that
yielding the proof.
This can be shown from the definition of Q-norm:
Let π τ = P ξt (•|F t-τ ) be a distribution over Ξ. By Assumption 2,
Furthermore, we know that Ξ W ij (ξ)dπ(ξ) = 0.
Thus,
The second inequality also follows similarly.
Note that x t = xt + x * and y t = ȳt -J -1 22 J 21 xt + J -1 22 J 21 x * . Plugging these to w x t = W 11 (ξ t )x t + W 12 (ξ t )y t + u 1 (ξ t ) and similarly to w y t yields the expressions.
By the recursion in (5),
Similarly, we have ∥ȳ t+1 ∥ 2 ≤ (1 + βJ max )∥ȳ t ∥ 2 + β(J max κ y ∥x t ∥ 2 + σ y ) + κ y α(κ y J max ∥x t ∥ 2 + J max ∥ȳ t ∥ 2 + σ x ), ∥ȳ t+1 -ȳt ∥ 2 ≤ β(J max ∥ȳ t ∥ 2 + J max κ y ∥x t ∥ 2 + σ y ) + κ y α(κ y J max ∥x t ∥ 2 + J max ∥ȳ t ∥ 2 + σ x ).
Adding two equations, κ y ∥x t+1 ∥ 2 + ∥ȳ t+1 ∥ 2 ≤ (1 + 2βJ max )(κ y ∥x t ∥ 2 + ∥ȳ t ∥ 2 ) + βσ y + ακ y σ x .
3. For the last term, similarly, ∥E[W 11 (ξ t )x t (ȳ t -ȳt-τ ) ⊤ ]∥ 1 ≤ W max E[∥x t ∥ 2 ∥ȳ t -ȳt-τ ∥ 2 ] ≤ O(1)βτ J 2 max • E[(κ y ∥x t ∥ 2 + ∥ȳ t ∥ 2 )∥x t ∥ 2 ] + O(1)τ J max E[∥x t ∥ 2 ](ακ y σ x + βσ y ) ≤ O(1)βτ J 2 max E[κ y ∥x t ∥ 2 2 + ∥ȳ t ∥ 2 2 ] + O(1)τ ((α 2 /β)κ y σ 2 x + βσ 2 y ). Combining these inequalities, and given that βτ ≪ κ -2 y /J max in Assumption 4, we can conclude that ∥E[W 11 (ξ t )x t ȳ⊤ t ]∥ 1 ≤ µ y 16κ y E[∥ȳ t ∥ 2 2 ] + O(1)βτ J 2 max κ y E[∥x t ∥ 2 2 ] + O(1)τ (α 2 /β)κ y σ 2 x + βσ 2 y . Similarly, we can also show that ∥E[W 12 (ξ t )ȳ t ȳ⊤ t ]∥ 1 ≤ µ y 16κ y E[∥ȳ t ∥ 2 2 ] + O(1)βτ J 2 max κ y E[∥x t ∥ 2 2 ] + O(1)τ ((α 2 /β)κ y σ 2 x + βσ 2 y ). For the last one, we proceed as ∥E[u 1 (ξ t )ȳ ⊤ t ]∥ 1 ≤ ∥E[u 1 (ξ t )ȳ ⊤ t-τ ]∥ 1 + ∥E[u 1 (ξ t )(ȳ t -ȳt-τ ) ⊤ ]∥ 1 ≤ ρ τ u max E[∥ȳ t-τ ∥ 2 ] + u max E[∥ȳ t -ȳt-τ ∥ 2 ] ≤ O(1)(ρ τ + βτ J max )u max E[κ y ∥x t ∥ 2 + ∥ȳ t ∥ 2 ] + τ u max (ακ y σ x + βσ y ) ≤ O(1)βτ J 2 max E[κ y ∥x t ∥ 2 2 + ∥ȳ t ∥ 2 2 ] + O(1)τ ((α 2 /β)κ y σ 2 x + βσ 2 y ),
where in the last inequality, we use u max ≤ σ y . Combining all the above inequalities yields the lemma.
To begin with, we start with unfolding the expression as
E[w x t x⊤ t ] = E[W 11 (ξ t )x t x⊤ t ] + E[W 12 (ξ t )ȳ t x⊤ t ] + E[u 1 (ξ t )x ⊤ t ]. To proceed, we first note that E[W 11 (ξ t )x t x⊤ t ] = E[W 11 (ξ t )x t-τ x⊤ t-τ ] + E[W 11 (ξ t )(x t -xt-τ )x ⊤ t-τ ] + E[W 11 (ξ t )x t (x t -xt-τ ) ⊤ ].
For each term in the above, we have the following inequalities:
1. Using the mixing-time assumption, we can show that
2. For the next term, we apply Lemma A.8 and Corollary A.9:
1)ατ J max W max • E[(κ y ∥x t ∥ 2 + ∥ȳ t ∥ 2 )∥x t-τ ∥ 2 ] + O(1)ατ W max (σ x + βJ max σ y )E[∥x t-τ ∥ 2 ] ≤ O(1) • ατ J 2 max E[κ y ∥x t ∥ 2 2 + ∥ȳ t ∥ 2 2 + τ 2 α 2 (σ 2 x + β 2 J 2 max σ 2 y ) ≤ µ x 32κ x E[∥x t ∥ 2 2 ] + O(1)ατ J 2 max E[∥ȳ t ∥ 2 2 ] + (ατ ) 3 (σ 2 x + β 2 J 2 max σ 2 y ), where we use the condition that ατ ≪ 1/(J max κ 2 x ) and W max ≤ J max . Combining these inequalities, and given that βτ ≪ κ -2 y /J max in Assumption 4, we can conclude that ∥E[W 11 (ξ t )x t x⊤ t ]∥ 1 ≤ µ x 16κ x E[∥x t ∥ 2 2 ] + ατ J 2 max E[∥ȳ t ∥ 2 2 ] + ατ σ 2 x . We also need to check the cross term: E[W 12 (ξ t )ȳ t x⊤ t ] = E[W 12 (ξ t )x t-τ ȳ⊤ t-τ ] + E[W 12 (ξ t )(ȳ t -ȳt-τ )x ⊤ t-τ ] + E[W 12 (ξ t )ȳ t (x t -xt-τ ) ⊤ ]. First term can be bounded similarly using the geometric mixing assumption. For the second term, ∥E[W 12 (ξ t )(ȳ t -ȳt-τ )x ⊤ t-τ ]∥ 1 ≤ W max E[∥(ȳ t -ȳt-τ )x ⊤ t-τ ∥ 1 ] ≤ O(1)βτ J 2 max • E[(κ y ∥x t ∥ 2 + ∥ȳ t ∥ 2 )∥x t-τ ∥ 2 ] + O(1)τ J max (ακ y σ x + βσ y )E[∥x t-τ ∥ 2 ] ≤ µ x 32κ x E[∥x t ∥ 2 2 ] + O(1) β 2 τ 2 J 4 max κ x µ x E[∥ȳ t ∥ 2 2 ] + O α 2 τ 2 J 2 max κ 2 y κ x µ x σ 2 x + κ x β 2 τ 2 J 2 max µ x σ 2 y . For the third term, similarly, we have ∥E[W 12 (ξ t )ȳ t (x t -xt-τ ) ⊤ ]∥ 1 ≤ W max E[∥ȳ t (x t -xt-τ ) ⊤ ∥ 1 ] ≤ O(1)ατ J 2 max • E[(κ y ∥x t ∥ 2 + ∥ȳ t ∥ 2 )∥ȳ t ∥ 2 ] + O(1)τ J max (ασ x + β 2 J max σ y )E[∥ȳ t ∥ 2 ] ≤ µ x 32κ x E[∥x t ∥ 2 2 ] + O(1)τ J 2 max (α + β 2 J max )E[∥ȳ t ∥ 2 2 ] + O ατ σ 2 x + β 2 τ J max σ 2 y .
For the last one, we proceed as
≤ (ρ τ + ατ J max )u max E[κ y ∥x t ∥ 2 + ∥ȳ t ∥ 2 ] + τ u max ασ x ≤ ατ J
2 max E[κ y ∥x t ∥ 2 2 + ∥ȳ t ∥ 2 2 ] + τ ασ 2 x ,
where in the last inequality, we use u max ≤ σ x . Combining all the above inequalities yields the lemma.