Randomized Linear Solvers

The use of randomization in asynchronous linear solvers allows for the possibility of statements concerning the rate of convergence to be made. A randomized Gauss-Seidel method was introduced by Leventhal and Lewis (2010) building off of the randomized Kaczmarz algorithm proposed by Strohmer and Vershynin (2009), whereby the decrease in the expected value of the error at each step is bounded. The analysis was generalized by Griebel and Oswald (2012) who also added a new parameter that allows for both over and under relaxation. Both of these studies weight the random selection of row 𝑖 by the size of the element 𝑎_𝑖𝑖 ∈ 𝐴. In the case that 𝐴 has unit diagonal this simplifies to a uniform distribution. More recently, Avron, Druinsky, and Gupta (2015) build upon the analysis by Leventhal and Lewis (2010) and Griebel and Oswald (2012) and explicitly analyze the case of asynchronous computation with a uniform distribution.

All of the methods select the vector component to update from a random distribution instead of either sequentially looping through the available components or by tying the updates for a single component to a particular processor (see Equation (2)). In a traditional parallelization of either a synchronous or asynchronous linear solver, processor 𝑗 is responsible for updating component 𝑗; the asynchronous variant allows processor 𝑗 to continue to compute relaxations for the component assigned to it regardless of the state of the other processors. The use of randomization in the selection of which component to update allows for the possibility of any processor updating any component. In a randomized asynchronous linear solver, when a processor finishes computing an update to a component, it writes the update to the shared memory and then randomly draws the next component to update from the list of all available components. In the randomized asynchronous linear solvers proposed by others to date, this random selection is always done using either uniform random number generation, or with a probability proportional to a row norm of the matrix 𝐴. Leventhal and Lewis (2010) cite Fourier analysis as an application area that can benefit from this type of weighting; however, there is no reason not to expect improved behavior for an arbitrary problem. The authors have proposed in Coleman, Jensen, and Sosonkina (2019) to use the non-uniform distributions in the asynchronous Jacobi iterative method. In this work, efficient implementations of such an iterative method are investigated.

Southwell Algorithm

The Southwell algorithm (Southwell, 1946) works similarly to Jacobi by relaxing a single equation at a time, but chooses the equation with the largest local contribution to the residual. For a given row 𝑖, this local contribution is defined to be

at iteration 𝑘. This difference allows the Southwell algorithm to often converge in fewer iterations than Jacobi, but raises the expense of computing an update since the local residuals need to be stored and ranked at each iteration. After a given iteration, the Southwell algorithm chooses the component that contributes the most to the global residual; thus, the algorithm ranks the residuals from largest to smallest. Using the insight from the Southwell algorithm, the idea behind the randomized linear solvers developed here is for each processor to select the next component for updating randomly, using a distribution that more heavily weights selection of components that contribute more to the global residual. Pseudo-code for a randomized variant is provided in Algorithm 2. The key difference of the present work is that here non-uniform distributions in Line 3 of Algorithm 2 are investigated.

In an effort to simulate the effect of the Southwell algorithm using randomized asynchronous solvers, the local residuals associated with each equation (or block of equations) are ranked and sorted, and the selection of the next equation (i.e., component) to update is performed using a non-uniform distribution that forces the random selection to pick components with larger local residuals more frequently. The goal behind the proposed modification is that relaxing the components with a more significant contribution to the global residual may reduce the total number of iterations required. Motivation for this comes from a myriad of different studies, see for instance Nutini et al. (2015) that shows that for some cases (Gauss-)Southwell selection can converge faster than uniform random selection for coordinate descent. In general, the improvement in convergence will have to be shown to be significant enough to offset the extra computational and communication cost associated with storing and ranking all of the local residuals. To help offset the increased computational expense, the periodicity with which the sorting and ranking procedures are done is experimented with since it contributes directly to the overall efficiency of the algorithm.

Problem Description

This work examines the asynchronous Jacobi relaxation algorithm for solving finite-difference discretizations of Partial Differential Equations (PDEs) on a regular grid. In science and engineering, PDEs mathematically model systems in which continuous variables, such as temperature or pressure, change with respect to two or more independent variables, such as time, length, or angle (Smith, 1985). The specific problem under study here is Laplace equation in two dimensions:

where the two-dimensional finite-difference discretization uses Dirichlet boundary conditions. This PDE, which is a fundamental equation for modeling equilibrium and steady state problems, is also used in more complex problems based on PDEs. Equation (3) may be discretized such that a finite difference operator computes difference quotients over a discretized domain. For example, the two-dimensional discrete Laplace operator

approximates the two-dimensional continuous Laplacian using a five-point stencil (Lindeberg, 1990). From this, a discretized version of the Jacobi algorithm

may be applied to solve a two-dimensional sparse linear system of equations (Strikwerda, 2004). Indices 𝑙, 𝑚, and 𝑘 define discrete grid nodes in two dimensions and the iteration number, respectively, for updating the discretized solution vector 𝑣.

In the particular instance of this 2D Laplacian problem, as solved with the Jacobi method here, the grid of 800 × 800 is used to obtain experimental results, the Dirichlet boundary conditions are 100, 0, 75, and 50 for the top, bottom, left, and right boundaries, respectively; the solution vector 𝑣 is initilalized to 0 in each non-boundary grid point, and the right-hand side vector 𝑏 is equal to the initial 𝑣.

Proof-of-Concept

Preliminary experiments are performed using MATLAB^® to demonstrate the improvement in convergence with Southwell and with non-uniform component selection, compared with Jacobi and with uniform component selection, for the problem tested in this work. As an example of potential convergence rates, Figure 1 shows the progression of the residuals over the first 10,000 iterations when solving the two- and three-dimensional finite-difference discretizations of the Laplacian over a 10 × 10 and 10 × 10 × 10 grids, respectively. Here, the four solution methods used are the traditional synchronous Jacobi algorithm, a traditional Southwell algorithm, and two randomized linear solvers: one choosing the component to update using a uniform random distribution, and another using an exponential random number distribution with the parameter λ = 2. Note that the convergence of the randomized linear solver using the uniform distribution is slightly inferior to traditional solvers and to the one with exponential distribution. The latter performs on par with the Southwell, both in the 2D (Figure 1a) and 3D (Figure 1b) cases.

Figure 1. Residual (𝑟 / 𝑟₀) progression for the first 10,000 iterations of four stationary methods solving the 2D (a) and 3D (b) Laplacian.

Block-based Algorithm

The following block-based algorithm design has been introduced in Coleman, Jensen, and Sosonkina (2019) and is provided here as the reference for a wider performance analysis and comparison with the novel, row-based, algorithm. In the task-based asynchronous solver, a thread chooses a block of grid rows to update by sampling from a distribution. The number it draws corresponds to an index in a list of blocks, ranked in order of descending component residuals. For example, if a thread draws the number zero from the distribution, it will update the block-row of components with the largest residual, assuming that block is not being updated by another thread. In the case that a thread selects a block that is already being worked on by another thread, the selecting thread searches sequentially either up or down in the rankings until it finds an available block.

Initially, block residual rankings are assigned via a natural ascending ordering (see Figure 2). A single thread, denoted the residual ranking thread, is tasked with computing the component residuals, sorting the residual rankings, and updating the global ranking list that all the threads use to select blocks for updating. Note that using a single thread leads to a more accurate global ranking list and does not result in a bottleneck for a moderate number of threads. For large-scale distributed implementations, a different ranking procedure has to be developed.

In this work, the residual ranking thread performs ranking and list-updating after every five iterations of the linear system solver. Essentially, Algorithm 2 may be modified to include ranking periodicity τ as shown in Algorithm 3. This ranking period needs to be chosen judiciously, depending on several factors, such as the number 𝑚 of relaxations performed, the number of threads used, and the number $\hat{n}$ of block-rows to rank. Here, τ = 5 was found experimentally to mitigate the ranking overhead for the obtained number of iterations to convergence, while the number of relaxations was varied. A more detailed investigation of the ranking periodicity is warranted and left as future work.

Row-based Algorithm

Algorithm 4 illustrates a novel row-based method. Similarly to Algorithm 3, the master thread periodically, every τ relaxations, ranks and sorts the rows (line 20). However, there are several important distinctions between the two algorithms, due to which Algorithm 4 exhibits better performance. In line 10, a thread uses a probability distribution function 𝑓 to select a single target row to relax instead of a block of rows shown in Algorithm 3, and then transitions from the current (start) row r̃ to the target row 𝑟_𝑣 by relaxing all the rows between r̃ and 𝑟_𝑣 in their natural ordering, instead of jumping to the target row to relax next as done in the block-based implementation (Algorithm 3). Furthermore, while making this transition, a thread may move inward the domain or toward its top or bottom boundary rows, depending on the direction of the shortest distance 𝑑_𝑣 from the current start row to the target (see Equation (4)).

where 𝑛 is the total number of rows in the subdomain, and the direction of progression to the target is toward and across the boundary if the first term in Equation (4) is taken as 𝑑_𝑣; otherwise, the boundary is not crossed. The former is also chosen when the terms are equal. Then, in line 13, the nextr function assigns the next row number to consider by decrementing or incrementing the row number r̃ for the boundary or non-boundary progression direction, respectively; and performing circular shift of the row numbers if they reach the boundary. Note that fewer than 𝑑_𝑣 rows may be relaxed if certain rows in the path towards the target row are not free, i.e., they are already being relaxed by another thread at the time of their consideration, as specified by the conditional statement in line 14. A shared array of size 𝑛 maintains row availability, in which a threads "locks" the row number while it relaxes that row and releases the lock upon finishing the operations in lines 15–20.

The use of the shortest distance is motivated by an attempt to adhere to the ranking order of rows while also relaxing in the neighborhood of the target row; thereby, making the transition to the target smoother. Additionally, in a distributed-memory environment, the ability to more frequently relax boundary rows may facilitate a better data movement among subdomains possibly leading to a faster convergence. Another distinction between the block-based implementation and the row-based one is that the row-based performs the ranking of rows using row-sum differences instead of residuals. In particular (see line 17), after every row r̃ relaxation, a thread performs a summation σ_r̃ of the absolute values of all the components in r̃ and updates the row-sum difference σ_r̃

where ${\sigma }_{\tilde{r}}^{-}$ is the sum taken after the previous relaxation of r̃. This difference σ_r̃ is assumed to be decreasing between the two adjacent relaxations and arbitrarily small when the algorithm has converged. A strong linear relationship has been observed between the row difference method rank and the row residual rank during the entire convergence process. Table 1 presents a small sample of representative correlation coefficients 𝑅 at regular intervals throughout a sample calculation, which quantify the magnitude and direction of this relationship. Of the hundreds of thousands of computed correlation coefficients, the minimum and mean coefficients are 0.77 and 0.96, respectively, with a standard deviation of 0.02. Using this difference instead of residuals decreases the computational overhead of ranking the rows. In particular, finding the row difference requires about 7 times fewer floating point operations per iteration than when using the row residual for this problem. Note, that, while it is shown that the difference-ranking method works for this sample problem, it has not been tested with other types of problems.

Experimental Design

The experiments using OpenMP^® are conducted on two computing platforms at Old Dominion University¹. The Rulfo system has an Intel^® Xeon Phi^™ Knight's Landing 7210 model processor with 64 cores running at 1.30 GHz and 112 GB of DDR4 physical memory used as DRAM in these experiments. One thread per core is utilized, with one core reserved for interfacing with the operating system, resulting in 63 computational threads for the experiments in Section 5.2. On the Wahab system, a single node of the cluster is utilized, containing two Intel^® Xeon Gold 6148 CPUs each with 20 physical cores and 376 GB of DDR4 memory. The code uses standard C++ routines for sorting residuals and generating random numbers, with the default parameters and the built-in distributions. Experimental parameters are presented in Table 2.

Block-based Implementation on Rulfo

For block selection, three different distributions are tested. The uniform distribution is used as a control; a thread may select any block with equal probability. The normal distribution is used to examine the effects of targeting different segments of blocks in the rankings, i.e., blocks with lower ranks and higher residuals versus blocks with higher ranks and lower residuals. This is achieved by varying the mean parameter μ while keeping the standard deviation σ fixed in the normal distribution. The exponential distribution, with the mode λ close to zero, will tend to sample lower-ranked blocks.

For both normal and exponential distributions, the algorithm convergence may be observed in Figure 3 and Figure 4, respectively. In the figures throughout Section 5, the term Recording Iteration points out that the data is recorded by a thread every 1,000 iterations. For the normal distribution, it may be observed in Figure 3a that the convergence rate depends strongly on μ: Its smaller values (up to μ = 46) lead to rapid convergence whereas, at μ = 46, the convergence sharply deteriorates. This may be also observed when considering the time-to-solution in Figure 3b. Due to very slow convergence, at large μ values, the normal distribution becomes extremely non-competitive with the uniform distribution, which timing is shown as red dashed line in Figure 3b. Figure 4a shows that the parameter λ for the exponential distributions does not have as much an impact on performance as the parameter μ does so for the normal distribution runs. As λ moves farther away from zero, however, it hinders convergence and the exponential distribution results in slower timings than those obtained with the uniform distribution as seen in Figure 4b. Once the best parameter choices are found for normal and exponential distributions, their performances compare favorably to the uniform distribution (Figure 5), and up to 10% and 13% fewer iterations are observed, respectively.

Figure 3. BBI convergence for normal distribution.

Figure 4. BBI convergence for exponential distribution.

Figure 6 provides a more detailed explanation for performance differences based on the selection of μ. In particular, Figure 6a and Figure 6b depict that the ordered component residual values for μ equal to 16 and 44 are nearly indistinguishable. However, when μ increases to 48 (Figure 6c) and then again to 52 (Figure 6d) residuals of the lowest-ranked blocks decrease slowly while the residuals of all other blocks decrease more quickly. Note that all the block-based implementation (BBI) experiments use 8 for σ, which is appropriate for all the chosen μ ranges of 16 to 54 on Rulfo and 16 to 40 on Wahab.

Figure 6. Block-row residuals for calculations using normal distributions.

Figure 7a and Figure 7b show that, for the minimum and maximum values of λ, respectively, the component residual decrease is balanced among the component ranks as iterations progress.

Figure 7. Block-row residuals for calculations using exponential distributions.

Row-based Implementation on Wahab

The results of the BBI show that non-uniform probability distribution functions may be used to efficiently select components to update, leading to convergence for the sample problem used in this work. However, relaxing blocks of rows asynchronously tends to cluster errors on block boundaries, and thereby hindering convergence. A row-based implementation (RBI) has been introduced to mitigate this problem. Here, the RBI solves the same sample problem in shared memory as BBI (see Section 5.2). Recall that RBI does not consider blocks of rows to be relaxed by a single thread. Instead, a thread selects only a single row to relax at a time.

For row selection, as with the BBI, the same three distributions are tested. Again, the uniform distribution is used as a baseline for comparison with the normal and exponential distributions. Similar to the BBI experiments, the normal and exponential distributions are geared to consider different ranges of row numbers by, respectively, keeping the standard deviation σ parameter fixed and the parameter λ close to zero. Figure 8a shows the diminishing row differences as the system converges, and the disparity between the rows with the least and greatest differences decreases. In Figure 8b, initially the lowest-index rows have the greatest differences since these are the boundary rows, and in effect, the greatest discontinuity initially is between the top boundary and the first row of grid points (see Section 4.1). Conversely, the least discontinuity initially is between the bottom boundary and the last row of grid points. These respective discontinuities are reflected in the row component differences of consecutive iterations, i.e., the top row initially changes quickly, while the bottom row changes slowly. However, as the calculation progresses, the change in the first rows decreases. For most of the calculation, the middle rows experience the most change.

Figure 8. Progression of row differences and rankings using a uniform distribution.

For the normal distribution, Figure 9 shows the effects of choosing appropriate and excessively large values of the normal distribution mean parameter μ, values of 80 and 400, respectively. Note that the normal distribution standard deviation parameter σ is kept at 40, which is appropriate for the range of μ values considered for RBI here. Compared with Figure 9a, Figure 9b shows increased iterations, greater row-difference disparity between bottom and top-ranked rows, and increasing row differences for ranks 300–400 during the first 1,000 iterations.

Figure 9. Progression of row differences using normal distributions.

Similarly to Figure 8b, Figure 10 shows the rank changes during the convergence processes albeit here for the normal distribution for the same parameters as in Figure 9. In Figure 10a with μ = 80, the middle rows are targeted so frequently that the ranks of the rows with the greatest differences are pushed outward, toward the first and last ranks, much more than what is observed for the uniform and normal with μ = 400 distributions (cf. Figure 8b and Figure 10b, respectively). For μ = 400, because the lower-difference rows are targeted more often, the group of high-ranked rows (shown as the middle yellow band) does not shift ranks to the extent seen with the uniform distribution in Figure 8b, and hence, is updated fewer times, which leads to inferior convergence. This pattern is expected to continue for μ > 400.

Figure 10. Progression of row rankings using normal distributions.

For the exponential distribution, Figure 11 and Figure 12 show the progression of row differences and rankings, respectively. Here, both small and large values of λ provide similar results and are equally effective, on par with good values of μ when using the normal distribution. The convergence history is presented in Figure 13 for the three distributions and their respective parameter choices considered for RBI. As expected, for the exponential distribution and the normal distribution with smaller μ of 80, the residual decreases more quickly than with the uniform distribution, whereas with the normal distribution parameter μ = 400, the residual decreases the most slowly (see Figure 13).

Figure 11. Progression of row differences using exponential distributions.

Figure 12. Progression of row rankings using exponential distributions.

Performance Comparison of Block- and Row-based Implementations

Here, block- and row-based implementations are mutually compared on the same platform, Wahab, as to their number of relaxations and time to converge for a range of non-uniform distribution parameters μ and λ, as shown in Table 2.

Note that the distribution parameters in the row-based implementation differ from those used by the block-based one, which reflects the sorted array sizes and different convergence behavior of the implementations. In particular, for the given test problem, the BBI has 160 entities (blocks) to sort, while there are 800 entities (rows) to sort in the case of RBI. The difference in convergence behavior is especially evident when comparing results from the two implementations when both use normal distributions to select components. In Figure 14a, for BBI, there is a distinct difference in results for μ = 44 and μ = 46. For RBI, Figure 14b shows a smoother transition between good and poor normal distribution parameters. Note that good and poor, respectively, are termed so because they yield the calculation times faster and slower than those for the uniform distribution test cases. In particular, the poor distribution parameters are those starting with the first μ that yields a significant jump in the calculation time; and this percentage increase for RBI is taken to be comparable with the one in BBI. By comparing the results for the μ values in Figure 14a and Figure 14b, it is seen that the RBI tolerates a much higher relative value for μ than the BBI does so before significantly degrading the performance. For example, while μ = 46 is already a poor choice for the BBI, μ = 230 (which is equal to 46 × 5 rows in a block) is still well within the range of good parameters for the RBI.

Figure 14. The total number of all the grid-component relaxations until convergence, for different probability distributions and parameters. The red lines refer to the number of component relaxations for serial Gauss-Seidel.

In addition, Figure 14a and Figure 14b compare the block- and row-based implementation iterations with the iterations of serial Gauss-Seidel (shown as red horizontal lines), respectively, to converge for the sample problem. The BBI cannot converge in fewer than serial Gauss-Seidel component relaxations even with the best distribution parameters. The RBI, however, converges in about 10% fewer component relaxations than serial Gauss-Seidel, using non-uniform distributions with appropriate parameters. This happens consistently, although it has been shown theoretically that more component relaxations may be required when threads update components asynchronously (Avron et al., 2015). A better convergence in the RBI compared with that in BBI may be attributed to the (fine-grained) ranking of rows rather than blocks and to relaxing all the rows on the path from the current and the selected target one. Such a relaxation process leads to a smoother transition between rows and possibly to relaxations of more rows by a thread at a time than those contained in a block of the BBI. Although the row-based implementation ranks and sorts more entries than the BBI does so, the former has faster time-to-solution (see Figure 18) and is not hindered at large scales—where distributed implementations are a must—because ranking and sorting will be performed within each node independently.

Complementing the convergence comparisons of BBI and RBI from Figure 14, Figure 15 demonstrates (as vertical lines in each bar) a greater variability in how often each block in BBI may be relaxed compared with each row relaxation in RBI. This metric bears significance for the non-uniform distributions since they may "neglect" certain components to relax often enough to hinder convergence, as has been shown earlier in Section 5, and thereby making a proof of convergence more difficult.

Figure 15. The average number of (a) block and (b) row relaxations required to converge for different probability distributions and parameters for the two implementations. The vertical lines in each bar show the standard deviation of the number of row relaxations among all rows.

Figure 16 and Figure 17 compare BBI and RBI as to which parts of the problem grid are relaxed more times when good or poor μ is used, respectively. For the former, Figure 16 shows not only that both implementations emphasize the relaxation of the middle rows, away from the fixed top and bottom boundaries, but also that the RBI places greater emphasis on the rows near the top and bottom boundaries, and less emphasis on the middle rows, compared to the BBI. In particular, about 15% of component selections result in a boundary-crossing event in the row based implementation, which provides for relaxing all the rows more uniformly. With poor distribution parameters, Figure 17 shows a different behavior of the RBI from the one in Figure 16. Now, the RBI relaxes boundary rows more frequently than it does so for the innermost rows. In particular, some of the inner rows are now relaxed about as many times as for good μ but the boundary rows are relaxed more frequently leading to an overall higher number of iterations to converge. Generally, the RBI permits more frequent relaxation of boundary rows, compared with the BBI. Note that the frequency of boundary-row relaxation stays low for BBI given either value of μ (cf. Figure 16 and Figure 17). Such a beneficial behavior of RBI is expressed in line 13 of Algorithm 4, in which the nextr function directs a thread to or from a boundary row according to the shortest distance (line 11) as determined by Equation (4).

Figure 18 shows that RBI decreases calculation time (Figure 18b) compared with BBI (Figure 18a) for all the distributions on the Wahab cluster. Furthermore, a 10% convergence-time reduction is observed for the row-based implementation using normal and exponential distributions with good parameter choices, as compared to a uniform distribution. Figure 18b shows a gradual increase in calculation time for increasing values of μ beyond 200, similar to the gradual increase in numbers of relaxations seen in Figure 14b and Figure 15b. For BBI on the Wahab platform (Figure 18a), the results show a jump in calculation time when the normal distribution is used, which is also observed on Rulfo (cf. Figure 3b) albeit at a larger μ value of 46. On Wahab, the BBI threshold μ is 40, which suggests that, for the normal distribution shared-memory implementation, good parameter selection is platform-dependent, as expected. In particular, having more threads results in smaller size blocks, which may mitigate poor μ selection in the BBI.

Figure 18. Wahab calculation times for each implementation and all three distributions. Note the λ-labeled axis pertains to the exponential distribution trajectory while the μ-labeled axis refers to the normal distribution trajectory.

In addition to the performance benefit seen with the row-based implementation, Figure 19 and Figure 20 illustrate that the RBI produces a solution with the residual values more uniformly dispersed among all components. For each implementation, the plots display the two runs with the smallest and largest maximum component residuals, out of a set of ten runs that use the exponential distribution with λ = 0.05 for BBI and λ = 0.01 for RBI. The BBI gives a mean maximum component residual of 4.3𝑒^-11, with a standard deviation of 2.2𝑒^-11, while the row-based implementation gives a mean of 1.0𝑒^-11 and a standard deviation of 1.6𝑒^-12. Note that the largest maximum component residual produced by the RBI, as seen in Figure 20b, is about half the size of the smallest component residual produced by the BBI, as seen in Figure 19a. Observe also that the variations between runs are less for RBI than they are for BBI.

Figure 19. BBI solution component residual values from the runs with the smallest and largest maximum component residuals, for exponential distribution, λ = 0.05.

Figure 20. RBI solution component residual values from the runs with the smallest and largest maximum components residuals, for exponential distribution, λ = 0.01.