Introduction

The study of concentration inequalities has now become textbook material, with a variety of applications [6]. Two of the most widely used and studied inequalities for scalar-valued random variables in the independent setting are Hoeffding's inequality [28], which operates under a boundedness assumption, and Bernstein's inequality [5], which operates under an additional bounded variance assumption.

With the wide usage of these inequalities, various generalizations have been made. One line of work has sought to relax the independence assumption, deriving concentration inequalities for sums of Markov-dependent random variables. Starting with the result of Gillman [19], improvements and refinements have been made to address the general setting of time-independent functions of a nonreversible Markov chain on a general state space [15,26,33,41,42,43,54], resulting in Hoeffding and Bernstein-type concentration inequalities that generalize nicely the independent setting; these results are largely spectral in nature.

Another line of work has sought to generalize the scalar setting to concentration inequalities of sums of independent random matrices, beginning with the work of Ahlswede and Winter [2]. Several authors have followed the Ahlswede-Winter approach to develop analogs of Hoeffding's and Bernstein's inequalities for sums of random matrices; some works include [8,24,52,53,57]. However, results using this framework often have variance proxies that are suboptimal. To this end, [63,64] developed techniques to circumvent this barrier, allowing for concentration inequalities that are much tighter in many cases. These results are powerful and easy to use [65,66], with various applications [13,45,47].

We combine the above two lines of study in our work to develop concentration inequalities -Hoeffding and Bernstein-type -for the largest eigenvalue of sums of Markov dependent Hermitian random matrices, by combining the Ahlswede-Winter style argument for random matrices with the spectral techniques used to study Markov dependent scalar-valued random variables. The starting point is a deep result of Garg et al. [16], which provides a multi-matrix Golden-Thompson inequality that supports a spectral approach. Our results are broad and general: we provide inequalities for timedependent matrix-valued functions of a nonreversible Markov chain on general (continuous) state spaces. In this way, our Hoeffding-type bounds vastly generalize previous work and provide sharper constants, and our Bernstein-type bounds are new in the literature (Section 6). Along the way we prove a novel Perron-Frobenius type limit lemma for the matrix setting (Lemmas 5.4, 5.5), necessitated by our assumptions of nonversibility of the chain on continuous state spaces, that may be of independent interest.

Thus, this work can be seen as a first systematic study of concentration inequalities in the Markovdependent matrix setting. Each of our results involve two bounds: one bounding a noncommutative moment generating function of the form

and one bounding the tail probabilities

for Hermitian random matrices X 1 ,... ,X n arising as functions of a Markov chain. We expect our results on the former to be optimal, directly generalizing the scalar setting. These results are readily applicable to a wide variety of settings. For example, the best known guarantees for offline principal component analysis (PCA) are derived from a combination of the standard matrix Bernstein's inequality along with Wedin's perturbation theorem [32,38,69]. Now, if the sample matrices are not independent, but instead arise from a Markov chain, we again have a bound for leading eigenvector estimation using our bounds. Our results are also directly applicable to sums of matrix-valued random variables sampled from a random walk on an expander; Theorem 2.5 improves on the best known results in this setting [16,70]. There is also a line of work in machine learning that realizes stochastic gradient descent (SGD) with constant step size as a Markov chain with a stationary distribution [4,11,48]. Additionally, Hessian information has been used widely to give improved algorithms and better understand model performance [18,44]. Our results are able to give precise bounds on the convergence of the largest eigenvalue of the empirical mean of the Hessian matrices along the path of SGD to the largest eigenvalue of the expected Hessian with respect to the stationary measure. This has many interesting and relevant applications in statistical learning and generalization. More broadly, our bounds are useful in the study of online algorithms, where matrix-valued data is received and processed in a streaming fashion from an underlying Markov chain. In Section 7 we more closely focus on an application of our main results to offline PCA for samples arising as a function of an underlying Markov chain; we present a bound that is a direct generalization of the best known bound for offline PCA in the i.i.d. setting. This result was first stated in [38], where to our knowledge it was the first in the literature to give this sort of bound for offline Markov PCA ; it directly uses our Theorem 2.6.

Related work

In terms of results, the literature on concentration inequalities for sums of Markov-dependent random matrices is fairly small, in contrast to the scalar setting. Recent progress in this area builds on the work of [16,60], which develop a powerful Golden-Thompson inequality [21,62] that allows for a spectral analysis of a useful moment generating function; the latter provides a corresponding Hoeffding's inequality for when the stationary distribution is uniform and the chain is stationary. The work of [56] extends this to arbitrary stationary distributions and initial distributions; these works are valid only for finite state spaces and reversible chains. Our work generalizes these results to general state spaces and nonreversible chains using significantly different techniques, which allows us to also improve the constants in the bounds. Furthermore, our Hoeffding-type result is in terms of both the largest and smallest eigenvalues, matching the type given in the scalar setting; and we give a Bernstein's inequality, which is new in the literature.

In terms of techniques, most related to ours are the works of [15,26,33,41,42,49,54], which give Hoeffding's and Bernstein's inequalities for sums of Markov-dependent scalar random variables. Broadly speaking, the idea is to bound the largest norm of a perturbed operator. We prove the corresponding spectral properties of a Markov operator that is the tensored form of the operator appearing in these works in Section 5; in particular, we generalize to continuous state spaces using related functional analytic ideas. However, previous techniques are not sufficient for us to prove our main results; our proof of the result on the limit of the moment generating function (Lemma 5.4) takes a detour through dynamical systems theory and Banach space theory.

In contrast to the Markov-dependent setting, the sums of independent random matrices have been studied in depth, and is the model for which we give our bounds. In particular, our results mirror the type given in [2,8,24] in terms of variance proxy, and the analysis proceeds from the same starting point: bounding a matrix moment generating function; in contrast, later works of [52,63] sharpen the variance proxy. Our analysis of the moment generating function arising from the matrix-valued Golden-Thompson inequality result in bounds resembling more of the former; we believe our analysis is likely sharp given the form of the moment generating function, so an improvement in terms of the variance proxy would likely require a different approach and is an interesting open problem.

Besides the spectral approach, there have been a myriad of ways to derive concentration inequalities relaxing the independence assumption, both in the scalar and matrix settings. The papers of [1,20] exploit regeneration-type minorization conditions to derive exponential concentration for ergodic sums; though these works do not assume a nonzero spectral gap they often have less explicit constants. The work of [54] uses Marton couplings to obtain concentration inequalities for a scalar-valued function of a single random sample. The papers of [46,55] leverage Efron-Stein inequalities and exchangeability to develop concentration inequalities for random matrices, also relaxing the independence assumption. Similarly, recent works of [3,30,36] demonstrate that matrix functional inequalities, i.e., Poincaré, directly translate to matrix concentration inequalities; [36] develops a Bernstein's inequality for strongly Rayleigh distributions; [39] gives a Chernoff bound for a similar setting. These largely address single sample concentration of (Lipschitz) functionals; however, some methods can be generalized to product measures [31], recovering some of the Efron-Stein results. And the work of [16] sketches a method to reduce studying concentration of random variables sampled from a Markov chain to concentration of sums of martingale random variables [9,63]; though the resulting bounds are suboptimal, they show that qualitatively concentration for Markov chains is more generic than just the scalar or matrix settings (see [40]). These techniques may be very useful in improving the variance proxy, better matching the independent setting, where a direct spectral analysis may prove insufficient; in particular, functional inequalities have emerged as a powerful general tool to derive concentration. It may be interesting to see if this can be demonstrated in the scalar setting as well.

Main Results

In this section we give our main results. Each will hold under a combination of the following assumptions.

Assumption 2.1 P is a discrete-time Markov chain on a continuous state space X with stationary distribution µ and absolute spectral gap λ (Definition 3.1). s 1 ,... ,s n is a sequence of states driven by P with initial distribution µ. Assumption 2.2 F j : X → R d×d , j = 1,... ,n, is a sequence of functions each mapping to real symmetric d × d matrices. Each F j is ℓ 2 (µ ⊗ 1) measurable (see Section 3). Assumption 2.3 For all j, E µ [F j (x)] = 0 and a j I F j (x) b j I for all x ∈ X . Assumption 2.4 For all j, E µ [F j (x)] = 0 and there exists a constant V j such that E µ F j (x) 2 ≃ V j for all x ∈ X . Moreover, there exists an absolute constant M such that F j (x) ≃ M for all j, x ∈ X .

The following two theorems each consist of two statements. The first bounds a certain type of moment generating function arising from the expected Frobenius norm of a product of matrix exponentials. The second statement turns this into a tail bound. We only give upper tail bounds; the corresponding statement for lower tail bounds are clear and the proof proceeds in almost the exact same way.

Theorem 2.5 (Markov Matrix Hoeffding, Real) Instate Assumptions 2.1, 2.2, 2.3. Then for any θ > 0,

and for any t > 0,

Here α(λ ) = (1 + λ )/(1 ⇐ λ ), where λ is the absolute spectral gap (Definition 3.1).

The first statement in the above result reveals the sub-Gaussian nature of the moment generating function; in the scalar case, it exactly reduces to the standard moment generating function. Using our methods, this bound matches the scalar case given in [15]. The second statement provides a large deviation tail bound; the bound on the lower eigenvalue follows analogously. Our results offer four main improvements to that of [56]: first, we give strictly better constants, both in terms of d and the absolute constants in the exponent (the α(λ ) differs slightly from their equivalent term for the mixing of the Markov chain, but is strictly better and is more classical); second, we allow for time-dependent functions F j , which is a strict generalization; third, our inequality is in terms of both a lower and upper eigenvalue bound (the a j , b j parameters), giving a result that is directly comparable to the classic Hoeffding's lemma; and fourth, our bounds apply to general state spaces.

Furthermore, if θ < (1 ⇐ λ )/(8M /π), then for any t > 0,

Here

where λ is the absolute spectral gap (Definition 3.1).

This is to our knowledge the first result giving a Bernstein-type inequality for sums of random matrices arising from markov chains. The first statement in the above result is made up of two terms; the (e M θ ⇐ M θ ⇐ 1)σ 2 /M 2 term coincides with the convex function that makes up the classic Bernstein's and Bennett's inequality. The second term in the exponent reflects the influence of the Markov chain; compared to previous scalar results [33,54], our techniques are able to recover a slightly improved version of this second term in terms of absolute constants. A linear algebraic perspective of the first statement above is that the first term bounds how much the product of matrix exponentials increases the magnitude of vectors in the direction of 1, which is the leading eigenvector of the operator P. The second term is a bound on how much the product of matrix exponentials increases the magnitude of orthogonal directions, hence involving λ .

We now give two corollaries of the above results, generalizing the assumptions above; the proofs are in the appendix and are classical. Our first corollary generalizes to complex Hermitian matrices via a complexification technique [12]: Corollary 2.7 (Extension to Complex Matrices) Instate Assumption 2.1. Now assume that F j : X → C d×d , j = 1,... ,n, is a sequence of functions each mapping to complex Hermitian d × d matrices.

Under Assumption 2.3 (resp. Assumption 2.4), the conclusion of Theorem 2.5 (resp. Theorem 2.6) holds with an extra multiplicative factor of 2 on the right-hand side.

Our second corollary generalizes to the case that the Markov chain starts at a distribution ν: Corollary 2.8 (Extension to Nonstationary Chains) Let P be a Markov chain on general state space X with stationary distribution µ and absolute spectral gap λ . Let s 1 ,... ,s n be a sequence of states driven by P with initial distribution ν, where ν ⇒ µ. Instate Assumption 2.2. Let ess sup dν dµ be the essential supremum of the Radon-Nikodym derivative dν dµ . Under Assumption 2.3 (resp. Assumption 2.4), the conclusion of Theorem 2.5 (resp. Theorem 2.6) holds with an extra multiplicative factor of ess sup dν dµ on the right-hand side.

Remark 2.9 In many applications these bounds are used to determine how many samples are needed (i.e., how long the chain has to run) in order for the tail probability to be less than some fixed value; in these cases the number of samples will depend only logarithmically on ess sup dν dµ . Alternatively, the probability of an event under the Markov chain driven by P initialized at ν can be bounded by the probability of an event under the Markov chain driven by P initialized at µ up to an extra multiplicative factor of 1 + χ 2 (ν µ); see Lemma 29 of [37].

We now give some remarks on the above theorems. First, for the first statements regarding the moment generating function in both Theorems 2.5 and 2.6, the variance proxies are asymptotically optimal for a class of Markov chains. Let µ be a distribution on X and let P be a transition kernel such that P(x, y) = λ x=y + (1 ⇐ λ )µ(y), so that µ is the stationary distribution. Let f : X → R be any scalar-valued function, and let F : X → R d×d be defined as F(x) = f (x)I, so that these are all real symmetric. It is straightforward that

since all F(x) commute. Under this setup, start with our Bernstein assumptions. Let [17] states that for a two state Markov chain driven by a P of our above form that

which is called the asymptotic variance. Now suppose that V is such that E µ [ f 2 ] ≃ V ; it is classical that any such variance proxy V that satisfies g(θ ) = nV θ 2 /2 + o(θ 2 ) for a function g that upper bounds the scalar moment generating function via

upper bounds Var 1 √ n ∑ n j=1 f (s j ) [33]. This lower bound is indeed obtained by our Markov kernel P and choice of function f , and therefore the same holds for our matrix-valued function F. Now assume f are Rademacher, i.e., the probability under µ that f (x) = 1 is 1/2 and f (x) = ⇐1 is 1/2. Define F the same way as above, and so we see that E µ [F] = 0 and ⇐I F(x) I for all x ∈ X . Via Equation 2.1 we again reduce to the scalar case since F commute, and through 2.2 we have that V asy = (1 + λ )/(1 ⇐ λ ), since the variance of f is 1. Now in this case we have

as can be seen from our results with d = 1. We see that nα(λ ) is the variance proxy for a sub-Gaussian random variable and so naturally upper bounds the asymptotic variance. This lower bound is again attained for P and f , and therefore the same holds for the matrix-valued function F.

Roadmap

We begin with the quantity

for φ ∈ [⇐π/2, π/2] and θ > 0, developed first in [16,60]. This will play the role of our moment generating function. Using properties of the Kronecker product, in Section 4 we show that this moment generating function is equal to

where P is the lifted Markov operator and E θ /2 j is a certain multiplication operator (Definition 3.4). This quantity is a matrix version of a classical quantity that appears in the analysis of sums of Markovdependent scalar random variables -indeed, when d = 1 these are equivalent. Thus, we can apply a spectral analysis to bound this quantity.

In Section 5 we proceed with this study; this section addresses the main challenges of extending to continuous state spaces. Using Cauchy-Schwartz we bound the above by

is a sort of limit for a corresponding moment generating function. More precisely, we show that

where s 1 ,... ,s n is driven by the Markov chain represented by the Leon-Perron operator P; this result appears in Lemma 5.4 and 5.5. The realness assumption on the F j is important for the proof of this result -we invoke an interesting generalization of the Perron-Frobenius theorem that applies to linear transformations leaving invariant a cone. This result is our main technical innovation and is crucial in extending existing results to the continuous state space setting; conveniently it also allows us to sharpen known bounds in discrete state spaces. Section 6 gives our final bounds. We first show how to transfer bounds on our moment generating function to tail bounds; the technique is straightforward and uses the multi-matrix Golden-Thompson inequality from [16]. We then give these tight bounds on the moment generating function under both Hoeffding and Bernstein-type assumptions. For Hoeffding-type assumptions we use a coupling technique to exhibit a two-state chain that acts as the "limit" of our chain. The leading eigenvalue of the corresponding operator can be solved exactly, giving optimal bounds. For Bernstein-type assumptions, we use a robust linear algebraic approach to bound the operator norm of a related matrix directly; this approach will also give optimal results.

Preliminaries

Notation

Lower case, unbolded a denotes scalars or scalar-valued functions, bold a denotes vectors or vectorvalued functions, and upper case A denotes matrices or matrix-valued functions.

The operator ⊗ will denote the Kronecker product, where for matrices A, B of size a × b, c × d, respectively,

Markov chains

Let X be a general state space with σ -algebra B := B(X ). Let P be a Markov kernel on X defined, for a sequence of random variables X 1 ,... ,X n , in the standard way as

with stationary measure µ so that

for h : X → C measurable. This is a a Hilbert space equipped with the following inner product:

The norm of a linear operator T is defined in the usual way:

A transition kernel P acts as an operator on ℓ 2 (µ):

The projection operator Π corresponding to the distribution µ is defined as

(1 is the function such that 1(x) = 1 for all x); this is a rank-1 operator by definition, and is indeed a projection onto

An important subset of elements of ℓ 2 (µ) will be the class of "mean-zero functions," denoted

Then we have the following definition:

Definition 3.1 (Absolute spectral gap) A Markov kernel P with stationary measure µ admits an absolute spectral gap 1 ⇐ λ (P) if

Note that λ (P) ≃ 1 always, as Ph µ ≃ h µ by Jensen's inequality, with equality for h = 1. When it is clear from context we will just use λ = λ (P). The absolute spectral gap characterizes the convergence of a Markov chain to its invariant measure [59]. For reversible finite-state chains, the value λ (P) corresponds to the second largest eigenvalue and the existence of the gap corresponds to ergodicity. A Markov chain driven by P is reversible if and only if P is a self-adjoint operator on ℓ 2 (µ).

We define what is known as a Leon-Perron operator [15,33,41] that will be important in the sequel: We can interpret Pc as a transition kernel such that at a state, it stays at that state with probability c or samples a new state independently from µ with probability 1 ⇐ c. Note that if P admits an absolute spectral gap then P does so as well with the same absolute spectral gap.

Behavior in tensored space

As we work with matrices, the above will have to be lifted to a tensored space as necessary. We first consider the product space

with inner product induced from the standard Euclidean product, i.e., a ⊗ c, b ⊗ d = a, b c, d , and extend by linearity. Oftentimes we will simply identify an element as y ∈ C d 2 ; it is straightforward that the standard Euclidean inner product on this larger space is equivalent to the tensored inner product. So it is no loss to move between one representation to the other. Now define ℓ 2 (µ ⊗ 1) as the following "lift" of ℓ 2 (µ): formally, define

measurable as a vector-valued function. This is equipped with the inner product

and corresponding norm f µ =  f, f µ . This gives the understanding of ℓ 2 (µ ⊗ 1) as a direct integral of Hilbert spaces indexed by the points of X ; decomposable elements in this space are understood to be of the form

Norms of linear operators are thus defined with respect to this norm in the usual way.

For a Markov operator P on ℓ 2 (µ) we can define the operator P := P ⊗ I d 2 as an operator on ℓ 2 (µ ⊗ 1). This operator acts as

note that ( Ph)(x) is a vector. Whereas for the operator P the leading eigenfunction was 1 and spans the entire eigenspace for the leading eigenvalue of 1 (assuming nonzero spectral gap), for P the leading eigenspace is

For the projection operator Π we can define the lifted version as Π := Π ⊗ I d 2 , again satisfying Π P = P Π = Π, that performs exactly this operation, i.e., Π(

. This can all be extended by linearity and convergence in ℓ 2 . Then there is a corresponding lift of "mean-zero functions" as

and can also be identified as

Similarly, define

The following lemma, Lemma 3 from [56], will imply that if P admits an absolute spectral gap, then so does P with the same absolute spectral gap. Thus, many of the essential spectral properties of P and P are the same. The proof is deferred to the appendix.

Lastly, we can define Leon-Perron operators of P as the corresponding lifts of Leon-Perron operators of P, which is tensorization by I d 2 . With some abuse of notation, we will often interchange the notation P to refer to both Leon-Perron operators of P and Leon-Perron operators of P, but the usage will be clear from context.

Recall for a scalar-valued function g, the multiplication operator M g is defined as (M g h)(x) = g(x)h(x). For a matrix-valued function G we have the following definition: Definition 3.4 (Multiplication operator) Let G be a matrix-valued function mapping X to C d 2 ×d 2 matrices. Define M G be the operator on ℓ 2 (µ ⊗ 1) such that for any h ∈ ℓ 2 (µ ⊗ 1), (M G h)(x) = G(x)h(x). We also call these block diagonal operators.

In other words, M G is multiplication by a function G in the direct integral of Hilbert spaces [50]. We will often use the multiplication operator E θ H defined as

, sometimes dropping the subscript H when it is clear from context.

Starting point

Our starting point is the following matrix moment generating function:

where φ ∈ [⇐π/2, π/2], θ > 0, and s 1 ,... ,s n is a sequence of states driven by the Markov chain with transition matrix P.

To simplify this moment generating function, we first use the property that tr[AB ] = vec(I d ) (A ⊗ B) vec(I d ), where ⊗ is the Kronecker product and I d is the d × d identity matrix to write

and then use successive applications of the property (AC)

Lastly, we use the fact that exp(A)

where H 1 ,... ,H n is a sequence of matrix-valued function on X , implicitly with respect to φ , defined as

for x ∈ X . Note that unfortunately H j (x) is not in general Hermitian for any x, though its real and imaginary parts are both symmetric. We emphasize the identity

as we will switch back and forth from these two forms as needed.

The above allows us to focus on the matrices H j , which are measurable and satisfy many of the same probabilistic properties as F j -a more precise statement is given in the next section. We then write

Recall that E θ j the multiplication operator defined as (E θ j h)(x) = exp(θ H j (x))h(x) for any vectorvalued function h -see Definition 3.4, and recall that P = P ⊗ I d 2 . E θ j and P act via

and

as operators. Then simplifying Equation 4.5, and using the definition of the inner product on ℓ 2 (µ ⊗ 1), we have

Bounding the operator norm

Having derived the above, in this section we give a series of bounds that shift our focus to bounding the norm of a perturbed operator; we address the main technical challenges of studying general continuous state spaces.

From Equation 4.6, we apply Cauchy-Schwartz and submultiplicativity of the spectral norm to obtain

(5.1)

Now P possesses many of the same spectral properties as P, most notably the result from Lemma 3.3. Recall that P is the "lifted" Leon-Perron operator for P, defined as P = (λ

The following lemma allows us to replace P with its Leon-Perron version P.

Lemma 5.1 The operator P and its Leon-Perron version P satisfy the following:

(1) For any g, h

(3) For any multiplication operator M G with respect to a measurable matrix valued function G, and for

Proof For (1), we have

where the first inequality follows because Π P = P Π = Π and because Π is a projection and is thus self-adjoint on ℓ 2 (µ ⊗ 1). Because of this we also have

For (2), we have

The result follows by noting that

where the second to last line follows because I ⇐ Π is a projection. We can express these terms more explicitly. First, for the denominator, we have

Now for any Hermitian matrix A and any vector

Rearranging finishes the proof.

Going back to Equation 5.1, applying Lemma 5.1 parts ( 2) and ( 3), and noticing that vec

Thus, we have transferred the study of the problem from P, which represents a general, nonreversible Markov chain, to P, which represents a reversible chain and thus is much simpler to analyze. Therefore, the operators

are self-adjoint on ℓ 2 (µ ⊗ 1). Now we focus on bounding the leading eigenvalues of these operators; we make one more simplification:

. Let E θ T j be the operator defined as

Then the norms of E θ /2 j PE θ /2 * j and E θ /2 T j PE θ /2 T j are the same. Proof The operator E θ /2 j PE θ /2 * j is similar to the operator PE θ /2 * j E θ /2 j . The operator E θ /2 * j E θ /2 j acts via

where the first equality holds because

T j have the same spectrum; as they are both self-adjoint, the spectral radius equals the operator norm.

At this point, we have reduced our problem to a considerably simpler form. We can now show that the operators T j satisfy many of the same probabilistic properties as F j ; the proof is in the appendix.

where all expectations are taken with respect to a measure on X , and F and T are measurable.

The leading eigenvalue as the limit

In this subsection we will prove a limit lemma, justifying the use of the leading eigenvalue as a characterization of the moment generating function, mirroring the scalar case as laid out in [10] essentially, what we will prove is that the log limit of the moment generating function is the logarithm of the leading eigenvalue. This will allow us to shift between time-dependent and time-independent functions as needed. The main difficulty in extending this theory from the scalar case is proving the limit statement in part (3) of Lemma 5.5, and in particular, a specific lower bound in the proof. This lower bound is not difficult to prove in the scalar setting; however, the matrix setting is more subtle and requires some more care.

Our main lemma is the following:

Lemma 5.4 Let F be a map from X to real symmetric d × d matrices with aI F(x) bI for all x ∈ X . Let T (x) = 1 2 (F(x) ⊗ I d + I d ⊗ F(x)), x ∈ X , and let E θ /2 T be the operator defined as

Let µ be a distribution, let s 1 ,... ,s n be a sequence of states driven by the stationary Markov chain with transition matrix λ I + (1 ⇐ λ )1µ , and let P be this transition matrix tensored with I d 2 . Then

.

In other words, we take s 1 ,... ,s n to be driven by a Leon-Perron operator, where at a given state the chain stays at that state with probability λ and samples a state from µ with probability 1 ⇐ λ . We first address the case of a simple matrix-valued function F.

Assume that for all B j , aI B j bI, and let these be tight. Let T be the operator defined as

, with exponential multiplication operator E θ T . Define the matrix-valued function

Then the following hold:

(1) Let r * be such that F (r * ) has an eigenvalue of 1, and let v be an eigenvector corresponding to the eigenvalue 1 of F (r * ). Then r * is an eigenvalue of E

with corresponding eigenfunction h such that

Conversely, every eigenvalue of E

The third term is equal to

Adding the first two and subtracting the third gives zero, and shows that r * is indeed an eigenvalue with eigenfunction h.

T , then solving for h in the equation

Multiplying by exp θ 2 T (x)  on both sides and taking expectations, we see that

. For (2), the function F (r) is continuous in r when r is large enough (so that the inverses exist); indeed, since F is simple the expectation decomposes as a finite sum, so F is continuous if all of the summands are continuous. In particular, since the eigenvalues of all B j are bounded above by b, when r > λ e θ b F (r) is continuous in r. Now when r → ∞ the entries of F (r) approach zero, implying that all eigenvalues of F (r) also approach zero as the eigenvalues are continuous in the entries. On the other hand, F (r) decomposes as a finite sum. Consider the summand corresponding to some B j such that b is its largest eigenvalue. As r → λ e θ b from above it is clear that the largest eigenvalue of this summand goes to ∞; then it is clear that the largest eigenvalue of this finite sum must also go to ∞; implying that at least one eigenvalue of F (r) goes to ∞. Thus, as F (r) is entrywise continuous in r, and as the eigenvalues are continuous in the entries, for some r * ∈ (λ e θ b , ∞) F (r * ) has an eigenvalue of 1. By part (1). this implies that r * is an eigenvalue of

T . Now let ρ be the largest of all such eigenvalues, and so we have shown that ρ > λ e θ b .

For (3), let σ (E

is self-adjoint on ℓ 2 (µ ⊗ 1), the spectrum is real and decomposes as the disjoint union of the essential spectrum

T ), the latter consisting of all eigenvalues with finite multiplicity. Define the essential spectral radius as the largest element of the essential spectrum and the spectral radius as the largest element of the spectrum, so the latter is always at least the former.

and thus is a compact perturbation of the operator λ E θ T , as (1 ⇐ λ )E θ /2 ΠE θ /2 is finite rank (rank at most d 2 ). Weyl's perturbation theorem [67] ensures that

is a block diagonal operator (Definition 3.4), its spectrum is the union of the spectra of λ exp(θ T (x)), x ∈ X [50], and so the essential spectral radius is bounded above by λ e θ b . Therefore, by part (2). the largest eigenvalue ρ is equal to the spectral radius of

To prove the second statement, from Lemma 5.1 we have

This implies lim sup

Let h be a unit norm eigenfunction corresponding to ρ, and let h be the projection of

by definition of h being an eigenfunction and it being a projection of E θ /2

T (1 ⊗ vec(I d )). Next, we have

is positive semidefinite on ℓ 2 (µ ⊗ 1) -this follows from E θ /2 T being self-adjoint and P being positive semidefinite on ℓ 2 (µ ⊗ 1). Then

Putting Equations 5.3 and 5.5 together proves the statement. We now show that

To do this we make a connection to dynamical systems and the theory of positive operators in Banach spaces, and which uses our realness assumption on F. For a cone K in a Banach space, an operator A is positive with respect to K if AK ⊆ K. In our case, let the Banach space under consideration be the space of d × d symmetric matrices equipped with the Frobenius norm, or equivalently, the space of length d 2 vectors whose reshaping to a d × d matrix is symmetric equipped with the Euclidean norm; we will show that it is no loss to restrict to this space. Let S d + be the cone of real positive semidefinite d × d matrices; it is a closed, convex cone. It is also reproducing in the space of real symmetric d × d matrices, in that every symmetric matrix X can be written as X = X 1 ⇐ X 2 for X 1 , X 2 positive semidefinite [27]. We can also identify S d + as the space of length d 2 vectors whose reshaping to a d × d matrix is positive semidefinite; we now make this identification. Now let K ⊂ ℓ 2 (µ ⊗ 1) be the infinite direct product of S d + indexed by elements of X ; in other words, for h ∈ ℓ 2 (µ ⊗ 1), h ∈ K if and only if h(x) ∈ S d + for all x ∈ X . This cone is is closed, convex, and reproducing in the space of all h ∈ ℓ 2 (µ ⊗ 1) such that h(x) can be reshaped to a symmetric d × d matrix for all x ∈ X .

Next, we show that we can indeed restrict ourselves to the Banach space of symmetric matrices equipped with the Frobenius norm, in the sense that any eigenfunction h of E is positive with respect to K, and the spectral radius ρ is strictly larger than the essential spectral radius, which is at most λ e θ b . Then there is a nonzero eigenfunction h in K ⊆ ℓ 2 (µ ⊗ 1).

This allows us to assert that in the derivation of Equation 5.5 we could have chosen h ∈ K. With this choice of h it is straightforward in showing

 .

Each of these matrices exp θ 4 F(x)  ĥ(x) exp θ 4 F(x)  are positive semidefinite; moreover, as h is nonzero in the ℓ 2 space ĥ(x) is nonzero on a set of positive measure (i.e., ĥ(x) cannot be nonzero only on a set of measure zero). Therefore, since exp θ 4 F(x)

 is full rank, by Sylvester's inertia theorem [61] we have that exp θ 4 F(x)  ĥ(x) exp θ 4 F(x)

 has at least one strictly positive eigenvalue for all x in a set of positive measure, meaning that the trace is strictly positive. Then

Using this lemma, we are now able to give the proof of Lemma 5.4.

Proof of Lemma 5. 4 We first reduce to the simple case. The range of F is compact, as F can be seen as a continuous map of the compact set O(d) × [a, b]. Then for ε > 0 let N ε be an ε-net of the range of F with respect to the metric induced by the operator norm. Define F ε as

.

In fact, we can give a more refined proof of part (3). of Lemma 5.5 to show that lim n→∞ a ε,n uniformly converges in ε. Let ρ ε be the leading eigenvalue of E θ /2

T ε . We first have by Lemma 5.1 that

Next, we showed for any unit norm eigenfunction h ε corresponding to the leading eigenvalue ρ ε that

In the proof of part (3). of Lemma 5.5 we showed that h ε can in fact be chosen such that for all x ∈ X , h ε (x) can be reshaped to a positive semidefinite d × d matrix. With such structure we lower bounded

We can actually give a more quantitative lower bound of this quantity that is independent of ε; in addition, we can give an upper bound ρ ε also independent of ε.

First note that under an ε-net of the range of F, since aI F(x) bI, that aI F ε (x) bI also. Then the upper bound is simple; we have

For the lower bound we first write

 .

Since we have established exp

 is positive semidefinite for all x ∈ X , the trace is equal to the nuclear norm, which upper bounds the Frobenius norm. Therefore,

Now since F ε takes only finite number of values, we have that h ε can also only take a finite number of values (see part (1) of Lemma 5.5). Letting A j = F ⇐1 (B j ) and denoting h j = h ε (x) for x ∈ A j , we have

µ(A j ) h j .

Now as h ε µ = 1 we must have h j ≃ 1, and so then

→ e θ a , and we thus have a ε,n ⇐ log ρ ε → (θ (a ⇐ b))/n. Putting this all together, we have shown

which implies uniform convergence, since d, θ , a, and b are all constants independent of ε. We can then swap limits: lim

ε→0 lim n→∞ a ε,n = lim n→∞ lim ε→0 a ε,n .

But as E θ T ε converges to E Then we have the chain of equalities

Final bounds

In this section we prove our main theorems regarding the moment generating function and give concentration inequalities. We first show how to go from bounds on the moment generating function of Equation 4.1 to tail bounds.

From the MGF to tail bounds

The reason we have studied the moment generating function of Equation 4.1 is due to the following result:

Theorem 6.1 (Multi-matrix Golden-Thompson inequality, [16]) Let H 1 ,... ,H k ∈ C d×d be Hermitian matrices. Then

where ξ is some probability distribution on ⇐ π 2 , π 2 .

From this, we can obtain tail bounds. Recall that F j : X → R d×d are a sequence of measurable functions from X to d × d real symmetric matrices, and let E[F j ] = 0. First,

We would like to apply Theorem 6.1; we follow the standard outline given in [16]. An immediate application of Jensen's inequality on the right-hand side of Theorem 6.1 furnishes an upper bound, and then taking exponents, we have Combining the above two lines, adding in θ , and taking expectation, we have

Now let M F (θ ) be such that

that is independent of φ , i.e., an upper bound on the moment generating function in Equation 4.1 valid for any φ ∈ [⇐π/2, π/2]. Then

where the last line holds because ξ is a probability distribution and M F does not depend on φ . We now give our final bounds on the moment generating function

under both Hoeffding and Bernstein-type assumptions. In Sections 4 and 5 we have reduced the problem to estimating the leading eigenvalues of the operators

for j = 1,... ,n and P the lifted Leon-Perron version of P. Lemma 5.4 has shown that each eigenvalue is a limit a corresponding moment generating function. Thus, the general idea will be to bound

for each j, and where s 1 ,... ,s n is a sequence of states driven by the Markov chain corresponding to P. In these bounds we will proceed by way of defining a function M F to bound these operators.

Hoeffding

In this section, fix a matrix-valued function F, and let E µ [F(x)] = 0 and aI F(x) bI for all x ∈ X . For a scalar convex functions Ψ, we have the following proposition: Proposition 6.2 Let H be a Hermitian d × d matrix such that aI H bI. Then for any scalar convex function Ψ,

Here Ψ acts spectrally; if H = UΛU * , then Ψ(H) = UΨ(Λ)U * , where Ψ acts on the diagonal matrix Λ diagonally, i.e., Ψ(Λ) is the diagonal matrix such that Ψ(Λ) ii = Ψ(Λ ii ).

Proof We simply need to show that for any unit vector u that

It suffices to show this for eigenvectors v of H. If λ is the corresponding eigenvalue, then

Since Ψ is a scalar convex function and a ≃ λ ≃ b,

For general u, let u = ∑ d k=1 c k v k be the representation of u in the basis of eigenvectors of H. Then

If in the above H is a random matrix such that E[H] = 0, then taking expectations we have This allows us to make the following argument about a convex ordering of distributions, with proof in the appendix: Lemma 6.3 Let F be a function from X to Hermitian d × d matrices, and assume E µ [F] = 0 and aI F(x) bI for all x ∈ X . Let s 1 ,... ,s n be driven by the Leon-Perron operator P. Let Q, µ, G be as above. Let y 1 ,... ,y n be driven by Q. Then for any scalar nonnegative convex function Ψ,

Remark 6.4 The condition on Ψ above does not seem to be very strict. Indeed, it applies to many natural matrix-valued functions, such as even matrix powers and the matrix exponential.

Applying this to the scalar convex function x ∀ → exp θ cos(φ )

2

x yields

Now let M θ ∈ R 2d 2 ×2d 2 be the operator on the two-state state space {0, 1} such that

Then Lemma 5.4 indicates that

where η is the largest eigenvalue of

Then by properties of the Kronecker product, the eigenvalues of the matrix on the left-hand side are equal to the eigenvalues of K θ , a 2 × 2 matrix, the largest of which can be solved for explicitly. Putting this all together,

With the above we can establish the following result, with proof in the appendix.

Lemma 6.5 Let K θ be the 2 × 2 matrix defined in Equation 6.8, and let η be its largest eigenvalue. Then

where α :

Putting our work above together, we can now prove the two statements of Theorem 2.5.

Proof of Theorem 2.5 The discussion in Lemma 5.1 (in particular Equation 5.2) and Proposition 5.2 reveal

where recall

T j is the multiplication operator with respect to the matrix-valued function T j and P = (λ I +(1 ⇐ λ )Π) ⊗ I d 2 . Then Lemma 5.4 and Lemma 6.3 allow us to bound each of these norms with the norms of corresponding matrices arising from a two state chain, each of which can be bounded using Lemma 6.5. In other words, for each j we have, verifying that the operators T j satisfy the assumptions of the lemmas using Proposition 5.3,

Then

and setting φ = 0 we obtain the first statement of the theorem.

To turn this into a tail bound, since cos(φ ) 2 ≃ 1, we can let M F (θ ), as described in Equation 6.3, be

Using Equation 6.4, we then have

Optimizing this quadratic over θ we have

as desired.

Bernstein

In this section, fix a matrix-valued function F, and let E µ [F(x)] = 0, E µ [F(x) 2 ] ≃ V , and F(x) ≃ M for all x ∈ X . We take a slightly different approach to bounding the norms of the operators E θ /2

T j . In particular, consider again Equation 4.1 for a time-independent function T and the Leon-Perron operator P. Then

with the last line following from Π P = Π so we can introduce another P. Now from Lemma 5.4 and the above we have

Although the above may seem a bit contrived, it allows for a general recursive method used in [16,26,56,68]; more importantly, we can extend the final result to time-dependent functions via our limit lemma. In particular, for an initial function z 0 , we can track the updates z k := ( PE θ T ) k z 0 recursively, and then provide a bound. This method is strong enough to give the sorts of inequalities we are after. Now we introduce some notation. For an element z ∈ ℓ 2 (µ ⊗ 1) let z = Πz and z ⊥ = (I ⇐ Π)z. In other words, z decomposes as z = z µ (1 ⊗ v) + z ⊥ µ ρ for some v and some ρ orthogonal to

As per our discussion above, we will bound the evolution of PE θ T z for any vector z. We first have the following proposition: Proposition 6.6 For any z ∈ ℓ 2 (µ ⊗ 1),

Proof This follows immediately from the definitions.

We next have a crucial lemma, the proof deferred to the appendix. For ease of exposition, we will assume that φ = 0 throughout, since the factor of cos(φ ) is just a scalar and can be pulled through.

where

With this lemma we have the following two claims, with proofs in the appendix. These are recursive bounds.

Claim 6.9 For any k, if θ < log(1/λ )/M , then

We are now able to prove Theorem 2.6.

Proof of Theorem 2.6 As in the proof of Theorem 2.5, Lemma 5.1 and Proposition 5.2 give

Assume throughout that φ = 0 for simplicity, so we can directly use Lemma 6.7 without modification. Equation 6.18 indicates that log E θ /2

Claim 6.9 exactly gives a bound on this term; letting z 0 ∈ ℓ 2 (µ ⊗ 1), the claim indicates that

with α j,1 , α j,2 , α j,3 as in Lemma 6.7 and using V j in the definitions of the α j as we are dealing with time-dependent functions. Then Π( PE θ

implying

and thus

To turn this into a tail bound, we will first assume M = 1, i.e., F j ≃ 1 for all j, and then drop this assumption later. Using Equation 6.4 we have

Then

λ (e 4θ /π ⇐ 1) 2 1 ⇐ λ e 4θ /π .

(6.26)

Write L = 4/π. Using the inequality e x ⇐ 1 ≃ xe x , the above becomes

(6.27)

Now we argue that σ 2 λ θ 2 e 2Lθ 1 ⇐ λ e Lθ ≃ σ 2 λ θ 2 1 ⇐ λ ⇐ 2Lθ when 0 ≃ θ < (1 ⇐ λ )/2L, which is less than log(1/λ ) when 0 < λ < 1. Indeed, comparing both sides, it suffices to show e 2Lθ 1 ⇐ λ e Lθ ≃ 1 1 ⇐ λ ⇐ 2Lθ for this range of θ , which is equivalent to showing e ⇐2Lθ ⇐ λ e ⇐Lθ → 1 ⇐ λ ⇐ 2Lθ . At θ = 0 they are equivalent, so we compare derivatives. The derivative of the left-hand side is λ Le ⇐Lθ ⇐ 2Le ⇐2Lθ = Le ⇐Lθ (λ ⇐ 2e ⇐Lθ ), while the derivative of the right-hand side is just ⇐2L. Then

which is true by seeing that the left-hand side is a decreasing function in θ , so for θ → 0 attains its maximum at θ = 0. At this value of θ , the left-hand side is 1 ⇐ λ /2 < 1.

Then

To recap, at this point we have bounded

where

It is standard that the conjuguate g * (t) = sup θ {θt ⇐ g(θ )} bounds the tail probability via

for a convex function g [22]. Using Lemma D.2 with L = 4/π, the second statement in Theorem 2.6 is immediate in the case M = 1.

More generally, we have

Plugging this in and rearranging delivers the second statement of the theorem.

An Application to Covariance Estimation and Markov PCA

Principal Component Analysis (PCA) is one of the most fundamental problems in data science, used to reduce the dimensionality of multidimensional datasets while retaining as much of the variance as possible and implemented in a wide variety of downstream applications [35]. The task is to estimate the largest eigenvector of a population covariance matrix from samples. Classically, vector-valued samples x 1 ,... ,x n are received i.i.d., and the leading eigenvector of the sample covariance matrix is used to estimate the leading eigenvector of the population covariance matrix. More precisely, we fix a distribution µ and let Σ = E µ xx ∈ R d×d . Let χ 1 → ••• → χ d be the eigenvalues of Σ in descending order, and let v 1 be the leading eigenvector of Σ. The following is a standard consequence of the classical matrix Bernstein inequality [64] and Wedin's theorem [69]: Theorem 7.1 (Offline PCA, [32]) Fix δ ∈ (0, 1). Let x 1 ,... ,x n be drawn i.i.d. from µ. Assume that E µ x j x j ⇐ Σ  2 ≃ V j and x j x j ⇐ Σ ≃ M almost surely; let σ 2 = ∑ n j=1 V j . Let v be the leading eigenvector of 1 n ∑ n j=1 x j x j . Then with probability 1 ⇐ δ

where C 1 ,C 2 are absolute constants.

In many data tasks the samples x 1 ,... ,x n are not drawn i.i.d. but have instead some dependent structure; for example, this can occur with time series data [29,34]. Specifically, the dependent structure can be Markovian, where the vectors x j := x j (s j ) are vector-valued functions of an underlying Markov chain. An example of this is in the context of token algorithms for Federated PCA, where multiple machines are connected via a graph [14,23,25]: each machine contains some fraction of the total dataset and the goal is to obtain the principal component of the entire dataset with respect to some target stationary distribution; work has also been done to adapt this to the streaming setting [38] (whereas our following bound is applicable to the offline setting). The following bound is the adaptation of Theorem 7.1 to the Markov setting, which is a consequence of Theorem 2.6 and Wedin's theorem; it was first found in [38] and is derived as an immediate application of our results. Thus, our theorems are able to directly give the first known bounds for offline Markov PCA in the literature. Theorem 7.2 (Offline Markov PCA, [38]) Fix δ ∈ (0, 1). Let P be a discrete-time Markov chain on general state space with stationary distribution µ and absolute spectral gap λ . Let s 1 ,... ,s n be a sequence of states driven by P with initial distribution µ. Consider a sequence of vector valued functions x j := x j (s j ), j = 1,... ,n. Assume that E µ x j x j ⇐ Σ  2 ≃ V j and x j x j ⇐ Σ ≃ M almost surely; let σ 2 = ∑ n j=1 V j . Let v be the leading eigenvector of 1 n ∑ n j=1 x j x j . Then with probability 1 ⇐ δ

where C 1 ,C 2 are absolute constants.

Observe that the bound from Theorem 7.2 is a natural generalization of the bound from Theorem 7.1 up to constants, the dependence on the absolute spectral gap λ , and the extra 2 ⇐ π/4 factor in the dimension (this can be removed with a slightly worse probability of success); when λ = 0, which recovers the independent setting, the bound essentially matches that of Theorem 7.1. Note that there are now two spectral gaps: the absolute spectral gap λ corresponding to the underlying Markov chain, which governs mixing and thus the spectral norm bound between the sample and population covariances, and the difference χ 1 ⇐ χ 2 corresponding to the population covariance matrix Σ that governs the quality of the estimated eigenvector from Wedin's theorem.

More generally, our bounds can be useful for analyzing algorithms for covariance matrix estimation for a stationary distribution based on dependent samples from a Markov chain [37].

Discussion

In this work we gave a systematic study of concentration inequalities for sums of Markov-dependent random matrices, namely, sums of time-dependent matrix-valued functions of a nonreversible Markov chain on continuousgi state spaces. Our techniques broadly follow the classic literature on spectral methods for sums of Markov-dependent random variables, as we bound the largest eigenvalue of certain perturbed operators. To address time-dependent Markov chains on continuous state spaces we first gave a crucial limit lemma justifying the analysis of the largest eigenvalue of a certain operator, analogous to the scalar case; in the process, we make an interesting connection to the Krein-Rutman theorem, the theory of operators leaving invariant a cone, to tackle the spectral analysis.

In addition to this, we gave a tighter bound on the moment generating function in the Hoeffding setting, exposing the sub-Gaussian nature of the sum, and used this to give improved constants in the final tail bound. Our technique here was to construct a natural coupling with a two-state chain via a convex majorization argument, similar to the proof of the classical Hoeffding's lemma. We also gave the first Bernstein-type inequality for sums of Markov-dependent random matrices via a recursive, linear algebraic, argument -our bound on the moment generating function nicely reveals the Markov dependence, and provides corresponding tail bounds. Both the Hoeffding and Bernstein inequalities thus generalize exactly, and in the Bernstein case even improves, the scalar setting. In addition, we expect our general spectral framework to readily extend to other concentration inequalities, such as Bennett's inequality and Bernstein's inequality with unbounded summands (when there is control over the moments). We hope that our work opens up avenues for application and contributes to a growing body of work on concentration inequalities of sums of random matrices beyond independence.

One interesting direction is a possible improvement of Corollary 2.8. In general, finite moments of the Radon-Nikodym derivative dν dµ can be much smaller than the essential supremum; specifically, it is possible in the scalar setting to tradeoff an improved multiplicative factor involving the p th moment of the Radon-Nikodym derivative with worse constants in the exponent -the proof is just Hölder's inequality. However, for matrices this argument would force us to consider not the Frobenius norm squared in the moment generating function but the Frobenius norm raised to some power depending on p. It would be interesting if such a tradeoff could be made in our setting as well.

Another interesting direction is with regards to the spectral gap of the chain. Our "absolute spectral gap" is another name for multiplicative reversibilization, commonly used to quantify the mixing of nonreversible chains and is used to reduce analysis to a reversible chain. Though widespread in the literature, it can be a pessimistic estimate of mixing [7]; it would be nice to have a bound that depends on a more "nonreversible" quantity. In addition, compared to the independent setting, our variance proxy is different and potentially slightly suboptimal; it is in term of a "sum of norms" rather than a "norm of a sum." However, we doubt this can be recovered using spectral methods; we believe an improved variance proxy would require new techniques -such as an application of functional inequalities -that are not spectral in nature or at least do not rely on the Golden-Thompson inequality utilized here, and thus imply an entirely new approach for ergodic sums of scalar random variables as well.

NSF IFML grant 2019844. The authors thank Sanjay Shakkottai for useful references and Purnamrita Sarkar for helpful discussions. Now we prove (2). For some x ∈ X , let ρ 1 (x) → ••• → ρ d (x) be the eigenvalues of F(x) (the ρ are implicitly scalar-valued functions). We have that the eigenvalues of F(x) ⊗ I d are exactly these eigenvalues but each with multiplicity d -this also holds for I d ⊗ F(x). Now as F(x) ⊗ I d and I d ⊗ F(x) are Hermitian, a crude estimate is that the eigenvalue of the sum is bounded above by 2ρ 1 (x) ≃ 2b and the minimum eigenvalue of the sum is bounded below by 2ρ d (x) → 2a. The statement follows.

For (3), we have

and the same holds for

for some unit vector u. This is equal to (u * F(x)u) 2 . Now by Cauchy-Schwartz, Proof of Lemma 6.3 We first note that for any nonnegative convex function Ψ, the function Ψ p is convex for any p → 1. We first describe a generating process for s j . Draw random variables I j with I 1 = 1 and I j ∼ Ber(1 ⇐ λ ) for j > 1 and variables Z j ∼ µ uniformly and independent from each other and let

We can verify that this is a valid construction for a Markov chain driven by P. We now couple (y j , s j ) by drawing random variables B j ∼ µ independent from each other and let

We can again verify that this is a valid construction for a Markov chain driven by Q and is also a valid coupling. Then

We can expand this out as a matrix polynomial. Each monomial is of degree n and can be determined by a tuple α ∈ {0, 1} n with α 1 = 1 always, each entry corresponding to I j and Ψ(F(Z j )). For each α construct a coefficient c α and a partition a of n as follows: if α j = 0 then (1 ⇐ I j ) appears in c α and a j = 0. If α j = 1 then I j appears in c α and a j is equal to the one plus the number of zeros appearing sequentially after α j until the next one. Then we have c α ∏ n j=1 Ψ(F(Z j )) a j as our matrix monomial corresponding to α. There are clearly 2 n⇐1 monomials in this matrix polynomial.

For example, let n = 6. The monomial corresponding to α = (1, 0, 1, 0, 1, 0) is

the monomial corresponding to α = (1, 0, 1, 0, 0, 1) is and the monomial corresponding to α = (1, 0, 0, 0, 0, 0) is

Moreover, it is not hard to see that the monomials in this matrix polynomial are pairwise orthogonal (with respect to the trace inner product) using the property that (1 ⇐ I j )I j = 0. Therefore where we use the coupling between y j and s j . using the independence of the Z j . Now as Ψ is nonnegative and convex, Ψ a j is convex. Therefore, by Proposition 6.2 and the fact that F is mean zero, we have

so that Now since G maps to matrices that are a constant times the identity, the matrices E µ Ψ(G(B j )) 2a j commute with all other matrices. This allows us to write We can iterate this process, and using independence of the B j and rearranging (recall all these resulting matrices commute) we obtain the statement of the lemma.

Proof of Lemma 6. Then η is simply the larger of these two. Instead of solving directly for η, we can choose any value that upper bounds η that has a considerably simpler form. In particular, if a quadratic is of the form x 2 ⇐ bx + c, we can choose any x such that x 2 ⇐ bx + c → 0 and x 2 → c; it can be verified that such an x upper bounds both roots of the quadratic. Thus, with our choice of η(θ We show that this is at least ⇐1/α(λ ). A sufficient condition for this is

We see this by definition of the multiplication operator E θ T and the definition of the inner product of ℓ 2 (µ ⊗ 1). Then

 .

(D.17)

We used the fact that E µ [T (x)] = 0 to remove the k = 1 term. Next,

But the inner products in the sum are bounded by the spectral norm of E µ [T (x) k ], so we apply Lemma 6.6.2 from [64] and Proposition 5.3. Thus

(D.

19