Abstract

We give lower and upper bounds on both the Lyapunov exponent and
generalised Lyapunov exponents for the random product of positive and
negative shear matrices. These types of random products arise in
applications such as fluid stirring devices. The bounds,
obtained by considering invariant cones in tangent space, give
excellent accuracy compared to standard and general bounds, and are
increasingly accurate with increasing shear. Bounds on generalised
exponents are useful for testing numerical methods, since these
exponents are difficult to compute in practice.

1 Introduction

Random matrix products have applications in many disciplines such as statistical and nuclear physics [1], population dynamics [2] and quantum mechanics [3]. Their rigorous study began over sixty years ago, when Bellman [4] studied the asymptotic behaviour of products of random matrices with strictly positive entries, corresponding to a weak law of large numbers. The seminal work of Furstenberg & Kesten [5] and Furstenberg [6] strengthens this to a strong law for more general matrices. Oseledec [7] extended this further to matrix cocycles of dynamical systems.

Here we consider the random product with N terms of the two matrices {A1,A2},

MN=N∏k=1Aik,ik∈{1,2},

(1)

where the ik are i.i.d. and the two index values have equal
probability 1/2. It will often be convenient to write

A=A1andB=A2.

(2)

The Lyapunov exponent is defined by

λ=limN→∞1NElog∥MN∥

(3)

where ∥⋅∥ is some matrix norm. The Furstenberg–Kesten
theorem [5, 6] states
that the limit (3) exists, and is positive under fairly weak
assumptions on the Ai, satisfied by the matrices we will be using. The Lyapunov exponent can be equivalently defined using a vector norm rather than a matrix norm as

λ=limN→∞1NElog∥XN∥,XN=MNX0,

(4)

for an arbitrary vector X0. In this paper we use the definition given by (4), and the choice of initial X0 will be clear.

There is a paucity of exact results concerning Lyapunov exponents for random
matrices, as famously lamented by Kingman [8, p. 897]. One well-known upper bound (described by [9] as “the most popular upper bound in the literature”) is easily derived from the submultiplicativity of ∥⋅∥. For two matrices chosen with equal probability, let

Ek=1kElog∥C∥,

(5)

with C∈Ak, where Ak is the set of all 2k products of matrices of length k. The numbers Ek converge monotonically to λ from above as k→∞ for any choice of matrix norm, although according to [9] the Euclidean norm is usual. In [10] the bound is described as “easy, if not efficient”, since the number of matrix product calculations required increases exponentially with k.

Further progress in this direction has tended to be either for specific simple cases, or algorithmic procedures leading to (sometimes very accurate) approximations. For example, [11] and [12] discuss cases where the Lyapunov exponent can be computed exactly, in particular when matrices can be grouped in commuting blocks. Chassain et al[13] establish the distribution for the matrix product, in terms of a continued fraction, in the case that the matrices are 2×2 shear matrices, but observe that even for these simple matrices, the Lyapunov exponent is still unobtainable. A similar approach allowed Viswanath [14] to give a formula for the exponent in the case of matrices which give rise to a random Fibonacci sequence (this was extended by Janvresse et al[15]). An exact expression for λ as the sum of a convergent series in the case for which one matrix is singular was given by Lima & Rahibe [16]. Analytic expressions for λ have also been obtained for large, sparse matrices [17], and for classes of 2×2 matrices in terms of explicitly expressed probability distributions [18, 19]. Pollicott
[20] recently gave a cycle expansion formula that
allows a very accurate computation for a class of matrices.
Protasov & Jungers [9] obtain an efficient algorithm for Lyapunov exponent bounds
using invariant cones for matrices with non-negative entries, concentrating on
generality and efficiency (they test their algorithm on examples up to dimension 60).

The difficulty of calculating Lyapunov exponents for products of random matrices can be seen in the variety of approaches. All the above results and algorithms, with the exception of those for random Fibonacci sequences, hold only for matrices with non-negative entries. Moreover, analytic expressions tend to be given in terms of probability distributions, continued fractions, or convergent series. In the present work we aim to provide rigorous and explicit upper and lower bounds for λ for the same non-singular matrices as studied in [13].

The matrices in question are shear matrices, which are of particular interest in several fluid mixing problems [21, 22, 23, 24] and devices known as ‘taffy pullers’ [25, 26, 27]. The principle of mixing by chaotic advection can be summarised as repeated stretching in transverse directions[28, 29]. Many industrial mixing devices are designed on this basis, with the most fundamental model being periodic application of transverse shear matrices [21, 22, 23, 24]. Our work applies to devices where the mixing action is random.

For the problems of passive scalar decay and random taffy pullers, knowledge of the Lyapunov exponent is insufficient [30, 31, 32]. We require more refined information via the growth rate of the qth moment of the matrix product norm. In particular, define the generalised Lyapunov exponents[33, 1] as

ℓmat(q)=limN→∞1NlogE∥MN∥q.

(6)

Again this can be restated using a vector norm:

ℓ(q)=limN→∞1NlogE∥XN∥q,

(7)

where XN is as defined in (4), but here we must observe that these definitions are not equivalent, particularly for q<0. In the present paper we will use the vector norm definition (7).

2 Rigorous bounds for Lyapunov exponents

We derive rigorous and explicit bounds for Lyapunov exponents and generalised Lyapunov exponents by reformulating the problem, grouping the matrices together. Assume without loss of generality that the first
matrix in the product (1) is A1=A. By grouping A’s and B’s together into J blocks the random product (1) can be
written

MNJ=J∏j=1AajBbj,aj+bj=nj,J∑j=1nj=NJ,

(8)

with 1≤ai,bi<ni, so ni≥2. Now it is the ai and bi
that are the i.i.d. random variables, with identical probability
distribution P(x)=2−x, x≥1. Hence the length of each block is governed by the joint
distribution P(a,b)=P(a)P(b)=2−(a+b). We have the expected
values Ea=Eb=2, so En=4.

Let us now take the specific matrices

A=(10α1),B=(1β01),Kab\raisebox0.2967pt:=AaBb=(1bβaα1+aαbβ).

(9)

We consider first the case α,β>0, for which Kab is positive definite, and hyperbolic ∀a,b≥1. Although our technique holds for all positive α,β, we state our results for α,β≥1. This is partly due to ease of exposition, but also because in many applications α and β would be assumed to be integers, so that a map induced by Kab is continuous on the 2-torus. In particular, the algebraically simplest case α=β=1 corresponds to the generators of the 3-braid group seen in many studies of topological mixing [25, 26, 27]. We then allow negative entries; in particular we consider α<0<β (note that α>0>β is essentially similar, while α<0,β<0 is no different from the positive α,β case). Now hyperbolicity is only guaranteed when the product |αβ| is sufficiently large, and we require this property to obtain our results. We gain different bounds by considering different vector norms, a valid approach since the limit in (4) is independent of the choice of norm. In particular, we will consider the L1, L2 and L∞ norms. Which norm produces the most accurate bound depends on α and β. This is easily discerned by computation.

2.1 Lyapunov exponents

Our theorems are stated in terms of infinite sums of products of an exponentially decreasing term and a choice of (logarithm of) increasing algebraic function, and so all obviously converge.

Theorem 1.

The Lyapunov exponent λ(α,β) for the product MN for α,β≥1 satisfies

maxk∈{1,2,∞}Lk(α,β)≤4λ(α,β)≤mink∈{1,2,∞}Uk(α,β)

where

Lk(α,β)

=

∞∑a,b=12−a−blogϕk(a,b,α,β)

Uk(α,β)

=

∞∑a,b=12−a−blogψk(a,b,α,β),

and

ϕ1(a,b,α,β)

=

1+α1+α(a+bβ+aαbβ)

ϕ2(a,b,α,β)

=

min⎧⎨⎩((1+aαbβ)2+b2β2)1/2(11+α2(α2(1+a+aαbβ)2+(1+αbβ)2))1/2

ϕ∞(a,b,α,β)

=

1+aαbβ

ψ1(a,b,α,β)

=

1+bβ+aαbβ

ψ2(a,b,α,β)

=

(12(2+Caαbβ+√Caαbβ(Caαbβ+4)))1/2

ψ∞(a,b,α,β)

=

1+a+aαbβ

where Caαbβ=(aα+bβ)2+(aαbβ)2.

Losing a little sharpness, the L∞-norm bounds provide a pair of simpler expressions with no infinite sums, stated in:

Corollary 1.

The Lyapunov exponent λ(α,β) for the product MN for α,β≥1 satisfies

κ+logαβ≤4λ≤κ+log(√αβ+1/√αβ)+12log(1+αβ),

where

κ=∞∑a,b=12−a−blogab≈1.0157….

Theorem 1 is obtained by considering a cone in tangent space which is invariant for all a and b. We can improve these estimates by recognising that a smaller cone can be used for certain iterates of the map. In particular, we use the fact that since ai and bi are independent geometric distributions, P(a=b)=P(a>b)=P(b>a) to give

An immediate observation is that since all the functions ϕ,ψ,^ϕ(m)k and ^ψ(m)k are greater than 1 for all a,b≥1, α,β≥0, and since ∑∞a,b=12−a−b=1, the bounds for ℓ(q,α,β) grow from 0 for positive q and decay from zero for negative q. This apparently contradicts Proposition 2 of [34], which states that there is always a minimum in the curve for ℓ(q), and in particular states that ℓ(−2)=0 if the 2-dimensional matrices in question have determinant 1. The existence of the invariant cone for these shear matrices guarantees that a vector is expanded at every application of A or B, which forces ℓ(q) to be monotonic. In [34], the assumption is made that the linear operator corresponding to the generalised Lyapunov exponent has the same spectrum as its adjoint, a property precluded by the invariant cone. The fact ϕ,ψ,^ϕ(m)k and ^ψ(m)k≥1 is also the reason that ϕ and ψ exchange roles in upper and lower bounds for positive and negative q.

When q is a positive integer we can evaluate this easily by expanding the power q. Since

∞∑a,b=12−a−banbm=(∞∑a=12−aan)(∞∑b=12−bbm)

we require values of the polylogarithm Li−q(12), defined by

Lis(z)=∞∑k=1zkks.

For integer q=−s we have special values

∞∑a=12−aan=1,2,6,26,150,1082,9366,… for n=0,1,2,3,4,5,6,…

and so the L∞ norm, for example, gives

Corollary 2.

Generalised Lyapunov exponents in the case α=β=1 are bounded by:

14log5

≤

ℓ(1,1,1)≤14log7

14log45

≤

ℓ(2,1,1)≤14log79

14log797

≤

ℓ(3,1,1)≤14log1543

14log25437

≤

ℓ(4,1,1)≤14log50531

14log1290365

≤

ℓ(5,1,1)≤14log2578567.

Theorem 3 also allows explicit estimates on topological entropy for the random matrix product, given by the generalised Lyapunov exponent with q=1.

Corollary 3.

The topological entropy ℓ(1,α,β) in the case α,β≥1 is bounded by

log(1+4αβ)≤4ℓ(1,α,β)≤log(3+4αβ).

2.3 Matrices with negative entries

The case where the direction of one of the shears is reversed (that is, allowing negative entries in the matrix) can be tackled in an almost identical manner, with one important condition. Taking α<0<β (the case α>0>β is essentially identical), the matrix K11=AB is hyperbolic only when the product |αβ|>4. In the following, for simplicity, we will assume α<−2, β>2 to achieve this111In fact, we might assume that α=−β. Otherwise we change coordinates, as in [35] to (x,√|α/β|y). But here we retain α≠−β to show explicitly the dependence of the bounds on choosing unequal strengths of twists..

Theorem 4.

The Lyapunov exponent λ(α,β) for the product MN in the case α<−2,β>2 satisfies

maxk∈{1,2,∞}~Lk(α,β)≤4λ(α,β)≤mink∈{1,2,∞}~Uk(α,β)

where

~Lk(α,β)

=

∞∑a,b=12−a−blog~ϕk(a,b,α,β)

~Uk(α,β)

=

∞∑a,b=12−a−blog~ψk(a,b,α,β),

and

~ϕ1(a,b,α,β)

=

11−Γ(bβ+Γ−aαbβ−1−aαΓ)

~ϕ2(a,b,α,β)

=

(11+Γ2((Γ+bβ)2+(1+aαΓ+aαbβ)2))1/2

~ϕ∞(a,b,α,β)

=

−aαbβ−aαΓ−1

~ψ1(a,b,α,β)

=

bβ−aαbβ−1

~ψ2(a,b,α,β)

=

((−aαbβ−1)2+b2β2)1/2

~ψ∞(a,b,α,β)

=

−aαbβ−1

where

Γ=−β2+√(β2)2+βα.

Again we can straightforwardly improve on the lower bounds by considering separately the cases when either, or both, of a and b are equal to 1.

Theorem 5.

The Lyapunov exponent λ(α,β) for the product MN in the case α<−2,β>2 satisfies

maxk∈{1,2,∞}^~Lk(α,β)≤4λ(α,β)≤^~U∞(α,β)

where

^~Lk(α,β)

=

∞∑a,b=12−a−blog142∑ma,mb=1^~ϕ(ma,mb)k(a,b,α,β)

^~U∞(α,β)

=

∞∑a,b=12−a−blog144∑m=1^~ψ(m)∞(a,b,α,β)

and

^~ϕ(ma,mb)1(a,b,α,β)

=

11−Γma,mb(bβ+Γma,mb−aαbβ−1−aαΓma,mb)

^~ϕ(ma,mb)2(a,b,α,β)

=

(11+Γ2ma,mb((Γma,mb+bβ)2+(1+aαΓma,mb+aαbβ)2))1/2

^~ϕ(ma,mb)∞(a,b,α,β)

=

−aαbβ−aαΓma,mb−1

^~ψ(1)∞(a,b,α,β)

=

−aαbβ−1−aαβ/(1+αβ))

^~ψ(2)∞(a,b,α,β)

=

−aαbβ−1−a

^~ψ(3)∞(a,b,α,β)

=

^~ψ(4)∞(a,b,α,β)=~ψ∞(a,b,α,β)

with

Γma,mb=Γ+mbβmaαΓ+mambαβ+1

(10)

for ma,mb=1,2. Note that Γ1,1=L.

We could also write improved upper estimates using L1 and L2 norms, but since these produce worse bounds than the L∞ norm in all cases we study here, we do not give these explicitly.

As before we can use the same estimates to give explicit bounds for generalised Lyapunov exponents.

3 Accuracy of the bounds

3.1 Lyapunov exponents

For α=β=1 we have bounds on the Lyapunov exponent given in table 1. The lowest upper bound (U2) and largest lower bound (^L1) differ by about 2.5%. The true value (from explicit calculation via a standard algorithm [36]) is 0.39625…, so the upper bound here is rather tighter than the lower.

Invariant cone

Smaller cones

Norm

Lower bound

Upper bound

Improved bounds

L1

L1(1,1) = 0.36886

U1(1,1) = 0.43835

^L1(1,1) = 0.38561

L2

L2(1,1) = 0.36347

U2(1,1) = 0.40277

^L2(1,1) = 0.36864

L∞

L∞(1,1) = 0.34613

U∞(1,1) = 0.43835

^U∞(1,1) = 0.41350

Table 1: Bounds to five significant figures for the maximal Lyapunov exponent for the matrix product (1) in the case α=β=1, where the true value (from numerical computation) is 0.39625… .

Figure 1 shows the accuracy of each bound from Theorem 1 increasing as α increases, with α=β∈[1,10]. In figure 0(a) the bounds are all so close to the true value of λ that the details of the graph are difficult to resolve. It is clear however, that the standard bound given by (5) (plotted in cyan) is a worse upper bound than all others in the figure, despite being calculated from the expected value of matrix products of length 212, and decreases in accuracy for this fixed k for increasing α.

Figure 0(b) shows the difference between the bounds and the true (numerically calculated) Lyapunov exponent. In this and other figures we colour bounds originating from L1-, L2- and L∞-norms green, red and blue respectively. It shows that for increasing α, upper bounds (solid lines) appear tighter than lower bounds (dashed lines). In black is shown the upper and lower bounds given in Corollary 1, which are typically worse than those of Theorem 1, but have the advantage of being explicit, finite expressions rather than infinite sums.

Figure 0(c) shows the envelope formed from the difference between upper and lower bounds derived from each norm, which does not require the explicit numerical calculation of the Lyapunov exponent to compute. To this figure we add, in figure 0(d), the corresponding bounds from Theorem 2 which improve on Theorem 1 by considering the expected relation between the random variables ai and bi. In black is the envelope formed from taking the minimum upper bound, and maximum lower bound for each value of α. This improves on all bounds produced from a single norm.

(a)Numerical estimate of the Lyapunov exponent obtained by random
multiplication of the matrices (9), with bounds as given in Theorem 1. Only the standard bound is appreciably far from the true value.

(b)Errors in bounds from numerically calculated value.

(c)Difference between upper and lower bounds.

(d)Envelope of bounds when improved cone is included.

Figure 1: Four different views of the accuracy of the upper and lower bounds for the Lyapunov exponent for the random matrix product (1) for α=β∈[1,10]. In each case, bounds obtained from different norms are shown, in particular L1 (green), L2 (red) and L∞ (blue) are shown, with bounds from the global cone shown dashed, and the improved cone shown as a solid line. When shown, the standard bound is cyan.

The increase of accuracy of the bounds with increasing strength of shear can be understood geometrically, as the size of the cone narrows with increasing shear, and dynamically, as the approach to the unstable eigenvalue is more rapid for matrices whose largest eigenvalue is greater. In figures 1 and 3 it is clear that the upper and lower bounds approach the same curve for large α=β. A simple calculation (using the L∞-norm) gives

4λ

≥

∞∑a,b=12−a−blog(1+aαbβ)

≥

∞∑a,b=12−a−blog(aαbβ)

=

∞∑a,b=12−a−blog(ab)+∞∑a,b=12−a−blogαβ

≥

κ+logαβ.

Meanwhile, for large α,β we also have

4λ

≤

∞∑a,b=12−a−blog(1+a+aαbβ)

∼

∞∑a,b=12−a−blog(aαbβ)

∼

κ+logαβ,

and this indeed appears to be the asymptotic limit in the graphs shown for α=β→∞.

3.2 Generalised Lyapunov exponents

In figure 2 we show bounds for generalised Lyapunov exponents for α=β=1. Equivalent figures are increasingly accurate with increasing α,β. Figure 1(a) confirms that for this choice of matrices we do not have ℓ(−2)=0, and that there is no minimum in the curve of ℓ(q). Green, red and blue lines again show bounds originating from L1-, L2- and L∞-norms respectively, with the explicit expressions from corollary 2 shown as black circles. Figure 1(b) shows the envelope of the difference between upper and lower bounds, for each norm, and, in black, the envelope of best combined bounds.

(a)The bounds confirm that the curve of generalised Lyapunov exponents has no minimum at q=−2.

(b)Difference between upper and lower bounds. Dashed lines represent bounds from Theorem 1, while solid lines are those from 2. The black line represents the minimal envelope over all norms, and thus our best bounds.

Figure 2: Bounds for the generalised Lyapunov exponent for the matrix product 9 with α=β=1. As before, estimates originating from L1-, L2- and L∞-norms are given in green, red and blue respectively. For integer values of q>0, values from Corollary 2 are given as black circles.

3.3 Negative shears

Figure 3 shows the bounds for the case α<−2, β>2. In these figures we set α=−β, and observe that again, the increasing hyperbolicity from increasing |α| results in improving bounds. In this case the L∞-norm always gives the minimal envelope, as seen in figure 2(b). Generalised Lyapunov exponents for α=−3, β=3 are shown in figure 4.

(b)Envelope of bounds from Theorem 4, shown dashed, and from Theorem 5, shown as solid lines.

Figure 3: Bounds for the negative entry case, in which α<−2,β>2. In this case the L∞-norm always produces the most accurate bounds.

(a)The curve of generalised Lyapunov exponents for α<0, β>0.

(b)Difference between upper and lower bounds. Dashed lines represent bounds from Theorem 4, while solid lines are those from 5. The black line represents the minimal envelope over all norms, and thus our best bounds.

Figure 4: Bounds for the generalised Lyapunov exponent for the matrix product 9 with α=−3, β=3. As before, estimates originating from L1-, L2- and L∞-norms are given in green, red and blue respectively.

4 Bounds on the growth of matrix norms

4.1 Invariant cones

We obtain bounds for Lyapunov exponents by computing explicit bounds for the norm of tangent vectors under the action of Kab. The expression (3) is independent of the matrix norm used, and we give bounds derived from three standard norms.

(a)The α>0 case. Here we show the cone C+ for α>1, where it lies inside the line u/v=1. The expanding eigenvector v+ also lies inside the cone C+, and so the orthogonal contracting eigenvector v− lies outside C+, and hence Lemma 3.

(b)The α<0 case. For αβ<−4, when the matrix Kab is hyperbolic, the cone C− lies inside the line u/v=L∈(−1,0). Both eigenvectors v+ and v− both lie outside the invariant cone for all α<−2, β<−2, which produces Lemma 9.

Figure 5: The invariant cones C+ and C− in both the α>0 and α<0 cases, with expanding and contracting eigenvectors, v+ and v− respectively, of the matrix KTabKab.

Lemma 1.

The cone C+={(u,v):0≤u/v≤1/α} (shown in figure 4(a)) is invariant under Kab for all a,b≥1, and is the smallest such cone.

Proof.

The vector

(u′v′)=(1bβaα1+aαbβ)(uv)

is such that

u′v′=u+bβvaαu+(1+aαbβ)v≥0,

clearly, and since au+abv≥u+bv for a≥1, we also have u′/v′≤1/α. Setting (a,b)=(1,∞) gives u′/v′=1/α, while setting (a,b)=(∞,1) gives u′/v′=0, showing that the cone cannot be made smaller.

∎

Throughout this section, we will consider a vector X=(u,v)T∈C+, and assume without loss of generality that u,v≥1 (the calculations for u,v,<0 are entirely analogous). We will consider the norm of the vector KabX, given by