A Feasible Approach to Determine the Optimal Relaxation Parameters in Each Iteration for the SOR Method

The paper presents a dynamic and feasible approach to the successive over-relaxation (SOR) method for solving large scale linear system through iteration. Based on the maximal orthogonal projection technique, the optimal relaxation parameter is obtained by minimizing a derived merit function in terms of right-hand side vector, the coefficient matrix and the previous step values of unknown variables. At each iterative step, we can quickly determine the optimal relaxation value in a preferred interval. When the theoretical optimal value is hard to be achieved, the new method provides an alternative choice of the relaxation parameter at each iteration. Numerical examples confirm that the dynamic optimal successive over-relaxation (DOSOR) method outperforms the classical SOR method. Keyword: linear equations system, successive over-relaxation (SOR) method, maximal projection, dynamic optimal relaxation parameter


Introduction
In the paper, we derive a better realization of the successive over-relaxation (SOR) method [Quarteroni, Sacco & Saleri (2000)] to find unknown variables x ∈ R n from the linear equations system: where A ∈ R n×n is a given non-singular coefficient matrix and b ∈ R n is a given right-hand side vector.
In general, f 0 > 1 because x does not exactly satisfy Equation (1) during the numerical process. When f 0 = 1 is reached, x exactly satisfies Equation (1) and meanwhile the solution is obtained. Liu (2013cLiu ( , 2014b) employed a scaling invariant property of Equation (2) to derive a maximal projection solution in the Krylov space and proved that Equation (2) implies the least squares solution [Blais (2010)]: Liu (2014c), based on Equations (2) and (5), has developed a double optimal technique for the solution of Equation (1).
The SOR is a well-known and well-developed classical iterative method, and the formulation depends on a relaxation parameter, whose optimal value needs to compute the spectral radius, defined as the absolute value of the largest eigenvalue in magnitude of the iteration matrix. Therefore, its computational burden to involve the computation of optimal relaxation parameter is great, when large scale linear problem is considered. Only for limited cases, the theoretical value of the optimal relaxation parameter is known. Bai and Chi (2003) have chosen the optimal relaxation parameter by the asymptotically optimal successive over-relaxation method in a dynamic fashion. Wen, Meng & Wang (2013) have obtained the optimal parameter by an optimization technique based on the quasi-Chebyshev accelerated iteration method. Similarly, Meng (2014) has proposed another asymptotic approach of the optimal relaxation parameter. On the other hand, Miyatake, Sogabe & Zhang (2020) proposed an adaptive method based on the Wolfe condition to search the optimal relaxation parameter.

Successive Over-Relaxation (SOR) Method
It is known that the matrix A can be uniquely decomposed into where the components of these matrices are given by D is a diagonal matrix, U is a strict upper triangular matrix, while L is a strict lower triangular matrix.
From Equations (1) and (6) it follows an equivalent linear system: Multiplying the above equation by a nonzero constant w, it becomes Then adding a term Dx on both sides, we have By using an iterative method to solve the linear system, we suppose that the value x k at the kth-step is known. Now, we remove wDx and −wUx to the right-side and take the value of x in the the right-side to be x k , while in the left-side we take the value of x to be x k+1 , and then the SOR iterative method is created as follows [Hadjidimos (2000)]: of which 0 < w < 2 is a relaxation parameter to guarantee the convergence of iteration. If D ii 0, i = 1, . . . , n, from Equation (13) it is easy to find x k+1 by using the forward substitution method, because D − wL is a lower triangular matrix.
When A is a positive definite matrix, the following w is the optimal one [Quarteroni, Sacco & Saleri (2000)]: where ρ is the spectral radius of I n − D −1 A.

Maximizing Orthogonal Projection
In the classical SOR method, one needs to compute the spectral radius of I n − D −1 A as shown in Equation (14), which requires a huge amount of calculation work. In general, the approximate manner to try w and to observe the convergence effect is adopted [Milewski & Orkisz (2014); Yang & Gobbert (2009)].
The minimization in Equation (2) is equivalent to maximize which is the orthogonal projection of b/∥b∥ onto the y direction.
In order to obtain the optimal value of w in Equation (13), we temporarily take x k+1 in the left-hand side to be x k , and then at each iterative step, we have a linear system: and insert them into Equation (2), we can derive At each iterative step, we can search the optimal value of w k to minimize f 0 : We can quickly find the optimal value of w k in a preferred interval (a, b) for the convergence of the SOR. In this study, we use this process to pick up w k and the resulting iterative algorithm is very time saving and is named a dynamic optimal SOR (DOSOR) method. Instead of using the constant value of w in the SOR method, the present method in Equation (24) provides a dynamic value of w k at each iteration.
the iteration process is terminated for satisfying a given convergence criterion. As mentioned below Eq. (4), when f 0 −→ 1, the solution is obtained. In this regard, we can also take the following as a convergence criterion: where ε 1 and ε 2 are small positive numbers.

Numerical Verification
In order to compare the performance of the DOSOR to the SOR, we test some linear problems.

Example 1
Equation (1) is solved with From Equation (14) w opt is found to be w opt = 1.016288735, and the SOR is convergence with 28 iterations as shown in Figure 1(a) with the following initial guess: The maximum error (ME) is 5.71 × 10 −11 and the root-mean-square-error (RMSE) is 2.98 × 10 −11 . Because w opt is given with one calculation, the CPU time of the SOR is very short with 0.2 s.
By using the DOSOR, we take a = 0.9, b = 1 and N w = 10, and through 26 iterations it converges under ε 1 = 10 −10 as shown in Figure 1 (b). Because N w is small with N w = 10 calculations to determine w k at each iteration, the CPU time of the DOSOR is still very short with 0.25 s. In Figure 1(b), the values of f 0 are plotted, which fast tends to the minimal value f 0 = 1, and in Figure 1(c), the values of w are plotted, which are varying between [0.9, 1], where the peak value is close to the optimal value given in Equation (28). The ME obtained by the DOSOR is 2.41 × 10 −11 and the RMSE is 1.27 × 10 −11 , which are better than that obtained by the SOR.
If we raise N w to N w = 100, the CPU time of the DOSOR is slightly increased to 0.3 s, and the ME and the RMSE become 4.85 × 10 −11 and 2.35 × 10 −11 , respectively. Larger N w would increase the CPU time but not necessarily serves better accuracy. If we change the right-hand side to the number of iterations, the CPU time and the accuracy of the DOSOR are not affected apparently. Of course, for different coefficient matrix A, the number of iterations, the CPU time and the accuracy of the DOSOR are different as to be shown below.

Example 2
In this example, we apply the DOSOR to solve the following boundary value problem: The exact solution is supposed to be and f (x) = −π 2 sin πx.

Residual
If we raise the value of N w = 100, the CPU time increases to 2.3 s, but the accuracy is not affected with ME=2.71 × 10 −5 and RMSE=1.48 × 10 −5 . Larger N w would increase the CPU time but not necessarily offers better accuracy. So we prefer to use N w = 10 for saving CPU time. If we adopt u(x) = x 3 − x 2 as being another solution, the new right-hand side is obtained to be f (x) = 6x − 2. By using the DOSOR, we take n = 99, a = On the other hand, we have [Demmel (1997); Watkins (2002); Yang & Gobbert (2009)] which is inserted into the SOR to find the solution. It converges through 202 iterations, faster than the DOSOR; however, the maximum error obtained by the SOR with the above w opt is 9.27 × 10 −5 and the RMSE is 6.74 × 10 −5 , which are less accurate than that obtained by the DOSOR with ME=2.11×10 −5 and RMSE=1.4×10 −5 . Because w opt is given by Eq. (33) with one calculation, the CPU time of the SOR is very short with 0.5 s. The CPU times of the SOR and the DOSOR are competitive if N w = 10 is taken in the DOSOR. The accuracy of the presented method is enhanced for the use of the merit function (24).

Example 3
We consider the Hilbert matrix:  (1), which is a notorious example of highly ill-conditioned matrices. Also, x j = 1, j = 1, . . . , n is an exact solution, leading to By using the DOSOR, we take n = 99, a = 0, b = 2 and N w = 50, and through 691 iterations it converges under ε 1 = 10 −4 as shown in Figure 3(a). The CPU time of the DOSOR is quite short with 1.7 s. In Figure 3(b), the values of f 0 are plotted, which tends to the minimal value f 0 = 1 very fast, and in Figure 3(c) the values of w are plotted, which are varying between [0.04, 0.24] and tend to 0.04 very fast. The ME obtained by the DOSOR is 1.5 × 10 −2 and the RMSE is 4.56 × 10 −3 . With an ad hoc value of w = 1.5, the number of iterations becomes 3297 and the ME and RMSE raise to 7.3 × 10 −2 and 1.36 × 10 −2 , respectively.

Example 4
In this example, we apply the DOSOR to solve the following 2D boundary value problem of the Poisson elliptic equation: While is an exact solution, p(x, y) = 2(x 2 − x) + 2(y 2 − y) is derived.
As shown by Bai & Chi (2003), we take n 1 = n 2 = 31, ∆x = ∆y = 1/32 and ε 1 = ∆x 2 ∥r 0 ∥/5, where r k = b − Ax k is the residual vector. When is inserted into the SOR, it converges through 62 iterations as shown in Figure 4(a). The CPU time of the SOR is very short with 0.25 s. By using the DOSOR, we take a = 1.8, b = 2 and N w = 10, and through 61 iterations it converges as shown in Figure 4(a) with the CPU time being 6.86 s. The ME obtained by the SOR is 3.66×10 −6 and RMSE=1.42×10 −6 , which is less accurate than that obtained by the DOSOR with ME=1.52 × 10 −6 and RMSE=5.59 × 10 −7 . The accuracy of the presented method is enhanced owing to the use of the merit function (24).

Conclusions
There are two factors to evaluate the performance of a newly developed iterative scheme: accuracy and convergence speed.
In this paper, we proposed a dynamic optimal SOR method based on the maximal orthogonal projection technique, which is equivalent to the minimization in Equation (24). The features of the proposed method are summarized as follows: (a) The merit function is in terms of the coefficient matrix, right-side vector and the value of unknown vector at the previous step. (b) Searching the minimization in a preferred interval through a few operations is easily performed. (c) The accuracy is better than the classical SOR with optimal relaxation parameter about two to four times in the maximum error as well as in the root-mean-square-error. (d) The convergence speeds are competitive of both methods. The CPU time of the DOSOR with low number of N w is slightly increased than the SOR. (e) When the theoretical value of the optimal relaxation parameter is in general not available, the new method provided an alternative and feasible choice of the relaxation parameter at each iteration.