1 0BIntroduction

Lanczos-type algorithms are effective methods to solve systems of linear equations (SLE). One of their attractions is that, in exact arithmetic, they achieve an exact solution in n steps where n is the dimension of the problem [1]. However, they have a weakness linked to their need to estimate the coefficients of the orthogonal polynomials used in the recurrence relationships on which they are based. For instance, when the denominator of any of these coefficients is zero, then the corresponding orthogonal polynomial does not exist; the solution process stops and fails to reach the solution. This is known as the breakdown phenomenon. There are many forms of this, such as hard, soft, curable, incurable [2], true, and ghost breakdowns [1,3]. Here, we are concerned with true or hard breakdown.

Note that the more iterations we do, the more likely breakdown will occur. Since we do not work in exact arithmetic, to get good approximate solutions, a number of iterations exceeding the dimension of the problem is often required. Therefore, breakdown is even more likely to occur. That is why Lanczos-type algorithms suffer from breakdown and that is why approaches to dealing with this problem are continually being developed [4-6].

There are already many strategies to tackle breakdown, such as the look-ahead strategy [7], also called the method of recursive zoom (MRZ), and the lookaround strategy [8,9]. These typically try to get over and/or around the nonexisting orthogonal polynomials. Attempts based on different strategies, such as switching between Lanczos-type algorithms and restarting them, have also been considered [4-6]. These latter approaches have been shown to perform better than MRZ in terms of robustness [5,6].

In this paper, we will focus on restarting from three specific points which are attractive as starting points: the last iterate preceding breakdown, the iterate with the lowest residual norm, and the iterate whose entries are the median values of all entries of the iterates generated so far.

It is reasonable to restart from the point immediately before the breakdown occurred since this leads one to believe that the sequence that has been broken is continued somehow, albeit in a different Krylov space. However, there is no theoretical evidence for this. Indeed, it may be counter-productive since just before breakdown the quality of the iterates deteriorated from the point of view of the residual norm. Our numerical results support this view, as will be seen later.

Restarting from the iterate with minimum residual norm in the sequence of solutions generated so far is also reasonable since we now have a good starting point which is not just a guess. Numerical results included here, also support this idea.

Finally, restarting from an approximate solution whose entries are the median of the entries of all the iterates generated so far also sounds like a reasonable starting point. This is because one might think that the median values of the entries capture the intrinsic qualities of the sequence of solutions found so far. However, it is also clear that it will include information from points that are very bad since they lead to the process of blowing up (breakdown) in terms of their residual values. As one would expect, the results point to this artificially made-up iterate as being in between the other two in terms of quality.

We will also look at the theoretical aspect of restarting. This has been considered previously in [2,10]. Numerical results obtained on problems ranging from 1000 to 70000 dimensions will be provided.

2 Lanczos-type Algorithms for SLE's and the Issue of Breakdown

2.1 Lanczos-type Algorithms: Review

Consider the system of linear equations

\[A\mathbf{x} = \mathbf{b},\tag{1}\] where ∈ and vectors and ∈ . Lanczos-type methods for solving Eq. (1) can be derived by combining the method of Krylov subspace and the method of projections, [11]. It is started by choosing an initial approximate solution 0 and defining the sequence of vectors by two conditions:

\[\mathbf{x}_k - \mathbf{x}_0 \in K_k(A, \mathbf{r}_0), \tag{2}\] and

\[\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k \perp L_k = K_k(A, \mathbf{y}), \qquad (3)\] where 0 = – 0 is the corresponding residual vector, is an arbitrary nonzero vector, and (, 0) is a Krylov subspace. From Eq. (2) we have:

\[\mathbf{x}_k - \mathbf{x}_0 = -\alpha_1 \mathbf{r}_0 - \alpha_2 A \mathbf{r}_0 - \dots - \alpha_k A^{(k-1)} \mathbf{r}_0. \tag{4}\]

Thus, multiplying both sides by , adding and substracting , we obtain:

\[A(\mathbf{x}_{k} - \mathbf{x}_{0}) = A \left(-\alpha_{1}\mathbf{r}_{0} - \alpha_{2}A\mathbf{r}_{0} - \cdots - \alpha_{k}A^{(k-1)}\right)\mathbf{r}_{0},\] \[A\mathbf{x}_{k} - A\mathbf{x}_{0} = -\alpha_{1}A\mathbf{r}_{0} - \alpha_{2}A^{2}\mathbf{r}_{0} - \cdots - \alpha_{k}A^{k}\mathbf{r}_{0},\] \[\mathbf{b} - \left(A\mathbf{x}_{k} - A\mathbf{x}_{0}\right) = \mathbf{b} + \alpha_{1}A\mathbf{r}_{0} + \alpha_{2}A^{2}\mathbf{r}_{0} + \cdots + \alpha_{k}A^{k}\mathbf{r}_{0},\] \[\mathbf{b} - A\mathbf{x}_{k} = (b - A\mathbf{x}_{0}) + \alpha_{1}A\mathbf{r}_{0} + \alpha_{2}A^{2}\mathbf{r}_{0} + \cdots + \alpha_{k}A^{k}\mathbf{r}_{0},\] \[\mathbf{r}_{k} = \mathbf{r}_{0} + \alpha_{1}A\mathbf{r}_{0} + \alpha_{2}A^{2}\mathbf{r}_{0} + \cdots + \alpha_{k}A^{k}\mathbf{r}_{0}.\] (5)

From the last relation, we can write

\[\mathbf{r}_k = P_k(A)\mathbf{r}_0,\tag{6}\] where () = 1 + 1 + ⋯ + is of degree at most, and satisfies the normality condition (0) = 1.

Following the definition of the family of orthogonal polynomials given in [3], we set \(c_i = \langle \mathbf{y}, A^i \mathbf{r}_0 \rangle\), with c a linear functional on the vector space of polynomials defined by

\[c_i = c(x^i), \quad i = 1, 2, \cdots \tag{7}\]

The condition Eq. (6) can be written as

\[c(t^{i}P_{k}(t)) = 0, i = 0,1,...,k-1,\] (8)

which shows that \(P_k\) is the polynomial of degree k at most, normalized by \(P_k(0) = 1\), belonging to the family of orthogonal polynomials with respect to functional c. Thus, the Lanczos method consists of determining implicitly the polynomial \(P_k\), computing recursively the residual \(\mathbf{r}_k = P_k(A)\mathbf{r}_0\), and finding \(\mathbf{x}_k\) from \(\mathbf{x}_k = \mathbf{b} - A\mathbf{r}_k\) without inverting the matrix A.

Some Lanczos-type algorithms have been discussed by Brezinski and his team, which are expressed in a table containing formulae \(A_i\) and \(B_j\), [12-14]. Basically, the computations of the recurrence relations are based on \(P_k\) and \(P_k^{(1)}\) respectively, with degrees of the polynomials on the right and left hand sides differing by one or two degrees at most. The new formulae which involve the polynomials with a difference in degrees of at most two or three can be seen in [4,15].

2.2 The Breakdown Phenomenon

Breakdown in Lanczos-type algorithms can be caused by either the non-existence of some orthogonal polynomials or because the recurrence relationship used is not appropriate. Following the explanation in [16], consider how breakdown occurs when estimating the coefficients of orthogonal polynomials.

Consider the three term recurrence relationship of the monic polynomial \(P_{k+1}\) as follows:

\[P_{k+1}(t) = (\alpha t + \beta) P_k(t) + \gamma P_{k-1}(t), \tag{9}\] where \(k=0,1,\cdots\). This polynomial satisfies the conditions: (1) \(P_{-1}(t)=0\), (2) \(P_0(t)=1\), and (3) \(c(t^iP_k)=0\), for \(i=0,1,\cdots,k-1\), where c is a linear function. The last condition, which is also known as the orthogonality condition, gives:

\[c(t^{i} P_{-}k + 1) = \alpha c(t^{i} + 1 P_{-}k(t)) + \beta c(t^{i} P_{-}k(t)) + \gamma c(t^{i} P_{-}(k - 1)(t)),\] (10)

For \(i = 0,1,\dots,k\). In particular, for i = k-1 and i = k we get relations

\[0 = \alpha c(t^k P_k(t)) + \gamma c(t^k - 1 \{P_k - 1(t)\})\] (11)

\[0 = \alpha c \left( t^k + 1 P_k(t) \right) + \beta c \left( t^k P_k(t) \right) + \gamma c \left( t^k P_{k-1}(t) \right)\] (12)

respectively. The second condition leads to the equation,

\[1 = \beta + \gamma. \tag{13}\]

As we can see, we now have a system of linear equations in \(\alpha\), \(\beta\), \(\gamma\). Solving for \(\alpha\), \(\beta\), \(\gamma\) gives:

\[\alpha = 1\] \[\gamma = \frac{-c(t^{k}P_{k}(t))}{c(t^{k-1}P_{k-1}(t))},\] \[\beta = \frac{\gamma c(t^{k}P_{k-1}(t)) - c t^{k+1} P_{k}(t)}{c(t^{k}P_{k}(t))}.\] (14)

In order for relation Eq. (9) to be finite, the scalars \(c\left(t^{k-1}P_{k-1}\left(t\right)\right)\) and \(c\left(t^{k}P_{k}\left(t\right)\right)\) must not be zero, or sufficiently close to zero, i.e. small enough to cause an NaN output in the computation. Otherwise, a breakdown occurs.

3 Restarting Lanczos-type Algorithms

Here, we introduce the different implementations of restarting Lanczos-type algorithms from three different points. RLLastIt restarts from the last iterate; RLMinRes restarts from the iterate with the minimum residual norm; RLMedVal restarts from a point made up of the median values of all the entries of previous iterates in a sequence generated by the Lanczos-type algorithm under consideration. Restarting in all cases is of the pre-emptive type as described in [5].

3.1 Restarting Lanczos-type Algorithms from the Last Iterate

Assume that we stop the algorithm after k iterations. Also assume that breakdown does not occur within k iterations. Note that a high k means the algorithm may break down. Consider the iterate \(\mathbf{x}_k\), as the last iterate, with residual norm \(\|\mathbf{r}_k\|\). Initialize the algorithm as follows,

\[\mathbf{x}_0 = \mathbf{x}_k \tag{15}\]

\[\mathbf{y} = \mathbf{b} - A\mathbf{x}_0 \tag{16}\]

RLLastIt is described in Algorithm 1.

Algorithm 1 The RLastIt Algorithm.

• 1. Fix the number of iterations to, say, k, and the tolerance, \(\epsilon\), to 1E-13.
• 2. Run a Lanczos-type algorithm for k iterations and get the approximate solution \(sol_{last} = \mathbf{x}_k\), as well as the residual norm \(norm_{last} = ||\mathbf{r}_k||\).
• 3. while \(norm_{last} \ge \epsilon\) do
• 4. Initialize the algorithm with

\[\mathbf{x} = sol_{last},\]
\(\mathbf{y} = \mathbf{b} - A \mathbf{x}.\)

• 5. Run the Lanczos-type algorithm for k iterations.
• 6. endwhile.
• 7. Take \(sol_{last}\) as the approximate solution.
• 8. Stop.

3.2 Restarting Lanczos-type Algorithms from the Iterate with the Minimum Residual Norm

Because of the potential accumulation of errors, the last iterate after k iterations is not necessarily the one with the smallest residual norm. It is therefore reasonable to consider restarting from the point among the k iterates with the lowest residual norm. Figure 1 illustrates the residual norm behaviour of the iterates after running the Orthodir algorithm [12] when solving problems in 300 dimensions using 150 iterations. This figure describes the situation mentioned earlier that the last iterate is not always the best one. In this case, the iterate with the lowest residual norm is the \(29^{th}\).

To run RLMinRes, we collect all of the iterates generated in k iterations and their residual norms.

\[data_{sol} = \{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_k\},\tag{17}\]

\[data_{resnorm} = \{ \|\mathbf{r}_1\|, \|\mathbf{r}_2\|, \cdots, \|\mathbf{r}_k\| \}.\] (18)

We then find the minimum value of Eq. (18) as

\[norm_{\min} = \min(data_{resnorm}), \tag{19}\] and its index m. Our new solution is the iterate in Eq. (17) with index m.

Figure 1 The behaviour of the residual norms of an SLE in 300 dimensions.

As in RLLastIt, RLMinRes is initialized as follows

\[\mathbf{x}_0 = \mathbf{x}_m,\tag{20}\] and y as in Eq. (16). If we repeat this step, we will find a good approximate solution as explained previously. The procedures of LMinRes and RLMinRes are presented in Algorithm 2 and Algorithm 3 respectively.

Algorithm 2 The LMinRes Algorithm.

• 1. Fix the number of iterations to, say, k, and the tolerance, \(\epsilon\), to 1E-13.
• 2. Collect all k iterates as in Eq. (17).
• 3. Collect all of the residual norms as in Eq. (18).
• 4. Compute the minimum values as in Eq. (19) and specify the index of the minimum value as m.
• 5. Obtain the approximate solution as well as the residual norm as follows

\[sol_{min} = data_{sol}(m), (21)\]

\[norm_{min} = \|\mathbf{r}_m\|. \tag{22}\]

6. Stop.

Algorithm 3 TheRLMinRes algorithm.

• 1. Initialize a Lanczos-type algorithm with initial guesses \(\mathbf{x}_0\) and \(\mathbf{y}\).
• 2. Run LMinRes algorithm for k iterations and obtain \(sol_{min}\) and \(norm_{min}\).
• 3. while \(norm_{min} \geq \epsilon\) do
  • i. \(\mathbf{x} = sol_{min}\),
  • ii. \(\mathbf{y} = \mathbf{b} A \mathbf{x}\).
• 4. Run LMinRes algorithm for k iterations.
• 5. endwhile.
• 6. Take \(sol_{min}\) as the approximate solution.
• 7. Stop.

3.3 Restarting Lanczos-type Algorithms with the Vector of Median Values

Consider again relations Eq. (17) and Eq. (18). For each \(x_i\), i=1,2,...,n, we set

\[w_{1} = \left\{ x_{1}^{(1)}, x_{2}^{(1)}, \cdots, x_{k}^{(1)} \right\}\] \[w_{2} = \left\{ x_{1}^{(2)}, x_{2}^{(2)}, \cdots, x_{k}^{(2)} \right\}\] \[\vdots\] \[w_{n} = \left\{ x_{1}^{(n)}, x_{2}^{(n)}, \cdots, x_{k}^{(n)} \right\}\] (23)

where \(w_i\) is arranged such that it consists of the \(i^{th}\) entry of each vector \(\mathbf{x}_i\), for i=1,2,...,n. For instance, \(w_1\) is the set of all of the first entries of the iterates, \(w_2\) is the set of all of the second entries of the iterates, etc. Thus, if we calculate the mean, median, and mode values of each \(w_i\) and we use them as the first entry, the second entry, etc... of a vector, it leads to a new vector solution. In other words, the new vector contains the values of either mean, median, or mode of all the first entries, the second entries, etc... of all the iterates. Our investigation of these entries showed that the best results are obtained with the median values. The approximate solution based on the median values of entries is as follows,

\[\mathbf{x}_{med} = \begin{bmatrix} x_1^{med} \\ x_2^{med} \\ \vdots \\ x_n^{med} \end{bmatrix}\] (24)

and \(x_i^{med} = med \{w_i\}\), for i = 1, 2, ..., n.

The procedure described above for creating an approximate solution from the median values of the entries of iterates generated by the Lanczos-type algorithm used is presented in algorithmic form as Algorithm 4 or LMedVal.

Algorithm 4 The LMedVal algorithm.

• 1. Initialize a Lanczos-type algorithm with initial guesses \(\mathbf{x}_0\) and \(\mathbf{y}\). Run it for k iterations
• 2. Collect all of k iterates as in Eq. (17).
• 3. for i = 1, 2, ..., n
• 4. Set \(w_i\) as in Eq. (23)
• 5. Compute the median value of the entries of each \(w_i\) as follows \(med(i) = median(w_i)\).
• 6. endfor
• 7. Compute the approximate solution as follows

\[\mathbf{x}_{med} = med(1:n)^T\]

8. Compute the residual norm of the solution as follows

\[res_{med} = \mathbf{b} - A \mathbf{x}_{med},\]
\(norm\_med = norm(res_{med}).\)

9. Stop

The procedure described above to restart the Lanczos-type algorithm from the output of LMedVal is presented in algorithmic form as Algorithm 5 or RLMedVal. Its initialization is as follows,

\[\mathbf{x}_0 = \mathbf{x}_{med},\tag{25}\] and vector \(\mathbf{y}\) as in Eq. (16), where \(\mathbf{x}_{med}\) as in Eq. (24).

Algorithm 5 TheRLMedVal algorithm

• 1. Fix the number of iterations to k and the tolerance \(\epsilon\) to 1E-13.
• 2. Run LMedVal for k iterations and obtain \(sol_{med}\) as well as \(norm_{med}\).
• 3. while \(norm_{med} \ge \epsilon\) do
  • i. \(\mathbf{x} = \mathbf{x}_{med}\),
  • ii. \(\mathbf{y} = \mathbf{b} A\mathbf{x}\).
• 4. Run LMedVal for k iterations.
• 5. endwhile.
• 6. Take \(sol_{med}\) as the approximate solution.
• 7. Stop.

3.4 The Theoretical Point of View

Suppose we solve an SLE using RLLastIt. Denote \(x_k^{(j)}\), j=1,2,...,s, and s>0 integer, as the \(k^{th}\) iterate of the \(j^{th}\) run of a Lanczos-type algorithm. In other words, s is the number of cycles of RLLastIt. Assume that the residual norm of the first cycle is greater than the convergence tolerance, i.e.

\[\left\|\mathbf{r}_{k}^{(1)}\right\| > \epsilon_{1},\tag{26}\] for some \(\epsilon_1 > 0\). Obviously, if this is not the case, the algorithm would have stopped with an acceptable approximate solution. Consider the results of Section 3.1 and by the definition of the Krylov subspace given in Section 2.1, iterate

\(\mathbf{x}_k^{(2)}\) of the second cycle of RLLasIt satisfies the following two conditions:

\[\mathbf{x}_{k}^{(2)} - \mathbf{x}_{0} \in K_{k}(A, \mathbf{r}_{0}),\] (27)

\[\mathbf{r}_k^{(2)} = \mathbf{b} - A \mathbf{x}_k^{(2)} \perp K_k(A^T, \mathbf{y}), \tag{28}\] where \({\bf r}_0 = {\bf r}_k^{(1)}\).

Substitute Eq. (15) into Eq. (27) to get:

\[\mathbf{x}_{k}^{(2)} - \mathbf{x}_{k}^{(1)} \in K_{k}\left(A, \mathbf{r}_{k}^{(1)}\right).\] (29)

Following the results in Eq. (5), we have:

\[\mathbf{r}_{k}^{(2)} = \mathbf{r}_{k}^{(1)} + \alpha_{1} A \mathbf{r}_{k}^{(1)} + \alpha_{2} A^{2} \mathbf{r}_{k}^{(1)} + \dots + \alpha_{k} A^{k} \mathbf{r}_{k}^{(1)}\] \[= P_{k}(A) \mathbf{r}_{k}^{(1)}, \tag{30}\] where \(P_k(A) = 1 + \alpha_1 A + \cdots + \alpha_k A^k\).

Calculating the norm of \(\mathbf{r}_k^{(2)}\) in Eq. (30) yields

\[\left\| \mathbf{r}_{k}^{(2)} \right\| = \left\| P_{k}(A)\mathbf{r}_{k}^{(1)} \right\|\] \[\leq \left\| P_{k}(A) \right\| \left\| \mathbf{r}_{k}^{(1)} \right\|. \tag{31}\]

If this residual norm is still bigger than the tolerance \(\epsilon_1\), we restart again, in which case we will have the third cycle, giving

\[\mathbf{r}_k^{(3)} = P_k(A)\mathbf{r}_k^{(2)}.\] (32)

Calculating the norm of \(\mathbf{r}_k^{(3)}\) as

\[\|\mathbf{r}_{k}^{(3)}\| = \|P_{k}(A)\mathbf{r}_{k}^{(2)}\|\] \[\leq \|P_{k}(A)\| \|\mathbf{r}_{k}^{(2)}\| \leq \|P_{k}(A)\|^{2} \|\mathbf{r}_{k}^{(1)}\|, \text{ from Eq. (31) (33)}\]

If this procedure is continued, then for the \(s^{th}\) cycle we have:

\[\|\mathbf{r}_{k}^{(s)}\| = \|P_{k}(A)\mathbf{r}_{k}^{(s-1)}\| \le \|P_{k}(A)\|^{s} \|\mathbf{r}_{k}^{(1)}\|\] (34)

From Eq. (31), Eq. (33), and Eq. (34), and since \(||P_k(A)|| > 0\), we conclude that

\[\left\|\mathbf{r}_{k}^{(s)}\right\| \leq \left\|\mathbf{r}_{k}^{(s-1)}\right\| \leq \cdots \left\|\mathbf{r}_{k}^{(1)}\right\| \tag{35}\] where \(s \ge 1\) is the number of cycles. This leads to the following result:

Theorem 1 Suppose we solve an SLE using a Lanczos-type algorithm. Let \(x_k^{(1)}\) be the last iterate generated by this algorithm after k iterations. Given a tolerance \(\epsilon > 0\), we assume that the associated residual norm \(\|r_k^{(1)}\| > \epsilon\). Restarting the algorithm with \(x_k^{(1)}\) allows it to generate an iterate with a better residual norm than those of the previous iterates.

Practically, breakdown may still occur between cycles, which causes the algorithm to stop before reaching a good approximate solution. We cannot guarantee that a breakdown will not occur in the first cycle, or any cycle for that matter.

4 Numerical Results and Discussion

We solved different size problems ranging from dimension 1000 to 70000. The test problems were carried out in MatLab 2012b using two different systems, Windows with processor RAM 6GB for solving the medium-scale problems, and Unix0 with 256 GB for solving the large-scale problems. We implemented the RLLastIt, RLMinRes, and RLMedVal algorithms, as well as other supporting procedures such as LMedVal. The experiments were carried out via Algorithm 6, which calls RLLastIt, RLMinRes, RLMedVal.

Algorithm 6 Restarting Lanczos-type algorithms from three different points.

• 1. Fix the number of iterations to k and the tolerance \(\epsilon\) to 1E-13.
• 2. Run RLLastIt, RLMinRes, RLMedVal for k iterations.
• 3. WHILE \(norm_{last} \ge \epsilon\) do
• 4. Run RLLastIt for k iterations.

• 5. Run RLMinRes for k iterations.
• 6. Run RLMedVal for k iterations.
• 7. ENDWHILE.
• 8. Take \(l_{last}\), \(sol_{min}\) and \(sol_{med}\) as the approximate solution of each above algorithm, respectively.
• 9. Stop.

We solve several systems of linear equations \(A\mathbf{x} = \mathbf{b}\), with matrix A being of the form

\[A = \begin{pmatrix} B & -I & \cdots & \cdots & 0 \\ -I & B & -I & \cdots & \cdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ & & -I & B & -I \\ 0 & \cdots & \cdots & -I & B \end{pmatrix}\] \[(36)\] where

\[B = \begin{pmatrix} 4 & \alpha & \cdots & \cdots & 0 \\ \beta & 4 & \alpha & \cdots & \cdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ & \beta & 4 & \alpha \\ 0 & \cdots & \cdots & \beta & 4 \end{pmatrix}\](37)

with \(\alpha = -1 + \delta\), \(\beta = -1 - \delta\), and \(\delta\) takes values 0.0, 0.2, 0.5, 0.8, 5.0, and 8.0. Matrix A is the result of the descritization of a linear differential operator [1]. We use the Orthodir algorithm, or Algorithm 7, in this experiment as a representative of Lanczos-type algorithms. We also use 100 iterations to make up a cycle, which means that the intermediate solution is found after 100 iterations.

Algorithm 7 Orthodir Algorithm [12].

• 1. Initialization. Choose \(\mathbf{x}_0\) and \(\mathbf{y}_0\). Set \(\mathbf{x}=\mathbf{x}_0\), \(\mathbf{r}_0=\mathbf{b}-A\mathbf{x}_0\), \(\mathbf{y}_0=\mathbf{y}\), and \(\mathbf{z}_0=\mathbf{r}_0\).
• 2. for k = 0,1,2,... do
• 3. \(\mathbf{y}_k = A^T \mathbf{y}_k 1\)4. \(A_{k+1} = \frac{\langle -\mathbf{y}_k | \mathbf{r}_k \rangle}{\langle \mathbf{y}_k | A \mathbf{z}_k \rangle}\)
• 5. \(\mathbf{x}_{k+1} = \mathbf{x}_k A_{k+1} \mathbf{z}_k\).
• 6. \(\mathbf{r}_{k+1} = \mathbf{r}_{k} + A_{k+1} A \mathbf{z}_{k}.\)
• 7. Run RLMedVal for k iterations.
• 8. end for
• 9. \(sol_{last} = \mathbf{x}_k\).
• 10. \(norm_{last} = ||\mathbf{r}_k||\).
• 11. Stop.

Table 1 Results of RLLastIt, RLMinRes, and RLMedVal on problems with = 0.0.

Table 2 Results of RLLastIt, RLMinRes, and RLMedVal on problems with = 0.2.

Table 3 Results of RLLastIt, RLMinRes, and RLMedVal on problems with = 0.5.

Table 4 Results of RLLastIt, RLMinRes, and RLMedVal on problems with = 0.8.

Table 5 Results of RLLastIt, RLMinRes, and RLMedVal on problems with = 5.

Table 6 Results of RLLastIt, RLMinRes, and RLMedVal on problems with = 8.

Figure 2 Performance of RLLastIt, RLMinRes, and RLMedVal on SLE's of several dimensions, for δ = 0.2.

4.1 Discussion

We consider two cases, (i) the cases with δ = 0.0, δ = 0.2, δ = 0.5, and δ = 0.8, leading to SLE's with matrices A having condition numbers 119.9999, 98.6081, 62.2227, and 44.3210, respectively; (ii) the cases with δ = 5 and δ = 8, leading to SLE's with matrices A having condition numbers 27.4932 and 24.4970, respectively. Note that the matrices of case (i) have slightly larger condition numbers compared to those of case (ii). However, according to [17-19], the systems are categorized as well-conditioned.

For the first case, RLMinRes has the best performance in terms of accuracy and stability compared to RLLastIt and RLMedVal. As can be seen in Table 1, RLMinRes consistently meets the tolerance when solving several problems. In fact, it is also more robust than RLLastIt and RLMedVal, where there is no breakdown. These results are similar to those of Tables 2, 3 and 4, where they are related to the cases of = 0.2 , = 0.5, and , = 0.8 respectively. For the second case, the results of which are in Tables 5 and 6, RLMedVal is more accurate than RLMinRes, though it is slower. The worst performance in both cases is that of RLLastIt. In fact, it still suffers from breakdown. The behaviour of the restarting from three different points is represented in Figure 2.

5 Summary

In this study we have discussed restarting Lanczos-type algorithms from three different, reasonable points to combat breakdown when solving SLE's. Algorithms called RLLastIt, RLMinRes, and RLMedVal, corresponding to the three different restartings, have been implemented. The numerical results point to restarting from the iterate with lowest residual norm, i.e. RLMinRes, as the most successful. Restarting from the iterate whose entries are the median values of the entries of previous iterates, i.e. alorithm RLMedVal, is sometimes more robust than RLMinRes. Restarting from the last iterate, or RLLastIt, is the worst in both cases. These results confirm our intuitive assessment of each case. For instance, restarting from the best approximate solution found so far, as in RLMinRes, should perform well. And it does indeed. Note that we have solved problems with dimension up to 70000. These are not trivial problems to solve numerically with any method. Moreover, Lanczos-type algorithms can't solve them without restarting or switching.

Interesting questions following from this work could be to investigate whether the quality of the Krylov subspace in which approximate solutions are generated can be estimated from the starting point. Is there any way to tell whether one Krylov subspace is better than another? Can the length of the cycle be determined accurately to avoid stopping too early before breakdown? These issues are relevant in other contexts too and certainly in the context of switching.

Acknowledgements

This project is supported by the Indonesian Government through the Directorate General of Higher Education (DIKTI) Scholarship 2010.