Subspace Acceleration for a Sequence of Linear Systems and Application to Plasma Simulation

We present an acceleration method for sequences of large-scale linear systems, such as the ones arising from the numerical solution of time-dependent partial differential equations coupled with algebraic constraints. We discuss different approaches to leverage the subspace containing the history of solutions computed at previous time steps in order to generate a good initial guess for the iterative solver. In particular, we propose a novel combination of reduced-order projection with randomized linear algebra techniques, which drastically reduces the number of iterations needed for convergence. We analyze the accuracy of the initial guess produced by the reduced-order projection when the coefficients of the linear system depend analytically on time. Extending extrapolation results by Demanet and Townsend to a vector-valued setting, we show that the accuracy improves rapidly as the size of the history increases, a theoretical result confirmed by our numerical observations. In particular, we apply the developed method to the simulation of plasma turbulence in the boundary of a fusion device, showing that the time needed for solving the linear systems is significantly reduced.


Introduction
The numerical solution of time-dependent partial differential equations (PDEs) often leads to sequences of linear systems of the form where t 0 < t 1 < t 2 < • • • is a discretization of time t, and both the system matrix A(t i ) ∈ R n×n and the right-hand side b(t i ) ∈ R n depend on time.
Typically, the systems (1) are available only consecutively.Such sequences of linear systems can arise in a number of applications, including implicit time stepping schemes for the solution of PDEs or iterative solutions of non-linear equations and optimization problems.A relevant example is given by timedependent PDEs solved in presence of algebraic constraints.In this case, even when an explicit time stepping method is used to evolve the nonlinear PDE, the discretization of the algebraic constraints leads to linear systems that need to be solved at every (sub-)timestep.This is the case of the simulation of turbulent plasma dynamics [10], where a linear constraint (Maxwell equations) is imposed upon the plasma dynamics described by a set of non linear fluid or kinetics equations.The linear systems resulting from the discretized algebraic constraints may feature millions of degrees of freedom, hence their solution is often computationally very expensive.One usually expects that the linear system (1) changes slowly in subsequent time steps.This work is focused on exploiting this property to accelerate iterative solvers, such as CG [17] for symmetric positive definite matrices and GMRES [24] for general matrices.An obvious way to do so is to supply the iterative solver for the timestep t i+1 with the solution of (1) at timestep t i , as initial guess.As a more advanced technique, in the context of Krylov subspace methods, subspace recycling methods [26] such as GCROT [7] and GMRES-DR [21] have been proposed.Such methods have been developed in the case of a single linear system, to enrich the information when restarting the iterative solver.The idea behind is often to accelerate the convergence by suppressing parts of the spectrum of the matrix, including the corresponding approximate invariant subspace in the Krylov minimization subspace.GCROT and GMRES-DR have then been adapted to sequences of linear systems in [22], recycling selected subspaces from one system to the next.For this class of methods to be efficient, it is necessary that the sequence of matrices undergoes local changes only, that is, the difference A(t i+1 ) − A(t i ) is computationally cheap to apply.For example, one can expect this difference matrix to be sparse when time dependence is restricted to a small part of the computational domain, e.g., through timedependent boundary conditions.We refer to [26] for a more complete survey of subspace recycling methods and their applications.In [5], subspace recycling was combined with goal-oriented POD (Proper Orthogonal Decomposition) in order to limit the size of the subspaces involved in an augmented CG approach.Simplifications occur when the matrices A(t i ) are actually a fixed matrix A shifted by different scalar multiples of the identity matrix, because Krylov subspaces are invariant under such shifts.In the context of subspace recycling, this property has been exploited in, e.g., [27], and in [25] it is shown how a smoothly varying right-hand side can be incorporated.
When A(t i ) and b(t i ) in ( 1) are samples of smooth matrix/vector-valued functions, one expects that the subspace of the previously computed solutions contains a very good approximation of the current one.This can be exploited to construct a better initial guess, either explicitly through (polynomial) extrapolation, or implicitly through projection techniques.Examples of the extrapolation approach include polynomial POD extrapolation [14], weighted group extrapolation methods [30] and a stabilized, least-squares polynomial extrapolation method [1], for the case that only the right-hand side evolves in time.For the same setting, projection techniques have been introduced by Fischer [11].Following this first work, several approaches have been developed to extract an initial guess from the solution of a reduced-order model, constructed from projecting the problem to a low-dimensional subspace spanned by previous solutions.In [28], such an approach is applied to fully implicit discretizations of nonlinear evolution problems, while [20] applies the same idea to the so called IMPES scheme used for simulating two-phase flows through heterogeneous porous media.
In this paper, we develop a new projection technique for solving sequences of linear systems that combines projection with randomized linear algebra techniques, leading to considerably reduced cost.Moreover, a novel convergence analysis of the algorithm is carried out to show its efficiency.This is also proved numerically by applying the algorithm to the numerical simulation of turbulent plasma in the boundary of a fusion device.
The rest of this paper is organized as follows.In Section 2, we first discuss general subspace acceleration techniques based on solving a projected linear system and then explain how randomized techniques can be used to speed up existing approaches.In Section 3, a convergence analysis of these subspace acceleration techniques is presented.In Section 4 we first discuss numerical results for a test case to demonstrate the improvements that can be attained by the new algorithm in a somewhat idealistic setting.In Section 5 our algorithm is applied to large-scale turbulent simulation of plasma in a tokamak, showing a significant reduction of computational time.

Algorithm
The algorithm proposed in this work for accelerating the solution of the sequence of linear systems (1) uses randomized techniques to lower the cost of a POD-based strategy, such as the one proposed in [20].Recall that we aim at solving the linear systems A(t i )x(t i ) = b(t i ) consecutively for i = 0, 1, • • • .We make no assumption on the symmetry of A(t i ) ∈ R n×n and thus GMRES is an appropriate choice for solving each linear system.Supposing that, at the ith timestep, M previous solutions are available, we arrange them into the history matrix where the notation on the right-hand side indicates the concatenation of columns.Instead of using the complete history, which may contain redundant information, one usually selects a subspace S ⊂ span(X) of lower dimension m ≤ M .Then, the initial guess for the ith linear system is obtained from choosing the element of S that minimizes the residual: where the columns of Q ∈ R n×m contain an orthonormal basis of S. We use • 2 to denote the Euclidean norm for vectors and the spectral norm for matrices.The described approach is summarized in Algorithm 1, which is a template that needs to be completed by an appropriate choice of the subspace S, in Sections 2.1 and 2.2.
Algorithm 1 Solution of ith linear system If the complete history is used, S = span(X), then computing Q via a QR decomposition [13], as required in Step 2, costs O(M 2 n) operations.In addition, setting up the linear least-squares problem in Step 3 of Algorithm 1 requires M (sparse) matrix-vector products in order to compute A(t i )Q.The standard approach for solving the linear least-squares problem proceeds through the QR decomposition of that matrix and costs another O(M 3 + M 2 n) operations.This strong dependence of the cost on M effectively forces a rather small choice of M , neglecting relevant components of the solutions that could be contained in older solutions only.In the following, we discuss two strategies to overcome this problem.

Proper Orthogonal Decomposition
An existing strategy [20] to arrive at a low-dimensional subspace S ⊂ span(X) uses a POD approach [19] and computes the orthonormal basis Q for S through a truncated SVD (Singular Value Decomposition) of X; see Algorithm 2. Note that only the first m left singular vectors Thanks to basic properties of the SVD, the basis Q POD enjoys the following optimality property [29]: where • F denotes the Frobenius norm and In words, the choice Q POD minimizes the error of orthogonally projecting the columns of X onto an m-dimensional subspace.
The relation to the singular values of X established in (2) also allows one to choose m adaptively, e.g., by choosing m such that most of the variability in the history matrix X is captured.At every time step, the history matrix X gets modified by removing its first column and appending a new last column.The most straightforward implementation of Algorithm 2 would compute the SVD needed of Step 2 from scratch at every time step, leading to a complexity of O(nM 2 ) operations.In principle, SVD updating techniques, such as the ones presented in [4] and [6], could be used to reduce this complexity to O(mn + m 3 ) for every time step.However, in the context of our application, there is no need to update a complete SVD (in particular, the right singular vectors are not needed) and the randomized techniques discussed in the next section seem to be preferable.

Randomized range finder
In this section, an alternative to the POD method (Algorithm 1) for generating the low-dimensional subspace S ⊂ span(X) is presented, relying on randomized techniques.The randomized SVD from [16] applied to the n × M history matrix X proceeds as follows.First, we draw an M ×m Gaussian random matrix Z, that is, the entries of Z are independent and identically distributed (i.i.d) standard normal variables.Then the so-called sketch is computed, followed by a reduced QR decomposition Ω = QR.This only involves the n × m matrix Ω, which for m ≪ M is a significant advantage compared to Algorithm 2, which requires the SVD of an n × M matrix.The described procedure is contained in lines 2-4 and 11 of Algorithm 3 below.
According to [16,Theorem 10.5], the expected value (with respect to Z) of the error returned by the randomized SVD satisfies where we partition m = r + p for a small oversampling parameter p ≥ 2. Also, the tail bound from [16,Theorem 10.7] implies that it is highly unlikely that the error is much larger than the upper bound (4).Comparing (4) with the error (2), we see that the randomized method is only a factor √ 2 worse than the optimal basis of roughly half the size produced by POD.As we also see in our experiments of Section 4, this bound is quite pessimistic and usually the randomized SVD performs nearly as good as POD using bases of the same size.
Algorithm 3 Method 2 (Randomized Range Finder) to generate basis Q , matrices Ω and Z from previous time step (see ( 5)) Draw new Gaussian random vector z M ∈ R m 9: Instead of performing the randomized SVD from scratch in every timestep, one can easily exploit the fact that only a small part of the history matrix is modified.To see this, let us consider the sketch from the previous timestep: Comparing with (3), we see that the sketch Ω of the current timestep is obtained by removing the contribution from the solution x(t i−M−1 ) and adding the contribution of x(t i−1 ).The removal is accomplished in line 6 of Algorithm 3 by a rank-one update: By a cyclic permutation, we can move the zero column to the last column, , updating Z as in line 7 of Algorithm 3. Finally, the contribution of the latest solution is incorporated by adding the rank-one matrix x(t i−1 )z T M , where z M ∈ R m is a newly generated Gaussian random vector that is stored in the last row of Z.Under the (idealistic) assumption that all solutions are exactly computed (and hence deterministic), the described progressive updating procedure is mathematically equivalent to computing the randomized SVD from scratch.In particular, the error bound (4) continues to hold.
Lines 6-9 of Algorithm 3 require O(nm) operations.When using standard updating procedures for QR decomposition [13], line 11 has the same complexity.This compares favorably with the O(nM 2 ) operations needed by Algorithm 2 per timestep.
When performing the progressive update of Ω over many timesteps, one can encounter numerical issues due to numerical cancellation in the repeated subtraction and addition of contributions to the sketch matrix.To avoid this, the progressive update is carried out only for a fixed number of timesteps, after which a new random matrix Z is periodically generated and Ω is computed from scratch.

Convergence Analysis
We start our convergence analysis of the algorithms from the preceding section by considering analytical properties of the history matrix After reparametrization, we may assume without loss of generality that each of the past timesteps is contained in the interval [−1, 1]: For notational convenience, we define where x(t) satisfies the (parametrized) linear system that is, each entry of A and b is a scalar function on the interval [−1, 1].Indeed, for the convergence analysis, we assume that each linear system of the sequence in ( 1) is obtained by sampling the parametrized system in (7) in In many practical applications, like the one described in Section 5, the time dependence in (7) arises from time-dependent coefficients in the underlying PDEs.
Frequently, this dependence is real analytic, which prompts us to make the following smoothness assumption on A, b.

Compressibility of the solution time history
The effectiveness of POD-based algorithms relies on the compressibility of the solution history, that is, the columns of X can be well approximated by an m-dimensional subspace with m ≪ M .According to (2), this is equivalent to stating that the singular values of X decrease rapidly to zero.Indeed, this property is implied by Assumption 1 as shown by the following result, which was stated in [18] in the context of low-rank methods for solving parametrized linear systems.
Theorem 2 ([18, Theorem 2.4]).Under Assumption 1, the kth largest singular value σ k of the history matrix X(t) from (6) satisfies Combined with (2), Theorem 2 implies that the POD basis Q POD ∈ R n×m satisfies the error bound

Quality of prediction without compression
Algorithm 1 determines the initial guess s * for the next time step t i > t i−1 = 1 by solving the minimization problem In this section, we will assume, additionally to Assumption (1), that S = span(X(t)), that is, X(t) is not compressed.Our analysis focuses on uniform timesteps Note that the next timestep The following result shows how the quality of the initial guess rapidly improves (at a square root exponential rate, compared to the exponential rate of Theorem 2) as M , the number of previous time steps in the history, increases.Theorem 3.Under Assumption (1), the initial guess constructed by Algorithm 1 with S = span(X) satisfies the error bound

Proof of Theorem 3
The rest of this section is concerned with the proof of Theorem 3. We establish the result by making a connection to vector-valued polynomial extrapolation and extending results by Demanet and Townsend [8] on polynomial extrapolation to the vector-valued setting.
Let P R ⊂ R n [t] denote the subspace of vector-valued polynomials of length n and degree at most R for some R ≤ M − 1.We recall that any v ∈ P R takes the form v(t Equivalently, each entry of v is a (scalar) polynomial of degree at most R.In our analysis we consider vector-valued polynomials of the particular form for a vector-valued polynomial y(t) of length M .A key observation is that the evaluation of p in the next timestep t i satisfies p(t i ) ∈ span(X(t equi )) = S.
According to (8), s * minimizes the residual over S. Hence, the residual can only increase when we replace s * by p(t i ) in Thus, it remains to find a polynomial of the form ( 9) for which we can establish convergence of the extrapolation error p(t i ) − x(t i ) 2 .For this purpose, we will choose p R ∈ P R to be the least-squares approximation of the M function samples contained in X(t equi ): (11) We will represent the entries of p R in the Chebyshev polynomial basis: where c k,p ∈ R n and q k denotes the Chebyshev polynomial of degree k, that is, we can express (12) more compactly as p R (t) = C p q R (t).Thus, In view of (11), the matrix of coefficients C p is determined by minimizing

has full row rank and thus the solution of this least-squares problem is given by
In summary, we obtain that which is of the form (9) and thus contained in span(X(t equi )), as desired.
In order to analyze the convergence of p R (t), we relate it to Chebyshev polynomial interpolation of x.The following lemma follows from classical approximation theory, see, e.g., [18,Lemma 2.2].Lemma 4. Let q R (t) ∈ R R+1 be defined as in (13), containing the Chebyshev polynomials up to degree R.Under Assumption 1 there exists an approximation of the form Following the arguments in [8] for scalar functions, Lemma 4 allows us to estimate the extrapolation error for p R (t) if R ∼ √ M .
Proof.Letting x R be the polynomial from Lemma 4, we write To treat the second term in (15), first note that, by definition, we have where we used Lemma 4 in the last inequality.Applying, once more, Lemma 4 to the first term in (15) gives Because Inserted into (16), this gives The proof is completed by inserting the lower bound Using Theorem 5 with t = t i and inserting the result in (10), we have proven the statement of Theorem 3.

Optimality of the prediction with compression
When the matrix X(t) is compressed via POD (Algorithm 2) or the randomized range finder (Algorithm 3), the orthonormal basis Q ∈ R n×m used in Algorithm 1 spans a lower-dimensional subspace S ⊆ span(X).Corollary 6. Suppose that Algorithm 1 is used with an orthonormal basis satisfying (QQ T − I)X(t equi ) 2 ≤ ε for some tolerance ε > 0. Under Assumption 1, the initial guess s * constructed by the algorithm satisfies the error bound be the polynomial constructed in ( 14).
Using that s * satisfies the minimization problem ( 8) and QQ T x(t i ) ∈ S = span(Q), we obtain: The first term is bounded using Theorem 5 with t = t i .For the second term, we use the bound in ( 17) on q R (t M+1 ) to obtain ).The proof is completed using the lower bound (18) on xσ.

Numerical results: Test Case
To test the subspace acceleration algorithms proposed in Section 2, we first consider a simplified setting, an elliptic PDE with an explicitly given time-and space-dependent coefficient a(x, t) and source term g(x, t): We consider the domain Ω = [0, 1] 2 ⊂ R 2 and discretize ( 19) on a uniform twodimensional Cartesian grid using a centered finite difference scheme of order 4.This leads to a linear system for the vector of unknowns f (t), for which both the matrix and the right-hand side depend on t: We discretize the time variable on the interval [t 0 , t f ] with a uniform timestep ∆t on N t points, such that t f = t 0 +N t ∆t.Evaluating (20) in these N t instants, we obtain a sequence of linear systems of the same type as (1).We set a(x, t) = exp [−(x−0.5) 2 −(y−0.5) 2 ] cos(tx) + 2.1 and choose the righthand side g(x, t) such that f (x, t) = sin(4πyt) sin(15πxt) 1 + sin(15πxt) cos(3πyt) exp [(x−0.5) 2 +(y−0.5) 2 −0.25 2 ] is the exact solution of (19).The tests are performed using MATLAB 2023a on an M1 MacbookPro.We employ GMRES as iterative solver for the linear system, with tolerance 10 −7 and incomplete LU factorization as preconditioner.We start the simulations at t 0 = 2.3 s and perform N t = 200 timesteps.The results reported in Figure 1 use a spatial grid of dimension 100 × 100, leading to linear systems of size n = 10000.Different values of M , the number of previous solutions retained in the history matrix X, and m, the dimension of the reduced-order model, were tested.We found that the choices M = 20, m = 10 and M = 35, m = 20 lead to good performance for ∆t = 10 −5 and ∆t = 10 −3 , respectively.The baseline is (preconditioned) GMRES with the previous solution used as initial guess; the resulting number of iterations is indicated with the solid blue line ("Baseline") in Figure 1.This is compared to the number of iterations obtained by applying GMRES when Algorithm 1 is employed to compute the initial guess, in combination with both the POD basis in Algorithm 2 ("POD" in the graph) and the Randomized Range Finder in Algorithm 3 ("RAND" in the graph).For the Randomized Range Finder algorithm, the matrix Ω is computed from scratch only every 50 timesteps, while in the other timesteps is updated as described in Algorithm 3, resulting in a computationally efficient version of the algorithm.Both the POD and Randomized versions of the acceleration method give a remarkable gain in computational time with respect to the baseline.
When employing ∆t = 10 −5 , in Figure 1a, the number of iterations computed by the linear solver vanishes most of the time, since the initial residual computed with the new initial guess is already below the tolerance, set to 10 −7 in this case.It is worth noticing that the new randomized method gives acceleration comparable to the existing POD one, but it requires a much lower computational cost, as described in Section 2.
The results obtained for larger timesteps, in Figure 1b, are slightly worse, as expected, since it is less easy to predict new solutions using the previous ones when they are further apart in time.Nevertheless, the gain of the acceleration method is still visible, obtaining always less than half iterations with respect to the baseline and adding the solution of a reduced-order system of dimension m = 20 only, compared to the full solution of dimension 10000.The resulting advantage of the new method can indeed be observed in Figure 2, which compares the computational time needed by the solver using the baseline approach with the one obtained by using the new guess (this includes the time employed to compute the guess).The timings showed are the ones needed to produce the results in Figure 1.The time employed by the POD method has not been included since it is significantly higher than the baseline, as predicted by the analysis in Section 2.1.

Numerical results: Plasma Simulations
In this Section, we apply the subspace acceleration method to the numerical simulation of plasma turbulence in the outermost plasma region of tokamaks, where the plasma enters in contact with the surrounding external solid walls, resulting in strongly non-linear phenomena occurring on a large range of time and length scales.In this work, we consider GBS (Global Braginskii Solver) [12,23], a three-dimensional, flux-driven, two-fluid code developed for the simulation of the plasma dynamics in the boundary of a fusion device.GBS implements the Braginskii two-fluid model [3], which describes a quasi-neutral plasma through the conservation of density, momentum, and energy.This results in six coupled three-dimensional time-evolving non-linear equations which evolve the plasma dynamics in Ω, a 3D toroidal domain with rectangular poloidal cross section, as represented in Figure 3.The fluid equations are coupled with Maxwell equations, specifically Poisson and Ampére, elliptic equations for the electromagnetic variables of the plasma.In the limit considered here the elliptic equations reduce to a set of two-dimensional algebraic constraints decoupled along the toroidal direction, therefore to be satisfied independently on each poloidal plane.The differential equations are spatially discretized on a uniform Cartesian grid employing a finite difference method, resulting in a system of differential-algebraic equations of index one [15]: for each kth poloidal plane (21) where Y(f (t), x(t)) is a non-linear, 6-dimensional differential operator and are the vector of, respectively, the electromagnetic and fluid quantities solved for by GBS, where the solutions of all the N Z poloidal planes are stacked together.More precisely, the time evolution of the fluid variables, f , is coupled with the set of linear systems A k (f (t))x k (t) = b k (f (t)) which result from the discretization of Maxwell equations.Indeed, the matrix A k ∈ R NRNZ ×NRNZ and right-hand side b k ∈ R NRNZ depend on time through f .In GBS, system (21) is integrated using a Runge-Kutta scheme of order four, on the discrete times {t i } Nt i=1 , with step-size ∆t.Given f i and x i , the value of f and x at time t i , the computation of f i+1 , requires performing three intermediate substeps where the quantities f i+1,j for j = 1, 2, 3 are computed.To guarantee the consistency and convergence of the Runge-Kutta integration method [15], the algebraic constraints are solved at every substep, computing x i+1,j k for j = 1, 2, 3 and for each k−th poloidal plane.As a consequence, the linear systems A k (f (t))x k (t) = b k (f (t)) are assembled and solved four times for each of the N ϕ poloidal planes, to advance the full system (21) by one timestep.Since the timestep ∆t is constrained to be small from the stiff nature of the GBS model, the solution of the linear systems is among the most computationally expensive part of GBS simulations.
In GBS, the linear system is solved using GMRES, with the algebraic multigrid preconditioner boomerAMG from the HYPRE library [9], a choice motivated by previous investigations [12].The subspace acceleration algorithm proposed in Section 2 is implemented in the GBS code and, given the results shown in Section 4, the randomized version of the algorithm is chosen.The results reported are obtained from GBS simulations on one computing node.The poloidal planes of the computational domain are distributed among 16 cores, specifically of type Intel(R) Core i7-10700F CPU at 2.90GHz.GBS is implemented in Fortran 90, and relies on the PETSc library [2] for the linear solver and Intel MPI 19.1 for the parallelization.
We consider the simulation setting described in [12], taking as initial conditions the results of a simulation in a turbulent state.We use a Cartesian grid of size of N R = 150, N Z = 300 and N ϕ = 64, with additional 4 ghost points in the Z and R directions.Therefore, the imposed algebraic constraints result in 64 sequences of linear systems of dimension N R N Z × N R N Z = 46816 × 46816.The timestep employed is ∆t = 0.7 × 10 −5 .The sequence of linear systems we consider represents the solution of the Poisson equation on one fixed poloidal plane, but the same considerations apply to the discretization of Ampére equation.In Figure 4a the number of iterations obtained with the method proposed in Section 2, denoted as "RAND" is compared with the ones obtained using the previous step solution as initial guess, depicted in blue as "Baseline".We notice that, employing the acceleration method, the number of GMRES iterations needed for each solution of the linear system is reduced by a factor 2.9, on average, at the cost of computing a solution of an m × m reduced-order system.In Figure 4b the wall clock time required for the solution of the systems is shown.The baseline approach is compared to the accelerated method, where we also take into account the cost of computing the initial guess.Thanks to the randomized method employed, the process of generating the guess is fast enough to provide a time speed up of a factor of 6.5 per iteration.
The employed values of M = 15, the number of previous solutions retained, and m = 10, the dimension of the reduced-order model, are the ones found to give a good balance between the decrease in the number of iterations and the computational cost of the reduced-order model.In Table 1 the results for different values of M and m are reported.It is worth noticing that an average number of GMRES iterations per timestep smaller than one implies that often the initial residual obtained with the initial guess is below the tolerance set for the solver.It is possible to notice that higher values of m lead to very small number of iterations, but the overall time speedup is reduced since the computation of the guess becomes more expensive.

M m
Average The iteration and time speedups are computed on the total of 180 linear systems, with respect to the baseline, that has an average number of 5 GMRES iterations and an average time per timestep of 0.0828 s.The highlighted row corresponds to the best result obtained in terms of time speedup.

Conclusions
In this paper, we propose a novel approach for accelerating the solution of a sequence of large-scale linear systems that arises from, e.g., the discretization of time-dependent PDEs.Our method generates an initial guess from the solution of a reduced-order model, obtained by extracting relevant components of previously computed solutions using dimensionality reduction techniques.Starting from an existing POD-like approach, we accelerate the process by employing a randomized algorithm.A convergence analysis is performed, which applies to both approaches, POD and the randomized algorithm and shows how the accuracy of the method increases with the history size.A test case displays how POD leads to a noticeable decrease in the number of iterations, but at the same time a nearly equal decrease is achieved by the cheaper randomized method, that leads to a time speedup per iteration of a factor 9. In real applications such as the plasma simulations described in Section 5, the speedup is more mod-est, given the stiff nature of the problem which constrains the timestep of the explicit integration method to be very small, but still practically relevant.

Figure 2 :
Figure 2: Computational time per timestep corresponding to Figure 1a and Figure 1b.The average speedup per iteration of the randomized method with respect to the baseline is a factor 9 for ∆t = 10 −5 and a factor 10 for ∆t = 10 −3 .

Figure 3 :
Figure 3: GBS computational domain.The toroidal direction is along ϕ, the radial direction is along R, and the vertical direction is along Z.The domain consists of N ϕ rectangular poloidal planes, each discretized on a N R × N Z Cartesian grid.

Figure 4 :
Figure 4: Performance of the algorithm applied to the solution of Poisson equation in GBS simulations.The time for the RAND algorithm is on average approximately one fourth of the time for the baseline.

Table 1 :
GBS simulations result corresponding to different values of M and m.