Fast iteration of cocyles over rotations and Computation of hyperbolic bundles

In this paper, we develop numerical algorithms that use small requirements of storage and operations for the computation of hyperbolic cocycles over a rotation. We present fast algorithms for the iteration of the quasi-periodic cocycles and the computation of the invariant bundles, which is a preliminary step for the computation of invariant whiskered tori.

The calculation of invariant bundles for cocycles is a preliminary step for the calculation of whiskered invariant tori. Indeed, these algorithms require the computation of the projections over the linear subspaces of the linear cocycle.
The main example of a cocycle in this paper is for K a parameterization of an invariant torus satisfying the invariance equation Other examples appear in discrete Schrödinger equations [Pui02].  Note that (2.3) and (2.4) are the exact analogues of (2.1) and (2.2) in a continuous context. Moreover, if M(T ℓ ) ⊂ G, where G is a sub-algebra of the Lie algebra of the Lie group G, then M(R, T ℓ ) ⊂ G.
The main example for us of a continuous in time cocycle will be where K is a solution of the invariance equation and X is a Hamiltonian vector field. In this case, the cocycle M(θ, t) is symplectic.
Definition 1. Given 0 < λ < µ we say that a cocycle M(n, θ) (resp. M(t, θ)) has a λ, µ− dichotomy if for every θ ∈ T ℓ there exist a constant c > 0 and a splitting depending on θ, or, in the continuous time case The notation E s and E u is meant to suggest that an important case is the splitting between stable and unstable bundles. This is the case when λ < 1 < µ and the cocycle is said to be hyperbolic. Nevertheless, the theory developed in this section assumes only the existence of a spectral gap.
In the application to the computation of tori, M(θ) = (DF • K)(θ) and x θ = K(θ). The existence of the spectral gap means that at every point of the invariant torus K(θ) one has a splitting so that the vectors grow with appropriate rates λ, µ under iteration of the cocycle. In the case of the invariant torus, it can be seen that the cocycle is just the fundamental matrix of the variational equations so that the cocycle describes the growth of infinitesimal perturbations.
It is well known that the mappings θ → E s,u x θ are C r if M(., ) ∈ C r for r ∈ N ∪ {∞, ω} [HL06c]. This result uses heavily that the cocycles are over a rotation.
A system can have several dichotomies, but for the purposes of this paper, the definition 1 will be enough, since we can perform the analysis presented here for each gap.
One fundamental problem for subsequent applications is the computation of the invariant splittings (and, of course, to ensure their existence). The computation of the invariant bundles is closely related to the computation of iterations of the cocycle.
The first algorithm that comes to mind, is an analogue of the power method to compute leading eigenvalues of a matrix. Given a typical vector (x θ , v) ∈ E u , we expect that, for n ≫ 1, M(n, θ)v will be a vector in E u . Even if there are issues related to the θ dependence, this may be practical if E u is a 1-dimensional bundle.

Equivalence of cocycles, reducibility
Reducibility is a very important issue in the theory of cocycles. We have the following definition.
Definition 2. We say that a cocycle M (θ) is equivalent to another cocycle M(θ) if there exists a matrix valued function Q : (4.1) It is easy to check that M being equivalent to M is an equivalence relation. If M is equivalent to a constant cocycle (i.e. independent of θ), we say that M is "reducible." The important point is that, when (4.1) holds, we have In particular, if M is a constant matrix, we have so that the iterations of reducible cocycles are very easy to compute. We will also see that one can alter the numerical stability properties of the iterations of cocycles by choosing appropriately Q. In that respect, it is also important to mention the concept of "quasi-reducibility" introduced by Eliasson [Eli01].

Algorithms for fast iteration of cocycles over rotations
In its simplest form, the algorithm for fast iteration of cocyles is: Algorithm 1 (Iteration of cocycles 1). Given M(θ), compute → ω and iterate the procedure.
We refer to M as the renormalized cocycle and the procedure as a renormalization procedure.
The important property is that applying k times the renormalization procedure described in Algorithm 1 amounts to compute the cocycle M(2 k , θ).
Then, if we discretize the matrix M(θ) taking N points (or N Fourier modes) and denote by C(N) the number of operations required to perform a step of Algorithm 1, we can compute 2 k iterates at a cost of kC(N) operations.
Notice that the important point is that the cost of computing 2 k iterations is proportional to k. Of course, the proportionality constant depends on N. The form of this dependence on N depends on the details on how we compute the shift and the product of matrix valued functions.
There are several alternatives to perform the transformation (5.1). The main difficulty arises from the fact that, if we have points on a equally spaced grid, then θ + ω will not be in the same grid. We have at least three alternatives: (1) Keep the discretization by points in a grid and compute M(θ+ω) by interpolating with nearby points. (2) Keep the discretization by points in a grid but note that the shift is diagonal in Fourier space. Of course, the multiplication of the matrix is diagonal in real space.
(3) Keep the discretization in Fourier space but use the Cauchy formula for the product.
The cost factor of each of these alternatives is, respectively, Besides their cost, the above algorithms may have different stability and roundoff properties. We are not aware of any study of these stability or round-off properties. The properties of interpolation depend on the dimension. In each of the cases, the main idea of the method is to precompute some blocks of the iteration, store them and use them in future iteration. One can clearly choose different strategies to group the blocks. Possibly, different methods can lead to different numerical (in)stabilities. At this moment, we lack a definitive theory of stability which would allow us to choose the blocks.
Next, we will present an alternative consisting of using the QR decomposition for the iterates. As described, for instance in [Ose68,ER85,DVV02], the QR algorithm seems to be rather stable to compute iterates. One advantage is that, in the case of several gaps, it can compute all the eigenvalues in a stable way.
Algorithm 2 (Computation of cocycles with QR). Given M(θ) and a QR decomposition of M(θ), perform the following operations: and iterate the procedure.
Since the QR decomposition is a fast algorithm, the total cost of the implementation depends on the issues previously discussed (see costs in (5.2)). Instead of using QR decomposition, one can also perform a SV D decomposition (which is somewhat slower).
In the case of one-dimensional maps, one can be more precise in the description of the method. Indeed, if the frequency ω has a continued fraction expansion ω = [a 1 , a 2 , . . . , a n , . . .], it is well known that the denominators q n of the convergents of ω (i.e. p n /q n = [a 1 , . . . , a n ]) satisfy q n = a n q n−1 + q n−2 , q 1 = a 1 , q 0 = 1.
As a consequence, we can consider the following algorithm for this particular case: Algorithm 3 (Iteration of cocycles 1D). Given ω = [a 1 , . . . , a n , . . .] and the cocycle over T ω generated by M(θ), define Then, for n ≥ 2 is a cocycle over ω n−1 = a n ω n−1 + ω n−2 and we have The advantage of this method is that the effective rotation is decreasing to zero so that the cocycle is becoming in some ways similar to the iteration of a constant matrix. This method is somehow reminiscent of some algorithms that have appeared in the mathematical literature [Ryc92,Kri99b,Kri99a].
6. The "straddle the saddle" phenomenon and preconditioning The iteration of cocycles has several pitfalls compared with the iteration of matrices.
The main complication from the numerical point of view is that the (un)stable bundle does depend on the base point.
In this section we describe a geometric phenomenon that causes some instability in the iteration of cocycles. This instability -which is genuine-becomes catastrophic when we apply some of the fast iteration methods described in Section 5. The phenomenon we will discuss was already observed in [HL06a].
Since we have the inductive relation, we see that we can think of computing M(n, θ) by applying M(n − 1, θ + ω) to the column vectors of M(θ).
The j th -column of M, which we will denote by m j (θ), can be thought geometrically as an embedding from T ℓ to R 2d and is given by M(θ)e j where e j is the j th vector of the canonical basis of R 2d . If the stable space of M(n − 1, θ + ω) has codimension ℓ or less, there could be points θ * ∈ T ℓ such that m j (θ * ) ∈ E s x θ * and such that for every θ = θ * we have m j (θ) / ∈ E s x θ . Clearly, the quantity M(n − 1, θ * + ω)m j (θ * ) decreases exponentially. Nevertheless, for all θ in a neighborhood of θ * such that θ = θ * M(n − 1, θ + ω)m j (θ) will grow exponentially. The direction along which the growth takes place depends on the projection of m j (θ) on E u x θ+ω . For example, in the case d = 2, ℓ = 1 and the stable and unstable directions are one dimensional, the unstable components will have different signs and the vectors M(n − 1, θ + ω)m j (θ) will align in opposite directions. An illustration of this phenomenon happens in Figure 1. The transversal intersection of the range of m j (θ) with E s is indeed a true phenomenon, and it is a true instability.
Unfortunately, if m j (θ) is very discontinuous as a function of θ, the discretization in Fourier series or the interpolation by splines will be extremely inaccurate so that the Algorithm 1 fails to be relevant.
This phenomenon is easy to detect when it happens because the derivatives grow exponentially fast in some localized spots.
One important case where the straddle the saddle is unavoidable is when the invariant bundles are non-orientable. This happens near resonances (see [HL07]). In [HL07], it is shown that, by doubling the angle the case of resonances can be studied comfortably because then, non-orientability is the only obstruction to the triviality of the bundle.
6.1. Eliminating the "straddle the saddle" in the one-dimensional case. Fortunately, once the phenomenon is detected, it can be eliminated. The main idea is that one can find an equivalent cocycle which does not have the problem (or presents it in a smaller extent).
In more geometric terms we observe that, even if the stable and unstable bundles are geometrically natural objects, the decomposition of a matrix into columns is coordinate dependent. Hence, if we choose some coordinate system which is reasonably close to the stable and unstable manifolds and we denote by Q the change of coordinates, then the cocycle is close to a constant. Remark that this is true only in the one-dimensional case. The picture is by far more involved when the bundles have higher rank.
This may seem somewhat circular, but the circularity can be broken using continuation methods. Given a cocycle which is close to constant, fast iteration methods work and they allow us to compute the splitting. Then if we have computed Q for some M, we can use it to precondition the computation of neighboring M.
The algorithms for the computation of bundles will be discussed next.
6.2. Computation of rank-1 stable and unstable bundles using iteration of cocycles. The algorithms described in the previous section provide a fast way to iterate the cocycle. We will see that this iteration method, which is similar to the power method, gives the dominant eigenvalue λ max (θ) and the corresponding eigenvector m(θ).
The methods based on iteration rely strongly on the fact that the cocycle has one dominating eigenvalue which is simple.
Consider that we have performed k iterations of the cocycle (of course we perform scalings at each step) and we have computed M(n, θ), with n = 2 k . Then, one can easily read the dominant rank-1 bundle from the QR decomposition of the cocycle M(n, θ), just taking the column of Q associated to the largest value in the diagonal of the upper triangular matrix R. One obtains a vector m(θ + 2 k ω) (and therefore m(θ) by performing a shift of angle −2 k ω) of modulus 1 spanning the unstable manifold. Since, we have then As it is standard in the power method, we perform scalings at each step dividing all the entries in the matrix M(θ) by the maximum value among the entries of the matrix.
Hence, for the simplest case that there is one dominant eigenvalue, the method produces a section m (spanning the unstable subbundle) and a real function λ max , which represents the dynamics on the rank 1 unstable subbundle, such that M(θ)m(θ) = λ max (θ)m(θ + ω).
Following [HL06b], under certain non-resonant conditions which are satisfied in the case of the stable and unstable subspaces, one can reduce the 1-dimensional cocycle associated to M and ω to a constant. Hence, we look for a positive function p(θ) and a constant µ such that λ max (θ)p(θ) = µp(θ + ω). (6.1) If λ max (θ) > 0, we take logarithms on both sides of the equation (6.1). This leads to log λ max (θ) + log p(θ) = log µ + log p(θ + ω), and taking log µ to be the average of log λ max (θ) the problem reduces to solve for log p(θ). The case λ max (θ) < 0 is analogous. Of course, p(θ) and µ can be obtained just exponentiating.
Stability and Instability in Mechanical Systems (SIMS08), for whose hospitality we are very grateful. YS would like to thank the hospitality of the department of Mathematics of University of Texas at Austin, where part of this work were carried out.