# Vector extrapolation applied to algebraic Riccati equations arising in transport theory.

1. Introduction. An important problem that arises in different areas of science and engineering is the computation of limits of sequences of vectors [x.sub.0], [x.sub.1],..., where [x.sub.i] are complex N-vectors with N large. Such sequences may result from iterative methods or perturbation techniques and may converge very slowly. This is the case, for example, when they result from finite-difference or finite-element discretizations of continuum problems, where rates of convergence become worse as the mesh size decreases. Vector extrapolation methods an be applied to such vector sequences. These methods transform a sequence of vectors gen rated by some process to a new sequence which converges faster, without requiring explicit knowledge of the sequence generator.

Our concern is a special kind of nonsymmetric algebraic Riccati equation (NARE) that irises in transport theory. We will apply a polynomial vector extrapolation method [7, 34] namely RRE, to a vector sequence produced by an iterative method for computing the mini mal positive solution of NARE.

This paper is organized as follows. In Section 2, we introduce a specific kind of NARE hat arises in transport theory. We will review some iterative methods which have been proven o be efficient for solving these kinds of equations. We will focus on the method of Lin  and we will make a small modification, similar to that of NBGS in  to accelerate it. Section 3 is devoted to polynomial vector extrapolation methods, namely the reduced rank extrapolation (RRE). RRE is an efficient convergence accelerator and a short review of this method will be given. An application of the RRE method to the iterative method of Lin  s considered followed by a numerical example in Section 4 which shows the effectiveness o our approach.

2. NARE. Nonsymmetric algebraic Riccati equations appear in transport theory, where a variation of the usual one-group neutron transport equation [3, 9, 12, 22] is formulated as

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [phi](x) is the neutron flux, [alpha] (0 [less than or equal to] [alpha] < 1) is an angular shift, and c is the average of the total number of particles emerging from a collision, which is assumed to be conserved, i.e., 0 < c [less than or equal to] 1.

The scattering function L : [-[alpha], 1] x [[alpha], 1] [right arrow] R for particle transport in the half-space can be derived from (2.1) and satisfies the integro-differential equation

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where ([mu], y) [member of] [-[alpha], 1] x [[alpha], 1] and ([alpha], c) [member of] [(0,1).sup.2] are pairs of nonnegative constants; see the appendix in .

When c = 0 or [alpha] = 1, the integro-differential equation (2.2) has a trivial solution . However, when 0 [less than or equal to] [alpha] < 1 and 0 < c [less than or equal to] 1, it has a unique, positive, uniformly bounded, and globally defined solution .

Discretization of the integro-differential equation (2.2) by a numerical quadrature formula on the interval [0,1] yields an algebraic matrix Riccati equation.

Nonsymmetric algebraic Riccati equations (NARE) are quadratic matrix equations of the general form

(2.3) XCX - XD - AX + B = 0.

where the coefficient matrices A, B, C, and D [member of] [R.sup.nxn].

Let [{[w.sub.i]}.sup.n.sub.i=1] and [{[c.sub.i]}.sup.n.sub.i=1] denote the sets of nodes and weights, respectively, of the specific quadrature rule that is used on the interval [0, 1]. This quadrature rule is obtained by dividing the interval into n/4 subintervals of equal length and applying a Gauss-Legendre quadrature with 4 nodes to each subinterval. These typically satisfy:

(2.4) [c.sub.1],...,[c.sub.n] > 0, [n.summation over (i=1)] [c.sub.i] = 1, and 1 > [w.sub.i] > [w.sub.2] > ... > [w.sub.n] > 0.

We are interested in a special case of NARE where

(2.5) A = [DELTA] - [eq.sup.T], B = [ee.sup.T], C = [qq.sup.T], D = [LAMBDA] - [qe.sup.T],

and

(2.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where diag(*) applied to a vector denotes a square diagonal matrix with the elements of the vector on its diagonal. Here, 0 [less than or equal to] [alpha] < 1 and 0 < c [less than or equal to] 1. As a consequence of (2.4), we have

0 < [[delta].sub.1] < [[delta].sub.2] < ... [[delta].sub.n] and 0 < [[gamma].sub.1] < [[gamma].sub.2] < ... < [[gamma].sub.n].

[[gamma].sub.i] [less than or equal to] [[delta].sub.i] for a [less than or equal to] 0,

[[gamma].sub.i] > [[delta].sub.i] for [alpha] [not equal to] 0,

for i = 1, 2,..., n.

The term nonsymmetric distinguishes this case from the widely studied continuous-time algebraic Riccati equations (CARE), defined by the quadratic matrix equation

XCX - AX - [XA.sup.T] + B = 0,

where B and C are symmetric. For a comprehensive analysis of CARE; see [5, 23, 27].

Nonsymmetric algebraic Riccati equations arise in many fields such as transport theory when dealing with particle transfer (or radiative transfer), Wiener-Hopf factorization of Markov chains, nuclear physics, applied probability, engineering, control theory, etc. In this paper, we are interested in special Riccati equations that arise in transport theory with A, B, C, and D given as in (2.5).

The existence of nonnegative solutions of (2.3) was shown by Juang  using degree theory. It is shown [18, 21] that the Riccati equation (2.3) has two entry-wise positive solutions X = [[x.sub.ij]] and Y = [[y.sub.ij]] in Rnxn which satisfy X [less than or equal to] Y, where we use the notation X [less than or equal to] Y if [x.sub.ij] [less than or equal to] [y.sub.ij] for all i, j = 1, ...,n. In applications from transport theory, only the smaller one of the two positive solutions is of interest and is physically meaningful. Some iterative procedures  were developed to find both nonnegative solutions.

2.1. Obtaining the minimal positive solution of NARE by a vector equation. It has been shown by Lu  that the minimal positive solution of a NARE can be obtained by computing the minimal positive solution of a vector equation. This vector equation is derived from a special form of solutions of the Riccati equation and by exploiting the special structure of the coefficient matrices of the Riccati equation.

Rewrite the NARE (2.3)-(2.5) as

[DELTA]X + XD = (Xq + e)([q.sup.T] X + [e.sup.T])

and let

(2.7) u = Xq + e and [v.sup.T] = [q.sup.T] X + [e.sup.T].

It has been shown in [18, 21] that any solution of (2.3) must be of the form

(2.8) X = T [omicron] ([uv.sup.T]) = ([uv.sup.T]) [omicron] T,

where o denotes the Hadamard product defined by A [omicron] B = [[a.sub.ij][b.sub.ij]] for any two matrices A = [[a.sub.ij]-] and B = [[b.sub.ij]], T is the Cauchy matrix, T = [[t.sub.i,j]] = + [1/([[delta].sub.i][gamma].sub.j])], X = [[x.sub.ij]], u = [[[.sub.u1], [u.sub.2],...,[u.sub.n]].sup.T] and v = [[[v.sub.1], [v.sub.2],..., [v.sub.n]].sup.T].

Finding the minimal positive solution of (2.3) requires finding two positive vectors u and v in (2.8). Substituting (2.8) into (2.7) gives the vector equation

(29) u = u [omicron] (Pv) + e;

v = v [omicron] (Qu) + e,

where

(2.10) P = [[P.sub.ij]] = [[q.sub.j]/[[delta].sub.i] + [[gamma].sub.j]] = T diag(q)

and

(2.11) Q = [[Q.sub.ij]] = [[q.sub.j]/[[delta].sub.i] + [[gamma].sub.j]] = [T.sup.T] diag(q)

The minimal positive solution of (2.3) can be obtained by computing the minimal positive solution of the vector equation (2.9). This vector equation can be expressed as

f(u, v) := u - u [omicron] (Pv) - e = 0,

g(u, v) := v - v [omicron] (Qu) - e = 0

which is a system of nonlinear equations having the vector pair (u, v) as a positive solution. Based on the above, Lu  defined an iterative method for the solution of (2.9).

2.1.1. The iterative method of Lu. Lu defined an iterative scheme to find the positive vectors u and v,

(2.12)

[u.sup.(k+1)] = [u.sup.(k)] [omicron] ([Pv.sup.(k)])+ e,

[v.sup.(k+1)] = [v.sup.(k)] [omicron] ([Qu.sup.(k)]) + e, for k = 0,1,...,

[u.sup.(0), [v.sup.(0)]) = (0,0),

satisfying a certain stopping criterion. For all 0 [less than or equal to] [alpha] < 1 and 0 < c [less than or equal to] 1, the sequence {([u.sup.(k)], [v.sup.(k)])} defined by (2.12) is strictly monotonically increasing, bounded above, and thus converging; see  for the proof. Each iteration costs approximately 4[n.sup.2] flops.

The minimal positive solution of (2.3) can be computed by

X* = T [omicron] (u*[(v*).sup.T]),

where (u*, v*) is the limit of ([u.sup.(k)], [v.sup.(k)]) defined by the iteration scheme (2.12).

2.1.2. A modified iterative method. Bao et al.  proposed a modified version of the iterative method of Lu . They noticed that since [u.sup.(k+1)] is obtained before [v.sup.(k+1)], it should be a better approximation to u* than [u.sup.(k)]. Consequently, [u.sup.(k)] is replaced with [u.sup.(k+1)] in the equation for [v.sup.(k+1)] of the iteration scheme (2.12).

The modified iteration is as follows

(2.13)

[[??].sup.(k+1)] = [[??].sup.(k)] [omicron] ([P[??].sup.(k)])+ e,

[[??].sup.(k+1)] = [[??].sup.(k)] [omicron] ([Q[??].sup.(k+1)]) + e, for k = 0,1,...,

[[??].sup.(0)] = [[??].sup.(0)] = 0.

The monotonic convergence of the modified iteration scheme (2.13) is illustrated in . Also, it is shown how the minimal positive solution X* of the NARE (2.3) can be computed from

X* = T [omicron] [(u* (v*).sup.T]) ,

where (u*, v*) is the limit of ([??](k), [??](k)) defined by the modified iteration (2.13).

Since it was also shown that [u.sup.(k)] < [[??].sup.(k)] and [v.sup.(k)] < [[??].sup.(k)] for k [greater than or equal to] 3, ([[??].sup.(k)], [[??].sup.(k)]) has the same limit (u*, v*) as ([u.sup.(k)], [v.sup.(k)]). Also, this proves how the modified iteration scheme

(2.13) is more efficient than (2.12). The cost of every iteration step is approximately 4[n.sup.2] flops.

Two accelerated variants of the iterative scheme (2.12) of Lu, proposed in , are the nonlinear block Jacobi (NBJ) iteration scheme

(2.14) [u.sup.(fc+1)] = [u.sup.(k+1)] [omicron] ([Pv.sup.(k)]) + e, [v.sup.(k+1)] = [v.sup.(k+1)] o ([Qu.sup.(k)]) + e, for k = 0,1,...,

and the nonlinear block Gauss-Seidel (NBGS) iteration scheme

(2.15) [u.sup.(k+1)] = [u.sup.(k+1)] [omicron] ([Pv.sup.(k)]) + e, [v.sup.(k+1)] = [v.sup.(k+1)] [omicron] ([Qu.sup.(k+1)]) + e, for k = 0, 1,.. ..

For both NBJ and NBGS, the sequence {([u.sup.(k)], [v.sup.(k)])} is shown in  to be strictly monotonically increasing and convergent to the minimal positive solution (u*,v*) of the vector equations (2.9). NBJ and NBGS have the same computational costs at every iteration step (approximately 4[n.sub.2] flops). Both are effective solvers of NAREs arising in transport theory, but NBGS is more efficient than NBJ in applications. In particular, in terms of the asymptotic rates of convergence, the NBGS method is twice as fast as the NBJ method; see [15, Theorem 5] which explains the numerical results presented in , where the number of iterations required for NBGS is half of that for NBJ.

In conclusion, the four fixed-point iterations (2.12), (2.13), (2.14), and (2.15) are easy to use and share the same low complexity at every iteration. However, NBGS was proved in  to be the fastest among these methods in terms of the asymptotic rate of convergence when ([alpha], c) [not equal to] (0,1), although being sublinear when ([alpha], c) = (0,1). The sublinear convergence, which takes place when the Jacobian at the required solution is singular, is transformed into a quadratic convergence by means of a Newton method proposed in .

2.1.3. The Newton method. Guo and Laub proposed in  an application of the Newton method to the Riccati equation (2.3). This method has an order of complexity O([n.sup.3]) per step and yields quadratic convergence. Lu  instead proposed an algorithm that consists of iteration (2.12) combined with a Newton iteration. This scheme is simple and more efficient than applying the Newton method directly to (2.3). Another approach was proposed by Bini et al.  who apply a fast Newton method to the NARE (2.3). It consists of an algorithm which performs the Newton iteration in O([n.sup.2]) flops. This approach relies on a modification of the fast LU factorization algorithm for Cauchy-like matrices proposed by Gohberg et al. . The same idea reduces the cost of Newton's algorithm proposed by Lu  from O([n.sup.3]) to O([n.sup.2]) while preserving the quadratic convergence in the generic case.

2.1.4. The iterative method of Lin. This method was proposed in . Let w = [[[u.sup.T], [v.sup.T]].sup.T] and reformulate the vector equation (2.9) as

(2.16) f(w) := w - w [omicron] Mw - e = 0,

where,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The Jacobian matrix J(f, w) of f(w) for any w [member of] [R.sup.2n] is given by

J(f,w)= [I.sub.2n] - N(w),

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Lin constructed a class of iterative methods to solve the vector equation (2.16) based on the following iterative scheme

(2.17) [w.sup.(k+1)] := [w.sup.(k)] - [T.sup.-1.sub.k] f([w.sup.(k)]), for k = 0,1, 2,...,

where [T.sub.k] is chosen to approximate J(f, [w.sup.(k)]) and [T.sub.k] [greater than or equal to] J(f, [w.sup.(k)]). Note that when [T.sub.k] = J(f, [w.sup.(k)]), the Newton method results. A disadvantage of the Newton method is that it requires an LU factorization of the Jacobian matrix J(f, [w.sup.(k)]) at each iteration step to obtain the new approximation to w*. This costs O([n.sup.3]) operations per step.

For any [w.sup.(k)] = [[[([u.sup.(k)]).sup.T], ([v.sup.(k)])T].sup.T] [member of] [R.sup.2n], Lin  proposed the choice

(2.18) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Substituting (2 18) into (2 17) gives

(2.19) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Convergence of the iterative method (2 19) was examined in  It was shown that the convergence is sublinear as ([alpha], c) tends to (0,1) and linear when ([alpha], c) = (0,1). Since ([I.sub.n] - diag([Pv.sup.(k)])) and ([I.sub.n] - diag(Q[u.sup.(k)])) are diagonal matrices, the computational cost of each iteration step is approximately 4[n.sup.2] flops. A numerical comparison shows how the iterative method of Lin converges faster than the modified iterative method of Lu (2.13) and the Newton method.

2.1.5. A modification of the iterative method of Lin. A modification of the iterative scheme (2.19) that is an analog of the NBGS method  is proposed. In the iterative scheme (2.19), [u.sup.(k+1)] is computed before [v.sup.(k+1)]. Therefore, [u.sup.(k+1)] should be a better approximation of u* than [u.sup.(k)]. Replacing [u.sup.(k)] by [u.sup.(k+1)] in the equation for [v.sup.(k+1)] in (2.19) leads to the modified iteration scheme

(2.20) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For a matrix G [member of] [R.sup.nxn], we define the following function

[PHI]G : [R.sup.n] [right arrow] [R.sup.n],

t [right arrow] [([I.sub.n] - diag(Gt)).sup.-1] e.

Then iteration (2.20) can be expressed as

(221) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

LEMMA 2.1. Let (u*,v*) be the minimal positive solution of (2.20). Then, u* > e and v* > e.

Proof. Let (u*, v*) be the minimal positive solution of (2.20). Then

u* = [([I.sub.n] - diag(Pv*)).sup.-1] e > e > 0,

v* = [([I.sub.n] - diag(Qu*)).sup.-1] e > e > 0,

since P and Q are positive matrices.

THEOREM 2.2. Given the sequence of vectors generated by (2.20) with initial vector ([u.sup.(0)],[v.sup.(0)]) = (0,0), let (u*, v*) be the minimal positive solution of (2.20). Then the sequence {([u.sup.(k)], [v.sup.(k)])} is strictly monotonically increasing, bounded above and thus convergent.

Proof. To show this, we have to prove component-wise that

(i) 0 [less than or equal to] [u.sup.(k)] < [u.sup.(k+1)] < u* and 0 [less than or equal to] [v.sup.(k)] < [v.sup.(k+1)] < v*, k [less than or equal to] 1;

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We start by proving (i) by induction. Let ([u.sup.(0)], [v.sup.(0)]) = (0, 0) be a starting point for this method. For k = 0, using (2.20) and Lemma 2.1, we get

[u.sup.(0)] = 0 < e = u(1) < u* .

Also,

[v.sup.(1)] = [[[I.sub.n] - diag([Qu.sup.(1)])].sup.-1]e = [[[I.sub.n] - diag(Qe)].sup.-1]e < [[[I.sub.n] - diag(Qv*)].sup.-1]e

and, therefore,

[v.sup.(0)] = 0 < e = [v.sup.(1)] < v* .

This shows that (i) holds for k = 0. Now, suppose that (i) holds for a positive integer k, i.e., we have

0 [less than or equal to] [u.sup.(k)] < [u.sup.(k+1)] < u* and 0 [less than or equal to] [v.sup.(k)] < [v.sup.(k+1)] < v*, for k [greater than or equal to] 1.

Using (2.20) and the fact that [v.sup.(k)] < [v.sup.(k+1)], we have

[[I.sub.n] - diag([Pv.sup.(k+1)])][u.sup.(k+2)] = e = [[I.sub.n] - diag([Pv.sup.(k)])][u.sup.(k+1)] > [[I.sub.n] - diag([Pv.sup.(k+1)])][u.sup.(k+1)].

Since [I.sub.n] - diag([Pv.sup.(k+1)]) > [I.sub.n] - diag(Pv*) > 0, using Lemma 2.1 and the fact that [v.sup.(k+1)] < v*, we get

0 [less than or equal to] [u.sup.(k+1)] <u(k+2).

Also,

[[I.sub.n] - diag([Pv.sup.(k+1)])][u.sup.(k+2)] = e = [[I.sub.n] - diag(Pv*)]u* < [[I.sub.n] - diag([Pv.sup.(k+1)])]u*.

Therefore, [u.sup.(k+2)] < u* holds, and consequently (i) holds for (k + 1). In conclusion, (i) is shown by induction.

The proof of (ii) can be done via (i) which provides the existence of two positive vectors ([??]*, [??]*), with 0 < [??]* < u* and 0 < [??]* < v*, satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[??]* = [[PHI].sub.P] ([??]*),

[??]* = [[PHI].sub.Q] ([??]*),

i.e., ([??]*, [??]*) is a positive solution of (2.20). Due to the minimal property of (u*, v*) and the comparative property of ([??]*, [??]*) with (u*, v*), it must hold that [??]* = u* and [??]* = v*.

We consider the case when [alpha] = 0 and c = 1, which is referred to as the critical case. Looking back at equations (2.6), it can be seen that when [alpha] = 0, we have [J.sub.i] = [[gamma].sub.i] = 1/[w.sub.i] and consequently [DELTA] = [LAMBDA]. We have T = [[t.sub.i,j]] = [1/[[delta].sub.1] + [[gamma].sub.j]] and since [[gamma].sub.i] = 7j, it follows that T = [1/[delta]+[[delta].sub.j]] and T is symmetric (T = [T.sup.T]). In view of (2.10) and (2.11), we get P = Q and consequently [[PHI].sub.P] = [[PHI].sub.Q].

Let K and L be two matrices in Rnxn and define the following transformation

[[PHI].sub.K,L] :[R.sup.n] [right arrow] [R.sup.n], t [right arrow] [[PHI].sub.K] ([[PHI].sub.L](t))

Then iteration (2.21) becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

And for [[PHI].sub.P] = [[PHI].sub.Q],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Starting with equal initial vectors [u.sup.(0)] = [v.sup.(0)] = 0, it can be easily seen by induction that the sequences [{[u.sup.(k)]}.sub.k[member of]N] and [{[v.sup.(k)]}.sub.k[member of]N] are equal. This implies that it is enough to compute one of the vector sequences, for example [{[u.sup.(k)]}.sub.k[member of]N], and the solution of NARE then can be calculated by

X = T [omicron] (u*(u*)t),

where u* denotes the limit of [u.sup.(k)]. Thus, only half of the computational work of (2.21) is needed.

Now, we consider the calculation of the Jacobian matrix of the modified iteration (2.21).

Let

[PHI] : [R.sup.n] x [R.sup.n] [right arrow] [R.sup.n] x [R.sup.n] ,

(u,v) [right arrow] ([[PHI].sub.P,Q](u), [[PHI].sub.Q,P](v))

be the transformation describing the iteration (2.21). Then the iterative scheme (2.21) can be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The Jacobian matrix is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

At the solution (u*, v*), the Jacobian is of the form

(2.22) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

because [PHI]Q(u*) = v* and [[PHI].sub.P](v*) = u*. We have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Therefore,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Define

K : [R.sup.nxn] x [R.sup.n] [right arrow] [R.sup.nxn],

(Y,t) [right arrow] [{[([I.sub.n] - diag(Yt)).sup.-1]}.sup.2]Y.

Then

(2.23) J ([[PHI].sub.P], v*) = K (P,v*).

Similarly, we have

(2.24) J ([[PHI].sub.Q],u*) = K (Q,u*).

Substituting (2.23) and (2.24) into (2.22), the Jacobian matrix at the solution becomes

(2.25) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

3. Vector extrapolation methods [7, 34]. The most popular vector extrapolation methods are the minimal polynomial extrapolation (MPE) method by Cabay and Jackson , the reduced rank extrapolation (RRE) method by Eddy  and Mesina , the modified minimal polynomial extrapolation (MMPE) method by Sidi et al. , Brezinski , and Pugachev , the topological e-algorithm (TEA) by Brezinski , and the vector e-algorithm (VEA) by Wynn [35, 36]. These methods do not require explicit knowledge of how the sequence is generated and can be applied to the solution of linear and nonlinear systems of equations. Several recursive algorithms for implementing these methods are presented in [6, 10, 11, 16]. A numerical comparison of vector extrapolation methods and their applications can be found in [17, 34].

We are interested in applying the RRE method to the iteration (2.20) to accelerate its convergence. See Table 4.1 for a comparison between RRE, MPE, and MMPE which explains the choice of using RRE. Throughout this paper, we denote by (*, *) the Euclidean inner product in RN and by [parallel] * [parallel] the corresponding norm.

3.1. Review of RRE . Let [{[s.sup.(k)]}.sub.k[member of]N] be a sequence of vectors in [R.sup.N]. Define the first and the second forward differences of [s.sup.(k)] by

[DELTA][s.sup.(k)] = [s.sup.(k+1)] - [s.sup.(k)] and [[DELTA].sup.2][s.sup.(k)] = [DELTA][s.sup.(k+1)] - [DELTA][s.sup.(k)], k = 0,1,...,

When applied to the vector sequence [{[s.sup.(k)]}.sub.k[member of]N], the RRE method produces an approximation [t.sup.(k)] of the limit or the antilimit of [{[s.sup.(k)]}.sub.k[member of]N]; see . The approximation [t.sup.(k)] is defined by

(3.1) [t.sup.(k)] = [k.summation over (j=0)] [eta][n.sup.(k).sub.j][s.sup.(j)],

where

(3.2) [k.summation over (j=0)] [[eta].sup.(k).sub.j] = 1

and

(3.3) [k.summation over (j=0)] [[beta].sub.i,j][eta][n.sup.(k).sub.j] = 0, for i = 0, 1, ..., k - 1.

The scalars [[beta].sub.i,j] [member of] R are defined by

[[beta].sub.i,j] = ([[DELTA].sup.2][s.sup.(1)], [DELTA][s.sup.(j)]),

for i = 0, 1, ..., k - 1 and j = 0, 1, ..., k.

It follows from (3.1), (3.2), and (3.3) that [t.sup.(k)] can be expressed as the ratio of two determinants

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Using the following notation

[[DETAL].sub.i][S.sup.(k)] = [[[DETAL].sup.i] [s.sup.(0)], ...,[[DETAL].sup.I][s.sup.(k-1)]]

and Schur complements, [t.sup.(k)] can be written as

[t.sup.(k)] = [s.sup.(0)] - [DELTA][S.sup.(k)][([[DETAL].sup.2][S.sup.(k)]).sup.+] [DETAL][s.sup.(0)],

where ([[DELTA].sup.2][S.sup.(k)])+ denotes the Moore-Penrose generalized inverse of A2S(k) defined by

[([[DELTA].sup.2][S.sup.(k)]).sup.+] = [([([[DELTA].sup.2][S.sup.(k)]).sup.T] [[DELTA].sup.2][S.sup.(k)]).sup.-1] [([[DELTA].sup.2][S.sup.(k)]).sup.T].

It is clear that [t.sup.(k)] exists and is unique if and only if det([([[DELTA].sup.2][S.sup.(k)]).sup.T][[DELTA].sup.2][S.sup.(k)]) [not equal to] 0; see . We assume this to be the case. The computation of the approximation [t.sup.(k)] can be carried out by one of the recursive algorithms in .

To give an estimate for the residual norm for nonlinear problems, introduce the new transformation

[[??].sup.(k)] = [k.summation over (j=0)] [[eta].sup.(k).sub.j][s.sup.(j+1)].

We use the generalized residual of [t.sup.(k)] defined in . It is given by

[??]([t.sup.(k)]) = [[??].sup.(k)] - [t.sup.(k)] = [k.summation over (j=0)] [[eta].sub.j][DELTA][s.sup.(j)]

and can be expressed as

[??]([t.sup.(k)]) = [DELTA][s.sup.(0)] - [[DELTA].sup.2][S.sup.(k)]([[DELTA].sup.2][S.sup.(k)]) + [DELTA][s.sup.(0)].

Note that the generalized residual [??]([t.sup.(k)]) is obtained by projecting [DELTA][s.sup.(k)] orthogonally onto the subspace spanned by [[DELTA].sup.2][s.sup.(0)],..., [[DELTA].sup.2][s.sup.(k-1)]. It is an approximate residual unless the sequence is generated linearly, in which case it is the true residual. Therefore, a stopping criterion can be based on [parallel][??]([t.sup.(k)])[parallel].

3.2. Implementation of the RRE method. We follow the description given by Sidi in . An important feature of this method is that it proceeds through the solution of least squares problems by QR factorization.

Introduce the following notation

[DELTA][S.sup.(k+1)] = [[DELTA][s.sup.(0)], ...,[DELTA][s.sup.(k)]] and [[eta].sup.(k)] = [([[eta].sub.0.sup.k]), ...,[[eta].sub.k.sup.(k)]).sup.T].

In view of (3.1), (3.2), and (3.3), [[eta].sub.j.sup.(k)] can be determined by solving the over determined linear system

[DELTA][S.sup.(k+1)][[eta].sup.(k)] = 0

by the least-squares method subject to the constraint [[SIGMA].sup.k.sub.j=0][[eta].sub.j.sup.(k)] = 1. This leads to minimizing the positive definite quadratic form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

subject to the same constraint.

Set [d.sup.(k)] = [[d.sup.(k).sub.0] ,...,[[d.sup.(k).sub.k])].sup.T], [lambda] = [([[SIGMA].sup.k.sub.i=0] [d.sup.(k).sub.i]).sup.-1], and [[eta].sup.(k)] = [lambda][d.sup.(k)] ([[eta].sup.(k).sub.i] = [lambda][d.sup.(k)i)].

Then [[eta].sup.(k)] can be computed by solving the linear system of equations

(3.4) [([DELTA][S.sup.(k+1)]).sup.T] [DELTA][S.sup.(k+1)][d.sup.(k)] = e.

Assume that [DELTA][S.sup.(k+1)] has full rank. Then it has a QR factorization [DELTA][S.sup.(k+1)] = [Q.sup.(k)][R.sup.(k)].

This leads to another form of the linear system (3.4),

[([R.sup.(k)]).sup.T] [R.sup.(k)] [d.sup.(k)] = e.

Finally, the approximation [t.sup.(k)] can be expressed as

[t.sup.(k)] = [s.sup.(0)] + [Q.sup.(k-1)] ([R.sup.(k-1)][[xi].sup.(k)]),

where

[[xi].sup.(k)] = [[[xi].sup.(k).sub.0], [[xi].sup.(k)1], ...,[[[xi].sup.(k).sub.k-1]].sup.T], [[xi].sup.(k).sub.0] = 1 - [[eta].sup.(k).sub.0], and [[xi].sup.(k).sub.j] = [[xi].sup.(k).sub.j-1] - [[eta].sup.(k).sub.j]

for j = 1, ...,k - 1.

Another expression of [t.sup.(k)] is given by

(3.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then, using (3.4) and (3.5), the generalized residual [??] ([t.sup.(k)]) can be written as

[??]([t.sup.(k)]) = [k.summation over i=0] [n.sup.(k).sub.i] [DELTA][s.sup.(i)] = [DELTA][S.sup.(k+1)][[eta].sup.(k)].

Note that the QR factorization of [DELTA][S.sup.(k+1)] is formed by appending one additional column to [Q.sup.(k-1)] to obtain [Q.sup.(k)], and a corresponding column to [R.sup.(k-1)] to obtain [R.sup.(k)]. This QR factorization can be computed inexpensively by applying the modified Gram-Schmidt process (MGS) to the vectors [s.sup.(0)], [s.sup.(1)], ...,[s.sup.(k+1)]; see .

The computations for the RRE method are described by Algorithm 1. This algorithm becomes increasingly expensive as k increases, because the work requirement grows quadratically with the number of iteration steps k. The storage requirement grows linearly with k. To avoid this, the RRE algorithm should be restarted periodically every f steps for some integer r > 1 .
```
Algorithm 1 Basic RRE algorithm.

Step 0   Input:Vectors [s.sup.(0)], [s.sup.(1)], ..., [s.sup.(k+1)].
Step 1   Compute [DELTA][s.sup.(i)] = [s.sup.(i+1)] - [s.sup.(i)], for
i = 0,1, ...,k.

Set [DELTA][S.sup.(k+1)] = [[DELTA]][s.sup.(0)], [DELTA]
[s.sup.(1)], ...,[DELTA][s.sup.(k)]].

Compute the QR factorization of AS(k+1), namely, AS(k+1)
= Q(k)R(k).
Step 2   Solve the linear system
[([R.sup.(k])).sup.T][R.sup.(k)][d.sup.(k)] = e; where
[d.sup.(k)] = [[d.sup.k).sub.0], [d.sup.(k).sub.1), ...
,[d.sup.(k).sub.k]].sup.T] and e = [[1, 1, ...,1].sup.T].

(This amounts to solving two upper and lower triangular
systems.)

Set [lambda] = [([k.summation over i=0] [d.sup.(k).sub.i])
.sup.-1], [lambda] [member of] [R.sup.+].

Set [[eta].sup.(k).sub.i] = [lambda][d.sup.(k).sub.i],
for i = 0,1, ...,k.

Step 3   Compute [[xi].sup.(k)] = [[[xi].sup.(k)], [[[xi].sup.
(k).sub.1], ...,[[xi].sup.(k).sub.k-1]].sup.T] where
[[xi].sup.(k).sub.0] = 1 - [[eta].sup.(k).sub.0]
and [[xi].sup.(k).sub.j] = [xi].sup.(k).sub.j-1]
[[eta].sup.(k).sub.j], 1 [less than or equal to]
j [less than or equal to] k - 1.

Compute [t.sup.(k)] by :[t.sup.(k) = [s.sup.(0)] +
[Q.sup.(k-1)]([R.sup.(k-1)][[xi].sup.(k)]).
```

3.3. Application of the restarted RRE method to the solution of NARE. We apply a restarted RRE algorithm to the vectors computed by the iteration scheme (2.20) to solve the NARE (2.3)-(2.5). Two approaches for restarted RRE are detailed in Algorithms 2 and 3. Algorithm 2 applies a restarted RRE method directly to the vector sequence [{[w.sup.(k)]}.sub.k[member of]N], where [w.sup.(k)] = [[[([u.sup.(k)]).sup.T], [([v.sup.(k)]).sup.T]].sup.T], while Algorithm 3 applies a restarted RRE method to the vector sequences [{[u.sup.(k)]}.sub.k[member of]N] and [{[v.sup.(k)]}.sub.k[member of]N] separately. Both approaches accelerate the convergence of the vector sequences but Figure 3.1 shows that the technique of Algorithm 2 works better, especially near the critical case. For later numerical experiments, we therefore use Algorithm 2.
```
Algorithm 2 The restarted RRE(r) applied to
[{[w.sup.(k)]}.sub.k], r is fixed.

Input: k = 0, [u.sup.(0)] = [v.sup.(0)] = 0, choose an integer r.
Then, [w.sup.(0)] = [[[([u.sup.(0)]).sup.T], [([v.sup.(0)])
.sup.T]].sup.T] = 0.
For k = 1, 2, ...,
[y.sup.(0)] = [u.sup.(k-1)]; [z.sup.(0)] = [v.sup.(k-1)];
[s.sup.(0)] = [([y.sup.(0)], [z.sup.(0)]).sup.T];
[y.sup.(j+1)] = [[PHI].sub.P] ([z.sup.(j)]), [z.sup.(j+1)]
= [[PHI].sub.Q] ([y.sup.(j+1)]), [s.sup.(j+1)] =
[([y.sup.(j+1)], [z.sup.(j+1)]).sup.T], j = 0, ..., r - 1.
Compute the approximations [t.sup.(r-1)] by applying the RRE
Algorithm 1 on the vectors ([s.sup.(0)],[s.sup.(1)], ...
,[s.sup.(r)]);
If t(r-1) satisfies accuracy test, stop.
Else
[w.sup.(k)] = [([u.sup.(k)],[v.sup.(k)]).sup.T] = [t.sup.(r-1)];
End

Algorithm 3 The restarted RRE(r) applied to [{[u.sup.(k)]}.sub.k] and
[{[v.sup.(k)]}.sub.k], r is fixed.

Input: k = 0, [u.sup.(0)] = [v.sup.(0)] = 0, choose an integer r.
For k = 1, 2, ...,
[y.sup.(0)] = [u.sup.(k-1)]; [z.sup.(0)] = [v.sup.(k-1)];
[y.sup.(j+1)] = [[PHI].sup.P]([z.up.(j)], [z.sup.(j+1)] =
[[PHI].sub.Q]([y.sup.(j+1])], j= 0, ...,r - 1.
Compute the approximations [t.sup.(r-1).sub.1]) and
[t.sup.(r-1).sub.2]) by applying the RRE of Algorithm 1 on
the vectors [[([y.sup.(0)], [z.sup.(0)]).sup.T],
[([y.sup.(1)], [z.sup.(1)]).sup.T], ...,[([y.sup.(r)]
[z.sup.(r)]).sup.T]];
If [t.sup.r-1.sub.1]) and [t.sup.(r-1).sub.2] satisfy accuracy
test, stop.
Else
[u.sup.(k)] = [t.sup.r-1].sub.1]); ,[v.sup.(k)] =
[t.sup.r-1].sub.2]);
End
```

3.4. The choice of r. The RRE algorithm should be restarted every r iterations, for some integer r > 1, to avoid the increase in computational work and storage as k increases.

If w* = [PHI] (w*) is a fixed point and J([PHI], w*) is the Jacobian matrix of [PHI] at w*, then

[w.sup.(k+1)] - w* = J([PHI], w*)([w.sup.(k)] - w*) + O[([parallel]([w.sup.(k)] - w*)[parallel].sup.2]).

It suffices to examine the eigenvalues of the Jacobian matrix Jw*). Going back to the iteration (2.21) and observing the eigenvalues of the Jacobian matrix (2.25) at the solution, it can be seen that these eigenvalues range between zero and one. Most of these eigenvalues are close or equal to zero, except for a few which are close to one. Therefore, choosing a small integer r is sufficient. In particular, r = 4 yields good results. Figure 3.2 shows the distribution of the spectrum of the Jacobian matrix at the solution for n = 512, and ([alpha], c) = (0.001,0.999) with spectral radius [rho](J([PHI], ([u.sup.*], [v.sup.*]))) [approximately equal to] 0.86.

4. Numerical experiments. A numerical example is presented in this section to illustrate the performance of the new approach for solving the vector equation (2.17). We consider a special kind of Riccati equation (2.3)-(2.5). The constants c and w are given by a numerical quadrature formula on the interval [0,1], which is obtained by dividing [0,1] into n/4 subintervals of equal length and applying a composite Gauss-Legendre quadrature rule with 4 nodes in each subinterval.

Computations are performed for different choices of the parameters ([alpha], c) and for different values of n using MATLAB 7.4 (R2007a) on an Intel-I5 quad-core 2.5 Ghz processor (4 GB RAM) with Fedora 18 (Linux) and with approximately 15 significant decimal digits.

The stopping criterion is given by

ERR = [[parallel][w.sup.(k+1) - [w.sup.(k)][parallel]/[parallel][w.sup.(k+1)[parallel]] [less than or equal to] tol

for tol = [10.sup.-10].

Table 4.1 compares the three vector extrapolation methods RRE,MPE, and MMPE for two values of r. One can observe that for r = 4, the three methods are comparable and give close results, while for a larger r (r = 10), the RRE method is the best.

Table 4.2 compares the iterative method proposed by Lin (2.12), its modified version (2.13), and the application of RRE. Denote by ILin the iterative method proposed by Lin  with the choice of [T.sub.k] given in (2.18), by MILin in its modified version, and by RRE the modified version with the application of restarted reduced rank extrapolation every 4 iterations of Algorithm 2. Table 4.2 shows how the three methods converge to the minimal positive solution of (2.20) for several a and c values. The RRE method is seen to outperform ILin and MILin. Figure 4.1 shows the performance of RRE in comparison with ILin and MILin for ([alpha], c) = (0.001, 0.999).

We now compare RRE and the fast Newton method proposed in . Computations of this table have been implemented in Fortran 90 on Linux with 2.5 Ghz and with about 15 significant decimal digits and for tol = 10-10. The Fortran code implemented in  is used. Denote by LuF the fast Newton method, which is based on a fast LU algorithm that reduces the cost of Newton's algorithm proposed by Lu  from O([n.sup.3]) to O([n.sup.2]); see Section 2.1.3. Table 4.3 compares the restarted RRE method with LuF in terms of CPU time (in seconds) for different n and for ([alpha], c) = ([10.sup.-8], 1 - [10.sup.-6]). It can be seen that RRE is faster than LuF also when the convergence is slow for ([alpha], c) close to (0,1) and for large n.

5. Conclusions. In this paper, we dealt with a special kind of Riccati equations, NARE, that arises in transport theory. We presented some efficient iterative methods which solve this kind of equations either by computing the minimal positive solution of a vector equation or directly to NARE. We applied a polynomial vector extrapolation method, namely reduced rank extrapolation (RRE), to the iterative method proposed by Lin  which computes the minimal positive solution of NARE as the minimal positive solution of a vector equation. Numerical experiments showed the advantage of using the RRE method, especially near the critical case when ([alpha], c) is close to (0,1) and for large n.

Acknowledgments. We would like to thank F. Poloni for providing us with his program for the LuF algorithm. We also thank the referees for their useful comments and valuable suggestions.

REFERENCES

 Z. Z. BAI, Y. H. GAO, AND L. Z. LU, Fast iterative schemes for nonsymmetric algebraic Riccati equations arising from transport theory, SIAM J. Sci. Comput., 30 (2008), pp. 804-818.

 L. BAO, Y. LIN, AND Y. WEI, A modified simple iterative method for nonsymmetric algebraic Riccati equations arisingin transport theory, Appl. Math. Comput., 181 (2006), pp. 1499-1504.

 R. BELLMAN AND G. M. WING, An Introduction to Invariant Embedding, John Wiley, New York, 1975.

 D. A. BINI, B. IANNAZZO, AND F. POLONI, A fast Newton's method for a nonsymmetric algebraic Riccati equation, SIAM J. Matrix Anal. Appl., 30 (2008), pp. 276-290.

 D. A. BINI, B. MEINI, AND B. IANNAZZO, Numerical solution of algebraic Riccati equations. Fundamentals ofAlgorithms, SIAM, Philadelphia, 2012.

 C. BREZINSKI, Generalisations de la transformation de Shanks, de la table de Pade et de l' e-algorithme, Calcolo, 12 (1975), pp. 317-360.

 C. BREZINSKI AND M. REDIVO ZAGLIA, Extrapolation Methods. Theory and Practice, North-Holland, Amsterdam, 1991.

 S. CABAY AND L.W. JACKSON, A polynomial extrapolation method for finding limits and antilimits of vector sequences, SIAM J. Numer. Anal., 13 (1976), pp. 734-752.

 F. CORON, Computation of the asymptotic states for linear half space kinetic problem, Transport Theory Statist. Phys., 19 (1990), pp. 89-114.

 R. P. EDDY, Extrapolating to the limit of a vector sequence., in Information Linkage between Applied Mathematics and Industry, P. C. C. Wang, ed., Academic Press, New York, 1979, pp. 387-396.

 W. F. FORD AND A. SIDI, Recursive algorithms for vector extrapolation methods, Appl. Numer. Math., 4 (1988), pp. 477-489.

 B. D. GANAPOL, An investigation of a simple transport model, Transport Theory Statist. Phys., 21 (1992), pp. 1-37.

 I. GOHBERG, T. KAILATH, AND V. OLSHEVSKY, Fast Gaussian elimination with partial pivoting for matrices with displacement structure, Math. Comp., 64 (1995), pp. 1557-1576.

 C. H. GUO AND A. J. LAUB, On the iterative solution of a class of nonsymmetric algebraic Riccati equations, SIAM J. Matrix Anal. Appl., 22 (2000), pp. 376-391.

 C. H. GUO AND W. W. LIN, Convergence rates of some iterative methods for nonsymmetric algebraic Riccati equations arising in transport theory, Linear Algebra Appl., 432 (2010), pp. 283-291.

 K. JBILOU AND H. SADOK, LU implementation of the modified minimal polynomial extrapolation method for solving linear and nonlinear systems, IMA J. Numer. Anal., 19 (4) (1999), pp. 549-561.

--, Vector extrapolation methods. Applications and numerical comparison, J. Comput. Appl. Math., 122 (2000), pp. 149-165.

 J. Juang, Existence of algebraic matrix Riccati equations arising in transport theory, Linear Algebra Appl., 230 (1995), pp. 89-100.

 J. JUANG, C. L. HSING, AND P. NELSON, Global existence, asymptotics and uniqueness for the reflection kernel of the angularly shifted transport equation, Math. Models Methods Appl. Sci., 5 (1995), pp. 239251.

 J. JUANG AND I-DER CHEN, Iterative solution for a certain class of algebraic matrix Riccati equations arising in transport theory, Transport Theory Statist. Phys., 22 (1993), pp. 65-80.

 J. JUANG AND W. W. LIN, Nonsymmetric algebraic Riccati equations and Hamiltonian-like matrices, SIAM J. Matrix Anal. Appl., 20 (1999), pp. 228-243.

 C. KENNEY AND A. J. LAUB, Rational iterative methods for the matrix sign function, SIAM J. Matrix Anal. Appl., 12 (1991), pp. 273-291.

 P. LANCASTER AND L. RODMAN, Algebraic Riccati equations, Oxford University Press, New York, 1995.

 Y. LIN, A class of iterative methods for solving nonsymmetric algebraic Riccati equations arising in transport theory, Comput. Math. Appl., 56 (2008), pp. 3046-3051.

 L. Z. LU, Solution form and simple iteration ofa nonsymmetric algebraic Riccati equation arising in transport theory, SIAM J. Matrix Anal. Appl., 26 (2005), pp. 679-685.

--, Newton iterations for a non-symmetric algebraic Riccati equation, Numer. Linear Algebra Appl., 12 (2005), pp. 191-200.

 V. L. MEHRMANN, The Autonomous Linear Quadratic Control Problem, Lecture Notes in Control and Information Sciences 163, Springer, Berlin, 1991.

 M. MESINA, Convergence acceleration for the iterative solution of the equations X = AX + f, Comput. Methods Appl. Mech. Engrg., 10 (1977), pp. 165-173.

 B. P. PUGATCHEV, Acceleration of the convergence of iterative processes and a method for solving systems of nonlinear equations, U.S.S.R. Comput. Math. and Math. Phys., 17 (1977), pp. 199-207.

 A. SIDI, Convergence and stability properties of minimal polynomial and reduced rank extrapolation algorithms, SIAM J. Numer. Anal., 23 (1986), pp. 197-209.

 --, Efficient implementation of minimal polynomial and reduced rank extrapolation methods, J. Comput. Appl. Math., 36 (1991), pp. 305-337.

 --, Extrapolation vs. projection methods for linear systems of equations, J. Comput. Appl. Math., 22 (1988), pp. 71-88.

 A. SIDI, W. F. FORD, AND D. A. SMITH, Acceleration of convergence of vector sequences, SIAM J. Numer. Anal., 23 (1986), pp. 178-196.

 D. A. SMITH, W. F. FORD, AND A. SIDI, Extrapolation methods for vector sequences, SIAM Rev., 29 (1987), pp. 199-233.

 P. WYNN, Acceleration techniques for iterated vector and matrix problems, Math. Comp., 16 (1962), pp. 301322.

 --, On a device for computing the [e.sub.m]([S.sub.n]) transformation, Math. Tables Aids Comput., 10 (1956), pp. 91-96.

ROLA EL-MOALLEM ([dagger]) AND HASSANE SADOK ([double dagger])

* Received July 25, 2013. Accepted for publication October 1, 2013. Published online December 28, 2013 Recommended by L. Reichel.

([dagger]) Laboratoire Paul Painleve, Universite Lille 1--Sciences et Technologies, Batiment M2, 59 655 Villeneuve d'Ascq Cedex, France (rola. el-moallem@math.univ-lille1.fr).

([double dagger]) Laboratoire de Mathematiques Pure et Appliquees, Universite du Littoral, Centre Universitain le la Mi-Voix, Batiment H. Poincarre, 50 Rue F. Buisson, BP 699, 62228 Calais Cedex, Francs sadok@lmpa.univ-littoral.fr).
```
Table 4.1

Comparison of the three vector extrapolation methods
for different r and for n = 2048.

r=4                            RRE               MPE

CPU time (in seconds)          3.8               3.8
Residual norm           3.5.[10.sup.-11]    2.[10.sup.-11]

r = 10                         RRE               MPE

CPU time (in seconds)         6.59               7.57
Residual norm           6.87.[10.sup.-12]   9.[0.sup.-11]

r=4                           MMPE

CPU time (in seconds)          3.9
Residual norm           3.9.[10.sup.-11]

r = 10                        MMPE

CPU time (in seconds)         10.36
Residual norm           1.0.[10.sup.-11]

Table 4.2
Numerical results for n=256 with different ([alpha],c).

Iteration   Residual    CPU
[alpha]       c                 Method     steps       norm     time

ILin      4732      9.98e-11   5.47
[10.sup.-8]   1 - [10.sup.-6]   MILin      2517      9.98e-11   2.87
RRE        20       5.03e-14   0.11

ILin      1813      9.97e-11   2.08
[10.sup.-5]   1 - [10.sup.-5]   MILin       955      9.93e-11   1.08
RRE         7       3.83e-14   0.05

ILin       674      9.98e-11   0.77
0.0001        0.9999            MILin       353      9.88e-11   0.4
RRE         7       9.8e-16    0.04

ILin       246      9.7e-11    0.3
0.001         0.999             MILin       129      9.04e-11   0.14
RRE         9       2.99e-14   0.05

ILin       12       4.01e-11   0.01
0.5           0.5               MILin        7       3.32e-11   0.008
RRE         3       9.72e-17   0.02

Table 4.3
Comparison in terms of CPU time in seconds for different n.

[alpha] = [10.sup.-8], c = 1 - [10.sup.-6]

n       LuF    RRE(4)

512     0.35   0.24
1024    1.37   0.96
2048    6.6    3.8
```
COPYRIGHT 2013 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.