Low-rank Parareal: a low-rank parallel-in-time integrator

In this work, the Parareal algorithm is applied to evolution problems that admit good low-rank approximations and for which the dynamical low-rank approximation (DLRA) can be used as time stepper. Many discrete integrators for DLRA have recently been proposed, based on splitting the projected vector field or by applying projected Runge–Kutta methods. The cost and accuracy of these methods are mostly governed by the rank chosen for the approximation. These properties are used in a new method, called low-rank Parareal, in order to obtain a time-parallel DLRA solver for evolution problems. The algorithm is analyzed on affine linear problems and the results are illustrated numerically.


Introduction
This work is concerned with the parallel-in-time integration of evolution problems for which the solution can be well approximated by a time-dependent low-rank matrix. In particular, we aim to solve approximately the evolution problem . X (t) = F(t, X (t)), t ∈ [0, T ], where X (t) is a matrix of size m × m. When the dimension m is large, the numerical solution of (1) can be very expensive since the matrix X (t) is usually dense. One way to alleviate this curse of dimensionality is to use low-rank approximations where, for every t, we approximate X (t) by Y (t) ∈ R m×m such that rank(Y (t)) = r m. The accuracy of this approximation will depend on the choice of the rank r . Here, X (t) is assumed to be square for notational convenience and all results can be easily formulated for rectangular X (t).
A popular paradigm to solve directly for the low-rank approximation Y (t) is the dynamical low-rank approximation (DLRA), first proposed in [25]. As defined later in Def. 5, DLRA leads to an evolution problem that is a projected version of (1). In the last decade, many discrete integrators for this projected problem have been proposed. One class of integrators consists in a clever splitting of the projector so that the resulting splitting method can be implemented efficiently. An influential example is the projector-splitting scheme proposed in [29]. Other methods that require integrating parts of the vector field can be found in [4,22]. Another approach, proposed in [9,24,35], is based on projecting standard Runge-Kutta methods (sometimes including their intermediate stages). Most of these methods are formulated for constant rank r . Rank adaptivity can be incorporated without much difficulty for splitting and for projected schemes; see [3,6,9,35]. Finally, given the importance of DLRA in problems from physics (like the Schrödinger and Vlasov equation), the integrators in [7,29] also preserve certain invariants, like energy. However, none of these time integrators consider a parallel-in-time scheme for DLRA, which is particularly interesting in the large-scale setting.
Parallel computing can be very effective and is even necessary to solve large-scale problems. While parallelization in space is well known, also the time direction can be parallelized to some extent when solving evolution problems. Over the last decades, various parallel-in-time algorithms have been proposed; see, e.g., the overviews [12,32]. Among these, the Parareal algorithm from [28] is one of the more popular algorithms for time parallelization. It is based on a Newton-like iteration, with inaccurate but cheap corrections performed sequentially, and accurate but expensive solves performed in parallel. This idea of solving in parallel an evolution problem as a nonlinear (discretized) system also appears in related methods like PFASST [8], MGRIT [11] and Space-Time Multi-Grid [17]. Theoretical results and numerical studies on a large numbers of cores show that these parallel-in-time methods can have good parallel performance for parabolic problems; see, e.g., [16,21,37]. So far, these methods did not incorporate a low-rank compression of the space dimension, which is the main topic of this work.

The Parareal algorithm
The Parareal iteration in Def. 1 below is given for constant time step h (hence T = N h) and for autonomous F. Both restrictions are not crucial but ease the presentation. The quantity X k n is an approximation for X (t n ) at time t n = nh and iteration k. The accuracy of this approximation is expected to improve with increasing k. Here, and throughout the paper, we denote dependency on the iteration index k as k , which should not be confused with the kth power.
Definition 1 (Parareal) The Parareal algorithm is defined by the following double iteration on k and n, (Initial approximation) Here, F h (X ) represents a fine (accurate) time stepper applied to the initial value X and propagated until time h. Similarly, G h (X ) represents a coarse (inaccurate) time stepper.
Given two time steppers, the Parareal algorithm is easy to implement. A remarkable property of Parareal is the convergence in a finite number of steps for k = n. It is well known that Parareal works well on parabolic problems but behaves worse on hyperbolic problems; see [18] for an analysis.

Dynamical low-rank approximation
Let M r denote the set of m × m matrices of rank r , which is a smooth embedded submanifold in R m×m . Instead of solving (1), the DLRA solves the following projected problem: Definition 2 (Dynamical low-rank approximation) For a rank r , the dynamical lowrank approximation of problem (1) where P Y is the l 2 -orthogonal projection onto the tangent space To analyze the approximation error of DLRA, we need the following standard assumptions from [23]. Here, and throughout the paper, · denotes the Frobenius norm.

Assumption 1 (DLRA assumptions)
The function F satisfies the following properties for all X , Y ∈ R m×m : -Lipschitz with constant L: In the analysis in Section 3, it is necessary to have < 0 for convergence. This holds when F is a discretization of certain parabolic PDEs, like the heat equation. In particular, for an affine function of the form [20,Ch. I.10]. The quantity ε r is called the modeling error and decreases when the rank r increases. For our problems of interest, this quantity is typically very small. Finally, the existence of L is only needed to guarantee the uniqueness of (1) but it will actually not appear in our analysis. We can therefore allow L to be large, as is the case for discretized parabolic PDEs.
Standard theory for perturbations of ODEs allows us to obtain the following error bound from the assumptions above: Theorem 1 (Error of DLRA [23]) Under Assumption 1, the DLRA verifies where φ h is the flow of the original problem (1) and ψ h r is the flow of its DLRA (5) for rank r .
The solution of DLRA (5) is quasi-optimal with the best rank approximation. This can be seen already in Theorem 1 for short time intervals. Similar estimates exist for parabolic problems [5] and for longer time when there is a sufficiently large gap in the singular values and when their derivatives are bounded [25].

Contributions
In this paper, we propose a new algorithm, called low-rank Parareal. As far as we know, this is the first parallel-in-time integrator for low-rank approximations. We analyze the proposed algorithm when the function F in (1) is affine. To this end, we extend the analysis of the classical Parareal algorithm in [14] to a more general setup where the coarse problem is different from the fine problem. We can prove that the method converges for big steps (large h) on diffusive problems ( < 0). In numerical experiments, we confirm this behavior. In addition, the method also performs well empirically with a less strict condition on h and on a non-affine problem.

Low-rank Parareal
We now present our low-rank Parareal algorithm for solving (1). Since the cost of most discrete integrators for DLRA scales quadratically 1 with the approximation rank, we take the coarse time stepper as DLRA with a small rank q. Likewise, the fine time stepper is DLRA with a large rank r . We can even take r = m, which corresponds to computing the exact solution as the fine time stepper since Y ∈ R m×m . Definition 3 (Low-rank Parareal) Consider two ranks q < r . The low-rank Parareal algorithm iterates (Initial approximation) where ψ h r (Z ) is the solution of (5) at time h with initial value Y 0 = Z , and T r is the orthogonal projection onto M r . The notations ψ h q and T q are similar but apply to rank q. The matrices E n are small perturbations such that rank(Y 0 n+1 ) = r + 2q and can be chosen randomly. 2 Observe that the rank of Y k n is at most r + 2q for all n, k. The low-rank structure is therefore preserved over the iterations. The matrices E n insure that each iteration has a rank between r and r + 2q. These matrices impact only the initial error but do not have any role in the convergence of the algorithm as is shown later in the analysis. An efficient implementation should store the low-rank matrices in factored form. In this context, the truncated SVD can be efficiently performed. The DLRA flows ψ h r and ψ h q can only be computed for relatively small problems. For larger problems, a suitable DLRA integrator must be used; see Sect. 4 for implementation details.

Convergence analysis
Let X n = X (t n ) be the solution of the full problem (1) at time t n . Let Y k n be the corresponding low-rank Parareal solution at iteration k. We are interested in bounding the error of the algorithm, for all relevant n and k. To this end, we make the following assumption: The function F is affine linear and autonomous, that is, with A : R m×m → R m×m a linear operator and B ∈ R m×m . 2 One could even take all initial Y 0 n+1 random, as is sometimes also done in standard Parareal. The important property is rank(Y 0 n+1 ) = r + 2q so that the low-rank Parareal iterations (9) are performed on the correct manifold.
The following lemma gives us a recursion for the Frobenius norm of the error. This recursion will be fundamental in deriving our convergence bounds later on when we generalize the proof for standard Parareal from [14].

Lemma 1 (Iteration of the error) Under the Assumptions 1 and 2, the error of low-rank Parareal verifies
The constants , ε q , ε r are defined in Assumption 1. Moreover, C r ,q and C q are the Lipschitz constants of T r − T q and T q .
Proof Our proof is similar to the one in [24] where first the continuous version of the approximation error of DLRA is studied. Denote by φ h (Z ) the solution of (1) at time h with initial value Y 0 = Z . By definition, the discrete error is We can interpret each term above as a flow from t n to t n+1 = t n + h. Denote these flows by Defining the continuous error as we then get the identity E k+1 n+1 = E(t n + h). We proceed by bounding E(t) . By definition of the flows above, we have (omitting the dependence on t in the notation) where the last equality holds since the function F is affine. Using Assumption 1 and Cauchy-Schwarz, we compute we therefore obtain the differential inequality Gr-nwall's lemma allows us to conclude From (12), we get Denoting T ⊥ r = I − T r and T r ,q = T r − T q , we get after rearranging terms Taking norms gives where C r ,q and C q are the Lipschitz constants of T r ,q and T q respectively. Combining inequalities (13) and (14) gives the statement of the lemma.
We now study the error recursion (11) in more detail. To this end, let us slightly rewrite it as with the non-negative constants Our first result is a linear convergence bound, up to the DLRA approximation error. It is similar to the linear bound for standard Parareal.
Proof Define e k = max n≥0 E k n . Taking the maximum for n ≥ 0 of both sides of (15), we obtain By assumption, 0 ≤ β < 1 and we can therefore obtain the recursion with solution By assumption, we also have 0 ≤ α 1−β < 1, which allows us to obtain the statement of the theorem.
Next, we present a more refined superlinear bound. To this end, we require the following technical lemma that solves the equality version of the double iteration (15). A similar result, but without the κ term and only as an upper bound, already appeared in [14, Thm. 1]. Our proof is therefore similar but more elaborate.
Lemma 2 Let α, β, γ , κ ∈ R be any non-negative constants such that α < 1 and β < 1. Let e k n be a sequence depending on n, k ∈ N such that Then, Proof The idea is to use the generating function ρ k (ξ ) = ∞ n=1 e k n ξ n for k ≥ 1. Multiplying (18) by ξ n+1 and summing over n, we obtain Since e k 0 = 0 for all k, this gives the relations We can therefore obtain the linear recurrence .

Its solution satisfies
It remains to compute the coefficients in the power series of the above formula since by definition of ρ k (ξ ) = ∞ n=1 e k n ξ n they equal the unknowns e k n . The binomial series formula for |z| < 1, together with the Cauchy product gives Hence, the first term in ρ k (ξ ) satisfies while the second term can be written as Putting everything together, we have Finally, we can identify the coefficient e k n in front of ξ n with those on the right-hand side. The coefficient for ξ n in the first term is clearly nonzero only when n ≥ k + 1. In the second term, there is only one m for every j such that m + j = n. Substituting m = n − j allows us to identify the coefficient of ξ n .
Using the previous lemma, we can obtain a convergence bound that is superlinear in k.
Proof Define e k n = E k n . By Lemma 1, the terms e k n verify the relation described in Lemma 2 with = replaced by ≤ in (18). Hence, the solution (19) from Lemma 2 will be an upper bound for e k n . Since 0 ≤ α + β < 1 and using the binomial series formula (20), we bound the first term in (19) as For 0 ≤ i ≤ n − k − 1 and n ≥ k + 1, observe that Since 0 ≤ β < 1, we can therefore bound the second term as The conclusion now follows by the definition of γ .
The proof above can be modified to obtain a simple linear bound that is similar but different to the one from Theorem 2: Theorem 4 (Another linear convergence bound) Under Assumptions 1 and 2, and if α + β < 1, the error of low-rank Parareal satisfies for all n, k ∈ N the bound where α, β, κ are defined in (16).
Proof We repeat the proof for the superlinear bound but this time, the second term is bounded as (20) is

Remark 1 In the proof above, yet another bound based on
This time we recover the linear bound from Theorem 2.

Summary of the convergence bounds
In the previous section, we have proven four upper bounds for the error of low-rank Parareal. The first is directly obtained from Lemma 2. It is the tightest bound but its expression is too unwieldy for practical use. The other three bounds can be summarized as with B n,k Rate of (23) in k

Superlinear
Each of these practical bounds describes different phases of the convergence, and none is always better than the others. In Fig. 1, we have plotted all four bounds for realistic values of α and β. We took κ = 10 −15 ≈ ε mach since it only determines the stagnation of the error and would interfere with judging the transient behavior of the convergence plot. Furthermore, the errors e 0 n = γ = 1 at the start of the iteration k = 0 were chosen arbitrarily since they have little influence on the results.
The bounds above depend on α = e h C r ,q and β = e h C q , where C q and C r ,q are the Lipschitz constants of T q and T r ,q respectively; see (16). While it seems difficult to give a priori results on the size C q and C r ,q , we can bound them up to first order in the theorem below. Note also that in the important case of < 0, the constants α and β can be made as small as desired by taking h sufficiently large.

Theorem 5 (Lipschitz constants) Let A,Ã ∈ R m×n . Then
where σ q is the qth singular value of A. Moreover, Proof For the first inequality, we refer to [2,Theorem2] and [10,Theorem24].The second inequality follows from the first by the triangle inequality, In many applications with low-rank matrices, the singular values of the underlying matrix are rapidly decaying. In particular, when the singular values decay exponentially like σ k ≈ e −ck for some c > 0, we have This last quantity decreases quickly to 1 when c grows. Even for c = 1, it is less than 1.6. We therefore see that the constants in Theorem 5 are not too large in this case.

Remark 2
In the analysis, a sufficiently large gap in the singular values is required at both the coarse rank and the fine rank. In our experiments, we observed that such a gap is indeed required at the coarse rank, but not at the fine rank. It suggests that the bound (25) can therefore probably be improved.

Numerical experiments
We now show numerical experiments for our low-rank Parareal algorithm. We implemented the algorithm in Python 3.10 and all computations were performed on a MacBook Pro with a M1 processor and 16GB of RAM. The complete code is available at GitHub so that all the experiments can be reproduced. The DLRA steps are solved by the second-order projector-splitting integrator from [29]. Since the problems considered are stiff, we used sufficiently many substeps of this integrator so that the coarse and fine solvers within low-rank Parareal can be considered exact.

Lyapunov equation
Consider the differential Lyapunov equation, where A ∈ R m×m is a symmetric matrix, and C ∈ R m×k is a tall matrix for some k ≤ m. This initial value problem admits a unique solution for t ∈ [0, T ] for any T > 0. The most typical example of (27) is the heat equation on a square with separable source term. Other applications can be found in [31].

Assumption 3 The matrix A ∈ R m×m is symmetric and strictly negative definite.
Under Assumption 3, the one-sided Lipschitz constant for (27) (27) as which can be easily verified by differentiation using properties of the matrix exponential e tA (Z ) = e t A Ze t A .
The following result shows that the solution of (27) can be well approximated by low rank. It is the analogue to a similar result for the algebraic Lyapunov equation A (X ) = CC T . The latter result is well known, but we did not find a proof for the former in the literature.
Proof The aim is to approximate the following two terms that make up the closed-form solution X (t) in (28): The first term X 1 (t) can be treated directly. By the truncated SVD, the initial value satisfies By Assumption 3, the operator A is full rank. We therefore obtain where rank(Y 1 (t)) = rank(Y 0 ) = r 0 and E 1 (t) 2 ≤ e t σ r 0 +1 (X 0 ) since = 2 max i λ i (A). Next, we focus on the second term X 2 (t). Like above, the source term can be decomposed as By linearity of the Lyapunov operator, we therefore obtain Denote M = e tA D − D. By definition of the Lyapunov operator A , we have As studied in [34] and then improved in [1], the singular values of the solution S are bounded as σ rank(M)ρ+1 (S) where κ A = A 2 A −1 2 and 0 ≤ ρ ≤ m. Since rank(M) ≤ 2 rank(D) = 2r by assumption on D, the bound (31) then implies that where rank(Y 2 (t)) ≤ 2r ρ and where the last inequality holds by properties of the logarithmic norm μ of A which equals ; see [36,Proposition 2.2]. Moreover, we can bound the last term in (30) as Putting (29) and (30) together, we obtained which proves the statement of the lemma.
The lemma shows that if X 0 and CC T have good low-rank approximations, then the solution X (t) of the differential Lyapunov equation has comparable low-rank approximations as well on [0, T ]. Since < 0, we can even take T → ∞ and recover essentially the low-rank approximability of the Lyapunov equation X (∞) = A −1 (CC T ). This is clearly visible when X 0 and CC T are exactly of low rank, which we state as a simple corollary for convenience.

Corollary 1
Under Assumption 3 and assuming that rank(X 0 ) = r 0 , rank(CC T ) = r , the solution X (t) of (28) has an approximation Y (t) of rank at most r 0 + 2r ρ for any 0 ≤ ρ ≤ m with error The corollary clearly shows that the approximation error decreases exponentially when the approximation rank increases linearly via ρ. Furthermore, we see that the condition number of the matrix A has only a mild influence due to log(κ A ). [26]. In that work, the authors solve (27) with exact low-rank X 0 = Z Z T and CC T using a Krylov subspace method. More specifically, with U k an orthonormal matrix that spans the block Krylov space K k (A, [C Z]), the projected Lyapunov equation

Remark 3 Corollary 1 can be compared to a similar result in
The approximation error of X k (t) is studied in [26,Theorem4.2]. Since rank(X k (t)) ≤ k(rank(Z ) + rank(C)), we therefore also get a result on the low-rank approximability of (27). This bound is, however, worse than ours since it does not give zero error for t = 0 and k = 1, for example. On the other hand, it is a bound for a discrete method whereas our Lemma 3 and Corollary 1 are statements about the exact solution.
We now apply the low-rank Parareal algorithm to the differential Lyapunov equation (27). Let A = Δ dx be the n × n discrete Laplacian with zero Dirichlet boundary conditions obtained by standard centered differences on [−1, 1]. The Lyapunov equation is therefore a model for the 2D heat equation on Ω = [−1, 1] 2 . In the experiments, we used n = 100 spatial points and the time interval [0, T ] = [0, 2]. The matrix C for the source is generated randomly with singular values σ i = 10 −5(i−1) where i = 1, 2, . . . so that its numerical rank is 4. In order to have a realistic initial value, X 0 is obtained as the exact solution at time t = 0.01 of the same ODE but with a random initial valueX 0 with singular values σ i = 10 −(i−1) . Figure 2 is a 3D plot of the solution over time on Ω with its corresponding singular values. As we can see, the solution becomes almost stationary at t = 1.0. In addition, it stays low-rank over time in agreement to Lemma 3. Moreover, the singular values suggest to take the fine rank r = 16 for an error of the fine solver of order 10 −12 .
The convergence of the error of the low-rank Parareal algorithm is shown in Fig. 3. The algorithm converges linearly from the coarse rank solution to the fine rank solution. Figure 3a suggests that the coarse rank does not influence the convergence rate and it only reduces the initial error. This is consistent with our analysis. Indeed, since the singular values are exponentially decaying, the singular gap is approximately constant; see (26). Hence, the constants α and β from (16) that determine the convergence rate do not depend on the coarse rank q; as is shown up to first order in Theorem 5. Figure 3b shows that, similarly, the convergence rate does not depend on the fine rank either, although it limits the final error.
In Fig. 4a, we investigate the convergence for several sizes n. Even though the problem is stiff, the convergence does not seem influenced by the size of the problem. Figure 4b shows the error of the algorithm applied to the problem with several step sizes. According to our analysis, the convergence is faster when the stepsize h is large; see (16).

Parametric cookie problem
We now solve a simplified version of the parametric cookie problem from [27]. Consider the ODE where the sparse matrices A 0 , A 1 ∈ R 1580×1580 , b ∈ R 1580 , and C 1 = diag(c 1 1 , c 2 1 , . . . , c p 1 ) are given in [27]. The aim of this problem is to solve a heat problem simultaneously with several heat coefficients, denoted by The singular values of the reference solution are shown in Fig. 5. The stationary solution has good low-rank approximations, as was proved in [27,Thm. 2.4]. The singular value decay suggests that a fine rank r = 16 leads to full numerical accuracy.
In Fig. 6, we applied the low-rank Parareal algorithm with several coarse ranks q and fine ranks r . Like for the Lyapunov equation, it seems that the convergence rate does not depend on the coarse rank q. In agreement to our analysis (see Fig. 1), the convergence is linear in the first iterations and superlinear in the last iterations. In addition, the convergence is not influenced by the fine rank r .    . 6 Convergence of the errorof low-rank Parareal for the parametric cookie problem (32). Influence of the coarse and fine ranks

Riccati equation
The Riccati differential equation is given by .
where X ∈ R m×m , A ∈ R m×m , C ∈ R k×m , and S ∈ R m×m . We note that this is no longer an ODE with an affine vector field and hence our theoretical results do not apply here. As already studied in [33], we take S = I and A is the spatial discretization of the diffusion operator on the spatial domain Ω = [0, 1]. Furthermore, we take α(x) = 2 + 2 cos(2π x) and λ = 1. The discretization is done by the finite volume method, as described in [15]. The tall matrix C ∈ R k×m is obtained from k independent vectors As for the other problems, the singular values of the solution (shown in Fig. 7) indicate that we can expect good low-rank approximations on [0, T ]. We choose the fine rank r = 18. The convergence of low-rank Parareal is shown in Fig. 8. Unlike the previous problems, the coarse rank q has a more pronounced influence on the behavior of the convergence. While our theoretical results do not hold for this nonlinear problem, we still see that low-rank Parareal converges linearly when the coarse rank q is sufficiently large (q = 6, q = 8). The convergence is slower (but still superlinear) when q = 4. This could be due to the non-constant gaps in the singular values. The influence of the fine rank r is more like for the linear problems.

Rank-adaptive algorithm
Since the approximation rank of the solution is usually not known a priori, it is more convenient for the user to supply an approximation tolerance than an approximation Fig. 8 Convergence of the error of low-rank Parareal for the Riccati problem (33). Influence of the coarse and fine ranks rank. Even though the rank can change to satisfy the tolerance during the truncation steps, Algorithm 3 can be easily reformulated for such a rank adaptive setting. The key idea is to fix the coarse rank to keep the cost of the coarse solver low, while the fine rank is determined by a fine tolerance.
Definition 4 (Adaptive low-rank Parareal) Consider a small fixed rank q and a fine tolerance τ . The adaptive low-rank Parareal algorithm iterates where the notation is similar to that of the previous Def. 3, except for T τ which represents the rank-adaptive truncation. In particular, T τ (Y ) is the best rank q approximation of Y so that the (q + 1)st singular value of Y equals the tolerance τ . The matrices E n are small perturbations, randomly generated such that rank(Y 0 n+1 ) = rank(Y 0 ) and its smallest singular value is larger than the fine tolerance τ . Figure 9 shows the numerical behavior of this rank-adaptive algorithm. As we can see, the algorithm behaves as desired. Figure 9a shows the algorithm applied with several tolerances and is comparable to Fig. 3b with several fine ranks. Figure 9b shows the rank of the solution over time. Already after two iterations, the rank is reduced to almost the numerical rank of the exact solution and the rank does not change much for the rest of the iterations.

Conclusion and future work
We proposed the first parallel-in-time algorithm for integrating a dynamical low-rank approximation (DLRA) of a matrix evolution equation. The algorithm follows the Fig. 9 Adaptive low-rank Parareal. On the left, the algorithm is applied with several tolerances. On the right, the rank of the solution over time is shown for several iterations traditional Parareal scheme but it uses DLRA with a low rank as coarse integrator, whereas the fine integrator is DLRA with a higher rank. Taking into account the modeling error of DLRA, we presented an analysis of the algorithm and showed linear convergence as well as superlinear convergence under common assumptions and for affine linear vector fields, up to the modeling error.
In our numerical experiments, the algorithm behaved well on diffusive problems, which is similar to the original Parareal algorithm. Due to the significant difference in computational cost for the fine and coarse integrators, it is reasonable to expect good speed-up in actual parallel implementations. A proper parallel implementation to verify this claim is a natural future work. It may however be more appropriate to first generalize more efficient parallel-in-time algorithms, like Schwarz waveform relaxation and multigrid methods [13], to DLRA.
Since DLRA can also be used to obtain low-rank tensor approximations [30], another future work is to extend low-rank Parareal to tensor DLRA. Finally, our theoretical analysis assumes that the ODE has an affine vector field. Since this assumption was only needed in one step of the proof of Lemma 1, it might be possible that it can be relaxed to include certain non-linear vector fields.