A NEW REPROJECTION OF THE CONJUGATE DIRECTIONS

. In the paper we apply a Parlett–Kahan’s “twice is enough” type algorithm to conjugate directions. We give a lower bound of the digits of precision of the conjugate directions, too.


1.
Introduction. In our earlier paper [2] we applied the Parlett-Kahan's "twice is enough" method [7] for the Gram-Schmidt type algorithm, for reprojection of conjugate directions in the ABS class. We confirmed the correctness and applicability of the developed theory by two test problems. The second problem was defined by the well-known difficult Pascal matrix. Not knowing the diagonal elements of the P T AP matrix, where P contains the conjugate directions with respect to the positive definite symmetric matrix A, we could give the relative accuracy of the conjugate directions in the S2 ABS subclass (see [3]).
In the meantime, Hegedűs [5] in his paper "Reorthogonalization Methods Revisited" developed an alternative variant of the Parlett-Kahan's "twice is enough" algorithm. This algorithm makes it possible to apply it to the reprojection of the conjugate directions in the S2 subclass. So we can give the accurate number of digits of the conjugate directions. As in case of any n dimensional problem there are n different conjugate directions, all of them with a different number of accurate digits of precision. We always present the worst precision for the sake of brevity.
In the next section we define the problem to be solved and give the theoretical background of the "twice is enough" algorithm for conjugate directions. In Section 3 we give three preliminary test results for a very simple problem when practically there is no need for reprojections and two famous examples of ill-conditioned matrix like the Hilbert and Pascal matrices which demonstrate the efficiency of the reprojections as the number of accurate digits of the conjugate directions is practically coincides with that of all digits representable in double precision, that is, the accuracy can not be improved upon.
where D is a full rank diagonal matrix and P = [p 1 , ..., p n ]. Determine the accurate digits of these conjugate directions without knowing matrix D. This means that the algorithm is to compute the number of accurate digits during its execution. Note that this is important because the matrix D can contain rounding errors accumulated during its computation.
In the paper [5] Hegedűs gave a new explanation of the Parlett-Kahan's "twice is enough" algorithm [7] for reorthogonalization of the Gram-Schmidt type methods. We have to note that both in the Parlett-Kahan and in the Hegedűs's approach it is assumed that all values are accurate and the effect of the cancellation of error in the orthogonal Gram-Schmidt type computations is examined. We suppose the same here below.
In this paper we apply the Hegedűs type new approach to the ABS type conjugate gradient methods. These methods are defined in Subclass S2. The scaled ABS algorithm goes as follows.
Step 1 Set x 1 ∈ n , H 1 = I ∈ n,n where I is the unit matrix, i = 1, and if lag = 0.
Step 2 Let v i ∈ n be arbitrary and let v 1 , ..., v i be linearly independent. Compute the residual error vector Step 4 Compute the search direction p i by Step 5 Update the approximation of the solution where the step size α i is given by If i = m, stop and x i+1 is the solution of the equations.
Step 6 Update the matrix H i by where w i is arbitrary but the denominator must be non-zero.
Step 7 Set i = i + 1 and go to Step 2.
The properties of this class of methods can be found in [1]. Let A ∈ n,n be a symmetric positive definite matrix. Then the S2 subclass is defined by In this case the vectors {p 1 , ..., p n } are A conjugates see Theorem 8.6 p. 123 of [1]. Define now the free vector w i in the following way and furthermore let With these definitions we can write where .., n, H 1 = I, where I is the unit matrix and from the fact that coefficient matrix A is symmetric positive definite.
Note that Theorem 7.3 is much more general than we need. We only require the vectors q i i = 1, ..., n to be orthonormal.
Consequently, with this algorithm we constructed conjugate directions p i i = 1, ..., n and orthonormal vectors q i i = 1, ..., n by (4) and the vectors z i , i = 1, ..., n are still free. We can apply formula (5) of [5] for a = z i to construct the next orthogonal vector to Q i−1 . That is, we get .., n.From this p i , i = 1, ..., n we compute the next q i , i = 1, ..., n by (4) . Note that vectors p i , i = 1, ..., n are A conjugate.
Introduce in every step see [5]. Then Hegedűs proved the following Theorem 2 [5], If there are accurate digits in the computation, then one may expect the fulfillment of condition η max ≤ η after the second orthogonalization step at most. The largest choice of such η max is 1/ √ 2 to fulfill the condition. Hence the resulting vector p i can be considered orthogonal {q 1 , ..., q i−1 } up to computational accuracy if η max is not less than 0.5.
The number of digits of precision can be calculated by (15) of [5] that is In the sequel we will take getting thereby the worst number of digits of the precision of the conjugate directions p k . Based on Theorem 3 he developed the so-called orthogonal base algorithm OBA see p. 12 of [5]. The next algorithm follows using the two theorems above. We suppose that the free parameters v i and w i are defined by (2) and ( Note that 4eps mach where eps mach is the machine epsilon is usually chosen because of the rounding error. We have to mention that in our earlier paper [2] for η max a different value suggested in Parlett's work [7] was chosen. We have written an other paper in which we tested the algorithms developed in S2 [3]. These preliminary tests show the usefulness of reprojections. In this paper the value of η max is the theoretically obtained 1/ √ 2 and we now can determine the lower bound of the accurate digits in the p i vectors, which are A conjugate.
3. Preliminary test results. In this section we present the results of a preliminary test of the CDRA algorithm in the S2 ABS subclass. We give after each figure the number of reprojections and the maximum absolute norm of residual. The algorithms were implemented in MATLAB version 2007b.
Both for the η dependent linear dependency and the ABS linear dependency were controlled by 4eps mach .
The theoretical result of the previous section determines exactly the free parameters v i and w i of the scaled ABS class, while the parameter z i remains free. Note that the projection matrix H i is symmetric. Therefore we have chosen the following algorithms which are defined by the free vector z i . We have to mention that we use the idempotent property of the projection matrices H i therefore we define the w i vectors too even if we know that they are theoretically equal. The name of the actual algorithm is given in parentheses.
1) z i = r i and w i = H i Ap i (S2HSsz).
2) z i = r i and w i = Ap i (S2rsz). Note that theoretically the same as S2HSsz.
3) z i = a i and w i = H i Ap i (S2asz). 4) z i = a i and w i = Ap i (S2a). Note that theoretically the same as S2asz.
Note that theoretically the same as S2LU.
The reason for these choices is the fact that we frequently observed, that if in the dyad of the projection matrices H i reprojection is done then the accuracy of the projected vectors is more accurate.
In the sequel we present for all algorithms a figure which shows the dimension versus the digits of precision of p after reprojection and in the same row the figure of the parallel algorithm can be found. Consequently it is easy to verify the effect of the reprojections in the dyad. As in every step of the algorithms reprojection is possible, after the figures we give the number of reprojections in a given dimension and the number of linear dependency too. We have to mention that the zero value of the precision of the conjugate direction p in the figure means no reprojections.
Golub and Van Loan [4] on p. 60 of their book say "...we cannot justifiably criticize an algorithm for returning an inaccurate x if A is ill conditioned relative to the machine precision e.g. uκ ∞ (A) ∼ = 1." where x is the approximate solution u is the unit roundoff defined in (3.2.2) (p. 33) and κ ∞ (A) is the condition number of A in infimum norm. The unit roundoff u in MATLAB is the eps value because eps returns the distance from 1.0 to the next double-precision number, that is eps = 2 −52 . Therefore in the sequel we compute uκ ∞ (A) for the different dimensions and these values determine the interval where the conjugate vectors should be computed.
As we are interested in the effect of the reprojections in the figures we only present the lower estimation of the numbers of digits of precision of p after reprojection computed by (8). Therefore where the number of digits is zero in the figures, it means that it was not necessary to reproject the conjugate vector p. So we get information from the reprojection algorithm. Also we give the maximum norm of the residuals which indicates how accurate the solution of the linear system of equations is.

3.1.
Test results with an easy problem. This problem is practically a FOR-TRAN test problem type matrix of LAPACK. The positive definite symmetric coefficient matrix A is generated randomly in the interval [0, 1] by MATLAB rand function in dimensions 3000-3010 and also the solution of the linear system is generated randomly in [0, 1] by MATLAB rand function. The uκ ∞ (A) values in these dimensions are 6.974e-016 7.012e-016 7.011e-016 7e-016 6.998e-016 6.984e-016 6.999e-016 6.992e-016 6.981e-016 6.998e-016 6.997e-016, which means really a very easy problem. Figure 1a and 1b based on computations without CDRA in order to demonstrate the effectiveness of the CDRA algorithm numerically as well. In Figure 1 vectors w i are reprojected in the dyads (see the definition above), while in the second one no reprojection was done. Figure 1a and 1b clearly show that the number of exact decimal digits is essentially zero, while the additional diagrams show the effectiveness of the CDRA algorithm that is the number of the correct digits is often all the representable decimal digits. At the end of the paper we choose the best ones from the presented algorithms.
The maximum norm of the residuals is 4.391e-010 and 4.428e-010, respectively. These good residuals are correct because we show the worst digits of precision of p. We also have to mention that in figure S2rsz the zero values are not really zeros. For example at dimension 3004 it is 2.5495 −9 . Note that the worst conjugate directions have zero accurate digits. This will be improved to near the double precision accuracy in the sequel. Now we turn to the cases computed with CDRA. We begin with the same two algorithms.
The    η or in the ABS control. The maximum norm of the residuals was found 4.473e-010, very good again.
In case of S2asz there was no reprojection. There was no linear dependence either in the reprojection η or in the ABS control. The maximum norm of the residuals was found 5.574e-010.
In case of S2a everything is the same as in S2asz except the maximum norm of the residuals is 4.571e-010.
The P matrices are full in both cases. Both in case of S2LU and S2esz everything is as above. In case S2LU the maximum norm of the residuals is 4.148e-010. and 3.96e-010 at S2esz. The P matrices are full in both cases.
In case of S2Hpsz the numbers of the reprojections are 6 5 9 6 11 11 8 6 7 8 2 and the maximum norm of the residuals is 4.159e-010.
In case of S2psz the numbers of the reprojections are 6 5 9 6 11 11 8 6 7 8 2, the same as before while the maximum norm of the residuals is 4.195e-010.  The P matrices are full in both cases. After these 3000 dimensional randomly generated easy problems we conclude the followings: (a) Practically there is no difference between the reprojection in the dyads and the lack thereof. (b) There were four methods where no reprojection was done at all. (c) The algorithms with this reprojection technique work very well.
From the table it follows that we need to consider the dimensions between [5,15] only. In addition, we present results in the whole interval [5,50] showing that accurate conjugate directions can even be obtained over dimension 15 as well. Again, the solutions are generated randomly here too.
In this case the effect of the reprojections in the dyads can be seen very clearly because the computation of the conjugate vectors p i does not produce any errors.
The results of S2psz as we can see in the figure are very similar to those in the case S2Hpsz. Only in the second part can any difference be seen. The numbers of the reprojections are small again 0 0 0 1 0 0 0 0 1 0 1 1 1 1 1 1 1 2 0 2 0 2 0 1 3 0 0 0 1 1 1 0 0 1 2 1 2 1 0 0 1 0 2 1 1 2. There was no linear dependence either in the reprojection η or in the ABS control. The maximum norm of the residuals was 9.247e+016.
In both cases the number of the exact digits are excellent in every dimensions. Practically all digits are accurate. These two algorithms are very interesting because in Theorem 7.10 of [1] the ith conjugate direction equals to which practically means that all the computations in these two algorithms strongly depend on the first residual vector. Clearly, we should study this phenomenon in the future.
We conclude the followings: (a) There is a marked difference between the versions when reprojection is done in the dyads and when it is not. We further narrow the number of the good methods with the Hilbert problem.

3.3.
Test results with the Hilbert problem. We examine the Hilbert problem in [5,50] [5,11]. Here again we show results for the whole interval [5,50].    We conclude that even if the solutions are very nice, the P matrices contain many zero vectors, therefore from our point of view these algorithms are not acceptable.
We conclude that both figures show conjugate directions with all accurate digits, the solution of the linear systems are very small, and the P conjugate direction matrices are full.
We conclude this section that definitely only the last two algorithms give acceptable results for the Hilbert problem according to our aim. 4. Conclusion. In this paper we presented a new algorithm for the reprojection of the conjugate directions. The theory developed make it possible to estimate the exact number of the accurate digits in the conjugate directions. The preliminary test results underline the usefulness of the new reprojection algorithm in the S2 subclass of the ABS methods. Our first easy problems show that in case of simple problems there is no need for any reprojections in case of 4 algorithms. The Pascal problems show that reprojections were always necessary. The very difficult Hilbert problem gave the final selection among the 8 algorithms. The S2Hpsz and S2psz algorithms give full P matrices with all accurate digits. Notice that the last two methods show another interesting problem, where except for the first step in the computation of the conjugate directions the previous direction was used in all steps and therefore everything directly depends on the first residual.
From the results of the test problems section we recall that the worst absolute residuals were obtained in the case of the Pascal problems but this does not render our algorithms inefficient because the problems are ill-conditioned, see p. 60 of [4]. For the easy 3000-3010 dimensional problems and the Hilbert's problem the residuals are acceptable.
Finally we underline that for all test problems we got practically accurate values for all decimal digit and full P matrices in case of the last two algorithms S2Hpsz and S2psz.
We encourage further research since many interesting problems have remained unsolved.