Lyapunov Exponents for the Random Product of Two Shears

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.


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 {A 1 , A 2 }, where the i k are i.i.d. and the two index values have equal probability 1/2.It will often be convenient to write A = A 1 and B = A 2 . ( The Lyapunov exponent is defined by 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 A i , 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 for an arbitrary vector X 0 .In this paper we use the definition given by (4), and the choice of initial X 0 will be clear.

arXiv:1706.03398v1 [math.DS] 11 Jun 2017
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 with C ∈ A k , where A k is the set of all 2 k products of matrices of length k.The numbers E k 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 [1,33] as Again this can be restated using a vector norm: where X N 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).

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 A 1 = A. By grouping A's and B's together into J blocks the random product (1) can be written with 1 ≤ a i , b i < n i , so n i ≥ 2. Now it is the a i and b i 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 We consider first the case α, β > 0, for which K ab 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 K ab 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 L 1 , L 2 and L ∞ norms.Which norm produces the most accurate bound depends on α and β.This is easily discerned by computation.

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 M N for α, β ≥ 1 satisfies where and Losing a little sharpness, the L ∞ -norm bounds provide a pair of simpler expressions with no infinite sums, stated in: where 2 −a−b log ab ≈ 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 a i and b i are independent geometric distributions, where

Generalised Lyapunov exponents
We can use the functions defined above to bound the generalised Lyapunov exponents for each q: Theorem 3. We have, for α, β ≥ 1, and the more accurate expressions k and ψ(m) k defined as above.
An immediate observation is that since all the functions φ, ψ, k and ψ(m) k are greater than 1 for all a, b ≥ 1, α, β ≥ 0, and since ∞ a,b=1 2 −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=1 we require values of the polylogarithm Li −q ( 1 2 ), defined by For integer q = −s we have special values Theorem 3 also allows explicit estimates on topological entropy for the random matrix product, given by the generalised Lyapunov exponent with q = 1.

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 K 11 = AB is hyperbolic only when the product |αβ| > 4. In the following, for simplicity, we will assume α < −2, β > 2 to achieve this ‡. where 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. where In 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.
We could also write improved upper estimates using L 1 and L 2 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.
Theorem 6.We have, for α < −2 and β > 2, and the more accurate, but more complicated expressions k defined as above.

Lyapunov exponents
For α = β = 1 we have bounds on the Lyapunov exponent given in table 1.The lowest upper bound (U 2 ) 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.Figure 1 shows the accuracy of each bound from Theorem 1 increasing as α increases, with α = β ∈ [1,10].In figure 1a 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 2 12 , and decreases in accuracy for this fixed k for increasing α.
Figure 1b shows the difference between the bounds and the true (numerically calculated) Lyapunov exponent.In this and other figures we colour bounds originating from L 1 -, L 2 -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

Invariant cone Smaller cones Norm
Lower bound Upper bound Improved bounds 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... .typically worse than those of Theorem 1, but have the advantage of being explicit, finite expressions rather than infinite sums.
Figure 1c 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 1d, the corresponding bounds from Theorem 2 which improve on Theorem 1 by considering the expected relation between the random variables a i and b i .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.
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 Meanwhile, for large α, β we also have and this indeed appears to be the asymptotic limit in the graphs shown for α = β → ∞.

Generalised Lyapunov exponents
In figure 2 we show bounds for generalised Lyapunov exponents for α = β = 1.Equivalent figures are increasingly accurate with increasing α, β. Figure 2a 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 L 1 -, L 2 -and L ∞ -norms respectively, with the explicit expressions from corollary 2 shown as black circles.Figure 2b shows the envelope of the difference between upper and lower bounds, for each norm, and, in black, the envelope of best combined bounds.
Corollary 1 (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.
Corollary 1 (b) Errors in bounds from numerically calculated value.Upper bound -Lower bound (c) Difference between upper and lower bounds.Upper bound -Lower bound

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

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 3b.Generalised Lyapunov exponents for α = −3, β = 3 are shown in figure 4.

Invariant cones
We obtain bounds for Lyapunov exponents by computing explicit bounds for the norm of tangent vectors under the action of K ab .The expression (3) is independent of the matrix norm used, and we give bounds derived from three standard norms.
Lemma 1.The cone C + = {(u, v) : 0 ≤ u/v ≤ 1/α} (shown in figure 5a) is invariant under K ab for all a, b ≥ 1, and is the smallest such cone.(q) The bounds confirm that the curve of generalised Lyapunov exponents has no minimum at q = −2.As before, estimates originating from L 1 -, L 2 -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.Proof.The vector 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 curve of generalised Lyapunov exponents for α < 0, β > 0.
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.
Lemma 2. The norm K ab X 1 for a vector X ∈ C + satisfies (i) the lower bound Proof.For any X ∈ C + we have With α, β ≥ 1 this has no local minima or maxima in the cone C + .Thus the lower and upper bounds are attained at the boundaries of C + , given by (u, v) = ( 1 1+α , α 1+α ) and (u, v) = (0, 1) respectively.
For the L 2 -norm X 2 = √ u 2 + v 2 , the calculations are more involved, but the following holds: where Proof.The real 2 × 2 matrix K ab is non-singular, and so , where v + (v − ) is the eigenvector corresponding to e + (e − ), the larger (smaller) eigenvalue of K T ab K ab , by the definition of the spectral matrix norm (and by singular value decomposition).Moreover the value of K ab X 2 X 2 varies monotonically between these extremes.Since K T ab K ab is symmetric, v − and v + are orthogonal.
The eigenvector v + (r, s) satisfies Clearly r > 0, while s = C aαbβ −2a 2 α 2 + C aαbβ (C aαbβ + 4) > 2C aαbβ −2a 2 α 2 = 2b 2 β 2 +4aαbβ +2(aαbβ) 2 > 0, and so v + lies in the positive quadrant of tangent space.Moreover, we have s > , and so v + ∈ C + , giving the upper bound.The orthogonality of the eigenvectors then implies that v − / ∈ C + , and the lower bound is given by the minimum of the value of the spectral norm on the boundaries of C + .
Proof.For α ≥ 1 we have This takes minimum and maximum values at minimum and maximum values of u/v, respectively.For the cone C + these are given by 0 and 1/α, and the bounds follow immediately.
We can now use these bounds and the invariant cone to prove Theorem 1.
Proof of Theorem 1. Taking i.i.d.copies of K ab and defining . ., J, we have for an initial vector X 0 ∈ C + , By Lemma 1, each term in the product is a vector ∈ C + , and so is bounded according to Lemmas 2, 3 and 4. Hence since En = 4, so using the probability distribution P (a, b) we have as required.
To obtain Corollary 1 we select the algebraically simplest bounds (the L ∞ bounds), and evaluate the infinite sums where possible.

Cone improvement
In this section we improve on the lower bound by considering the relationship between two identical geometric distributions.
Lemma 5. When a and b are both i.i.d.geometric distributions with parameter 1/2, we have Proof.We have Then the remaining two equalities follow by symmetry.
Proof of Theorem 2. This follows the same argument as the proof of theorem 1, except that whenever it happens that a = b, or a > b, on the following iterate the vector X i is bounded according to Lemma 6.Since by Lemma 5 these conditions occur on average 1/3 of the time, the result follows.
The expanding eigenvalue e − has eigenvector (u, v) T with In the case α < −2 the minimal cone is bounded by this eigenvector when a = b = 1, so setting we have: Lemma 7. The cone C − = {(u, v) : Γ ≤ u/v ≤ 0} is invariant under K ab for all a, b ≥ 1, and for all α < −2, β > 2, and is the smallest such cone.
Proof.Without loss of generality we will take an initial vector (u, v) with u < 0, v > 0 (an initial vector in the opposite sector proceeds exactly analogously) in C − , so that −u < v and u > −v.Then we consider Now we use the fact that |aαu/v| > |u/v| to replace two of these terms while respecting the inequality: Rearranging then gives u /v ≥ Γ.This is the smallest such invariant cone, since setting (a, b) = (1, 1) gives u /v = Γ when u/v = Γ, and setting (a, b) = (∞, 1) gives u /v = 0.
As before the L ∞ −norm gives bounds easily: Proof.Since Γ > −1, for any X ∈ C − we have X ∞ = |v| and since C − is invariant under K ab , we have K ab X ∞ / X ∞ = |aα u v + 1 + aαbβ|, which takes minimum and maximum values at the boundaries (u, v) = (0, 1) and (u, v) = (Γ, 1) of the cone C − , and the bounds follow immediately.
For this invariant cone, the L 2 -norm • 2 cannot attain the spectral maximum, and the following holds: .
Proof.As in Lemma 3, we consider eigenvectors of K T ab K ab .For α < −2, β > 2, the expanding eigenvector v + = (r, s) still lies in the northeast-southwest quadrant, outside C − .But since v − = (−s, r), we have and so −s/r < −1, and hence v − also lies outside C − .Since the norm in question increases monotonically between the two extremes, neither of which lie in the cone, the lower and upper bounds are achieved at the minimum and maximum values (respectively) at the boundaries of C − .At the boundary given by (u, v) = (0, 1), we have 2 , while at the other boundary, given by (ii) the upper bound Proof.With the L 1 -norm we have K ab Proof of Theorem 4. This follows exactly the argument of Theorem 1, using Lemma 7 to guarantee an invariant cone, and using Lemmas 8, 9 and 10 to bound each term in the matrix product.

Cone improvement
In the α < 0 case we can make a significant improvement on the bounds given by Theorem 4 by recognising that the boundary u/v = Γ of the cone C − can only be achieved when a = b = 1, which occurs on average P (a = b = 1) = 1/4 of the time.Whenever a or b (or both) is greater than 1, we can assume a smaller cone for the following iterate.More precisely, since Inserting the boundaries of C − , given by u v = 0 and u v = Γ into this expression produces the required inequalities.The bounding functions are then obtained using analogous arguments to Lemmas 8, 9 and 10, with the new cone boundaries.
Proof of Theorem 5. Again this follows the same argument as the proof of theorem 1, using improved bounds given by Lemma 11, each of which applies 1/4 of the time, on average.

Generalised Lyapunov exponents
The expressions for (q) can be obtained in largely the same way, bounding the expansion of vectors at each application of matrix A or B.
Proof of Theorems 3 and 6.Using properties of expectation, and the independence of X i , we have and so since the a i , b i are i.i.d., we have ∞ a,b=1 2 −a−b ψ q .
Then from the definition of (q) given in (7) the results follow immediately.

Conclusions and discussion
In this paper we addressed the question of obtaining rigorous bounds for Lyapunov exponents, generalised Lyapunov exponents, and topological entropy for randomised mixing devices.The matrices under discussion are 2 × 2 shear matrices, but the same technique will work for any set of matrices that share an invariant cone.This notion is proved formally in [9], who give a rapid algorithm involving unconstrained minimisation problems.Here the optimisation is achieved analytically, giving explicit upper and lower bounds.We also obtain bounds in the novel case of shear matrices with negative entries.A pair of hyperbolic matrices sharing an invariant cone was shown to enjoy exponential decay of correlations in [37], where the rate of decay depends on the Lyapunov exponent, but here the Lyapunov exponent is simply bounded from below by global expansion and contraction rates in the invariant cone.The method in this paper could be adapted to tighten their lower bound, and provide an upper bound.
The assumption that the matrices A and B should be chosen with equal probability at each iterate can be relaxed.Altering these probabilities does not change the invariant cone, or the resulting bounds on vector norms; only the probability distribution P (a, b) = 2 −a−b is changed.For example, replacing the geometric probability distribution with a Bernoulli distribution gives P (a, b) = p a q b , and then Ea = q −1 , Eb = p −1 , and En = (pq) −1 .Similarly, one may choose from k matrices A i with probability p i at each iterate.The crucial element is that the expected length of a block should be computable.
Theorem 2 improves on 1 by involving the relative values of a and b in one block to shrink the cone for the next, in the three cases a = b, a < b and a > b.Similarly, the nine cases comprising the relative values of a and b in two consecutive blocks can increase the tightness of bounds in the following block.This procedure could be extended to further improve bounds, but the number of cases increases exponentiallyin k blocks there are 3 k combinations of relative values of a and b.Our original explicit bounds are appealing in their simplicity and accuracy.

Figure 1 :
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 L 1 (green), L 2 (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.

Figure 2 :
Figure2: Bounds for the generalised Lyapunov exponent for the matrix product 9 with α = β = 1.As before, estimates originating from L 1 -, L 2 -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.
Envelope of bounds from Theorem 4, shown dashed, and from Theorem 5, shown as solid lines.

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

8 L 1
-norm L 2 -norm L ∞ -norm Combined norm best (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 :
Figure 4: Bounds for the generalised Lyapunov exponent for the matrix product 9 with α = −3, β = 3.As before, estimates originating from L 1 -, L 2 -and L ∞ -norms are given in green, red and blue respectively.

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