Combining Sparse Approximate Factorizations with Mixed-precision Iterative Refinement

We want to determine the convergence conditions for LU-IR (Algorithm 1). We can apply [18, Corollary 3.3] and [18, Corollary 4.2], respectively, for the convergence of the forward and normwise backward errors of the system \(Ax=b\), and for both, we need, respectively, a bound on the forward and backward errors of the computed solution \(\widehat{d}_i\) of the correction equation \(Ad_i=\widehat{r}_i\). Note that for LU-IR, the factorization (4.6) and the LU solves (4.7) are performed in precision \(u_f\).

Considering the solution of the linear system \(Ad_i=\widehat{r}_i\), (4.8) yieldsTaking norms, we obtainUsing [32, Lemma 9.6]where \(\rho _n\) is the growth factor for \(A + \Delta A^{(1)}\). Dropping second-order terms finally givesIn the same fashion, we can show that

Dropping constants and applying [18, Corollary 3.3] and [18, Corollary 4.2] using (4.10) and (4.11) guarantees that as long asthe forward and the normwise backward errors of the system \(Ax=b\) will converge to their limiting values (4.4) and (4.5).

As a check, if we set \(\epsilon =0\) (no approximation) and drop \(\rho _n\) (negligible element growth), then we recover (4.1).

Before commenting in Section 4.2.3 on the significance of these new LU-IR convergence conditions, we first similarly derive the GMRES-IR conditions.

We now determine the convergence conditions of GMRES-IR (Algorithm 2). We proceed similarly as for LU-IR, seeking bounds on the forward and backward errors of the computed solution \(\widehat{d}_i\) of the correction equation \(Ad_i = \widehat{r}_i\). One difference lies in the fact that the GMRES solver is applied to the preconditioned system \(\widetilde{A}d_i=\widehat{s}_i\) where \(\widetilde{A}= \widehat{U}^{-1}\widehat{L}^{-1}A\) and \(s_i = \widehat{U}^{-1}\widehat{L}^{-1}\widehat{r}_i\). Our analysis follows closely the analysis of [3, Section 3.2], so we mainly focus on the changes coming from the use of a more general approximate factorization model and refer the reader to that work for the full details.

We first need to bound the error introduced in forming the preconditioned right-hand side \(s_i\) in precision \(u_p\). Computing \(s_i\) implies two triangular solves (4.7), which differ from the original GMRES-IR analysis by having an error term on the right-hand side. Adapting [3, (3.11)–(3.13)] with (4.6) and (4.7) and using (4.9) provides the bound

As in [3], we compute the error of the computation of the preconditioned matrix–vector product \(z_j = \widetilde{A}\widehat{v}_j\) to use [3, Theorem 3.1]. We obtain \(z_j\) through a standard matrix–vector product with A followed by two triangular solves (4.7) with \(\widehat{L}\) and \(\widehat{U}\). The computed \(\widehat{z}_i\) satisfies \(\widehat{z}_i = z_j + f_j\), where \(f_j\) carries the error of the computation. With a very similar reasoning as for deriving [3, (3.14)], considering our new assumptions, we obtain the boundApart from the constants and the presence of the growth factor \(\rho _n\), which can be arbitrarily large without assumptions on the pivoting, (4.14) and (4.13) are similar to [3, (3.14)] and [3, (3.13)] and meet the assumptions of [3, Theorem 3.1], which can be used to compute a bound of \(\Vert s_i - \widetilde{A}\widehat{d}_i\Vert _{\infty }\).

We can finally bound the normwise relative backward error of the system \(\widetilde{A}\widehat{d}_i = s_i\) [3, (3.17)] byand the relative error of the computed \(\widehat{d}_i\) [3, (3.18)] byIn addition, the backward error of the original correction equation \(Ad_i=\widehat{r}_i\) can be bounded using \(\widehat{r}_i - A \widehat{d}_i = \widehat{L}\widehat{U}(s_i - \widetilde{A}\widehat{d}_i)\) and (4.15), yielding

It is essential to study the conditioning of the preconditioned matrix \(\widetilde{A}\) to express the convergence conditions according to the conditioning of the original matrix \(\kappa (A)\). Using the same reasoning as for [17, 3.2], we obtain

Dropping constants and applying [18, Corollary 3.3] and [18, Corollary 4.2] using (4.16) and (4.17) guarantees that as long asthe forward and the normwise backward errors of the system \(Ax=b\) will converge to their limiting values (4.4) and (4.5).

As a check, with \(\epsilon =0\) and dropping \(\rho _n\), we recover (4.2) and (4.3).

We summarize the analysis in the following theorem.

We now comment on the significance of this result. Compared with the original convergence conditions (4.1)–(4.3), the new conditions of Theorem 4.1 include two new terms. The first is the growth factor \(\rho _n\) that, without any assumption on the pivoting strategy, cannot be assumed to be small. This shows that a large growth factor can prevent iterative refinement from converging. The second is \(\epsilon\), which reflects the degree of approximation used by the factorization. The terms \(\rho _nu_f + \epsilon\) show that we can expect the approximation to impact the convergence of iterative refinement when \(\epsilon \gtrsim \rho _n u_f\) (ignoring the difference between the constants in front of each term). It is important to note that the instabilities introduced by element growth and numerical approximations combine additively, rather than multiplicatively (there is no \(\epsilon \rho _n\) term). In particular, this means that the usual wisdom that it is not useful to use a very high precision for an approximate factorization (\(u_f \ll \epsilon\)) is no longer true in presence of large element growth. This is a key property that we confirm experimentally in Section 5.

Note that the left-hand side quantities in the convergence conditions (4.20)–(4.23) also bound the convergence rate of the refinement: thus, smaller quantities will in general lead to faster convergence.

We now apply the above analysis to the use of BLR approximations and static pivoting.

The BLR format approximates the blocks of the matrix by replacing them by low-rank matrices. The ranks are determined by a threshold, \(\tau _b\) in this article, that controls the accuracy of the approximations. Higham and Mary [35] have carried out error analysis of the BLR LU factorization and obtained a backward error bound of order \(\tau _b\Vert A\Vert +u_f\Vert \widehat{L}\Vert \Vert \widehat{U}\Vert\). One issue is that their analysis derives normwise bounds, whereas our model (4.6) and (4.7) require componentwise bounds. However, we have checked that, at the price of slightly larger constants by about a factor \(O(n^{3/4})\) and a more complicated analysis, analogous componentwise bounds can be obtained. Therefore, using the componentwise version of [35, Theorem 4.2] or [35, Theorem 4.3] for (4.6) and of [35, Theorem 4.4] for (4.7), we conclude that Theorem 4.1 applies with \(\epsilon =\tau _b\).

We now turn to static pivoting, assuming a strategy that replaces pivots smaller in absolute value than \(\tau _s\Vert A\Vert _{\infty }\) by \(\tau _s\Vert A\Vert _{\infty }\), where \(\tau _s\) is a threshold that controls the accuracy of the factorization. With such a strategy, we are actually solving a perturbed system,where E is a diagonal matrix having nonzero entries equal to \(\tau _{s}\Vert A\Vert _{\infty }\) in the positions corresponding to pivots that were replaced. By applying [32, Theorem 9.3] to (4.24) we meet the condition (4.6) with \(\epsilon =\tau _s\), while condition (4.7) is met, since the triangular solves are standard. Therefore Theorem 4.1 applies with \(\epsilon =\tau _s\).

Finally, we can also derive convergence conditions for the case where BLR and static pivoting are combined. This amounts to using BLR approximations on the perturbed system (4.24), and so Theorem 4.1 applies with \(\epsilon =\tau _b+\tau _s\).

The potential gain in execution time of mixed-precision iterative refinement comes from the fact that the most time-consuming operation, the LU factorization, is carried in low-precision arithmetic and high precision is only used in iterative refinement steps, which involve low-complexity operations. For this gain to be effective, the cost of the refinement iterations must not exceed the time reduction resulting from running the factorization in low precision. This is very often the case. First, on current processors (such as the model used for our experiments), computations in single precision can be up to twice as fast as those in double precision. Additionally, operations performed in the refinement steps have a lower asymptotic complexity compared with the factorization. Nonetheless, in practice, the overall time reduction can vary significantly depending on a number of parameters. First, the ratio between the complexity of the factorization and that of the solution phase is less favorable on 2D problems than on 3D problems (see Section 2.1). Second, the single-precision factorization may be less than twice as fast as the double-precision one; this may happen, for example, on small problems where the overhead of symbolic operations in the factorization (data indexing, handling of data structures, etc.) is relatively high or, in a parallel setting, because the single-precision factorization is less scalable than the double-precision one due to the relatively lower weight of floating-point operations with respect to that of symbolic ones. It must also be noted that although the factorization essentially relies on efficient BLAS-3 operations, the operations done in the iterative refinement, in particular the LU solves, rely on memory-bound BLAS-2 operations and are thus less efficient. Finally, in the case of badly conditioned problems, iterative refinement may require numerous iterations to achieve convergence.

Figure 2 shows the execution time of variants encountered in Table 2 normalized with respect to that of DMUMPS for a selected subset of matrices from our test set; each bar also shows the time breakdown into LU factorization, LU solves and all the rest, which includes computing the residual and, for the GMRES-based variants, casting the factors, computing the Krylov basis, orthonormalizing it, and so on. The values on top of each bar are the number of LU solve operations; note that for GMRES-based variants, multiple LU solve operations are done in each refinement iteration.

In the first row of this figure, we report problems that behave well, in the sense that all the parameters mentioned above align in the direction that leads to a good overall time reduction. For these problems the single-precision factorization is roughly twice as fast as the double-precision one, the complexity of the solve is much lower than that of the factorization (three orders of magnitude in all cases, as reported in Table 1), and relatively few iterations are needed to achieve convergence. For all these problems, the complexity of the matrix–vector product is more than two orders of magnitude lower than that of the solve (see columns “NNZ” and “Slv.” of Table 1). As a result, the computation of the residual only accounts for a small portion of the total execution time—even for variants with \(u_r=\) q, for which it is carried out in slow quadruple precision arithmetic (which is not supported by our hardware). This is a very desirable property, since these variants greatly improve the forward error with only a modest overhead. The figure clearly shows, however, that despite their relatively low complexity, the operations in iterative refinement are relatively slow and, therefore, the gain is considerably reduced when many solves are necessary. This issue is exacerbated in the case of GMRES-IR variants, because the solves are carried out in double instead of single precision as for LU-IR variants (\(u_f=\) s).

In the second row of Figure 2, we report some cases where mixed-precision iterative refinement does not reduce execution time. Chevron4 is a relatively small 2D problem where the cost of the solve and the matrix–vector product relative to that of the factorization is high; as a result, even for a moderate number of refinement iterations, variant DMUMPS achieves the best execution time and all other variants are much slower. perf009ar is one where the single-precision factorization is only 1.6\(\times\) faster than the double-precision one and, additionally, it produces little fill-in (as shown by the small ratio Slv./NNZ in Table 1) and so the relative cost of computing the residual in quadruple precision is high. Finally, CarBody25M is badly conditioned and variants based on single-precision factorization either do not converge or require so many iterations that the execution time is higher than that of DMUMPS. It is, however, worth noting that on these particular matrices variants based on single precision factorization may be slower than DMUMPS but at a significantly reduced memory cost (as shown in Table 2).

One distinctive feature of the multifrontal method in comparison with left or right-looking ones is in the way it uses memory. In addition to the memory needed to store the factors, which grows monotonically throughout the factorization, the multifrontal method also needs a temporary workspace, which we refer to as active memory. As a result, the peak memory consumption achieved in the course of the factorization is generally higher than the memory needed to store the factors. It must also be noted that parallelism does not have any effect on the memory needed to store the factors but generally increases the size of the active memory; this is because more temporary data is generated at the same time. For a thorough discussion of the memory consumption in the multifrontal method, we refer the reader to the paper by Agullo et al. [1].

In the rest of this article, we refer to the difference between the peak memory consumption and the size of factors as active memory overhead. We assume that at the end of the factorization all the active memory is freed and only the factors are left in memory. It is only at this moment that the original problem matrix is cast to quadruple precision for computing the residual at each refinement iteration. Therefore, the active memory overhead and the memory required to store a quadruple precision version of the matrix do not accumulate. In our implementation, the GMRES-IR variants with \(u_p=u^2=\) d also require the factors to be cast to double precision, which we do upon completion of the factorization, when the active memory is freed. We also report the size of the Krylov basis in the GMRES solver; although in most of our experiments this is completely negligible, there might be cases (we will show one) where the number of GMRES iterations is sufficiently high to make the memory consumed by the Krylov basis relevant. Finally, we do not count the memory consumption of the solution, residual and correction vectors.

All these assumptions lead us to Figure 3 where we present the normalized memory consumption of certain LU-IR and GMRES-IR variants relative to that of variant DMUMPS for a selected subset of problems. We do not include variants using \(u_r=\) d, because they behave very similarly to variants with \(u_r=\) q. For each problem and variant the bar is split in two parts showing the memory consumption during and after the factorization, respectively.

In the first row, we report problems that behave well, which corresponds to the most common case as shown in Table 2. It shows, as expected, that LU-IR with single-precision factorization consumes half as much memory as DMUMPS, because the memory needed to store the problem matrix in quadruple precision does not exceed the active memory overhead. Thus, the highest memory consumption corresponds to the single-precision factorization peak. GMRES-based variant (\(u_p=\) d) casts the factors to double precision, which exceeds the peak of the single-precision factorization. Nonetheless, even if on top of this we have to add the memory needed to store the quadruple precision matrix, the overall consumption is lower than the double-precision factorization peak by a factor that can be almost up to 50% on this set, making the memory consumption of the GMRES-IR variant almost identical to that of the LU-IR one in a few cases (such as matrices nlpkkt80 and Serena). As for the LU-IR variant with \(u_f=\) d, it clearly does not bring any improvement with respect to DMUMPS but no loss either, because the memory for storing the matrix in quadruple precision is much lower than the active memory overhead.

In the second row of Figure 3, we report problems where the memory reduction is not as good, for four different reasons, the last two of which are exclusive to GMRES-IR.

In Table 3, we present the execution time and memory consumption of four iterative refinement variants using low rank BLR factorization for different values of the compression threshold \(\tau _b\). All variants provide a forward error on the solution equivalent to the one of DMUMPS. If \(\tau _b=\) “full-rank,” then the factorization is run without BLR, this is a standard factorization as in Section 5.3. It should be noted that in this case the double-precision factorization LU-IR variant is equivalent to DMUMPS, and we will refer to it as DMUMPS in the text. We denote by “—” the cases where convergence is not reached and, for each matrix, we highlight in bold the execution time and memory consumption that do not exceed by more than 10% the best execution time or memory consumption. We choose to work with the default BLR settings of MUMPS, in which case the data in the active memory is not compressed with low rank approximations. It should be noted, however, that MUMPS also offers the possibility to compress the data in the active memory; this would result in additional memory savings, but will badly affect the execution time, because it adds the compression overhead without reducing the operational complexity.

The experimental results of Table 3 are in good agreement with the theoretical convergence conditions of Theorem 4.1 developed in Section 4. We can clearly see how the robustness of the presented variants is related to both the condition number \(\kappa (A)\) of the matrix (specified for each matrix in Table 1) and the BLR threshold \(\tau _b\). Convergence is not achieved for excessively large values of the BLR threshold; the largest \(\tau _b\) value for which convergence is achieved depends on the matrix condition number and, in general, it is smaller for badly conditioned problems. In the case of GMRES-IR variants, which are more robust, the BLR threshold can be pushed to larger values without breaking convergence. Note that there is no situation where a variant does not converge when in theory it should (that is, when the convergence condition is met). However, as the theoretical convergence conditions in Theorem 4.1 can be pessimistic, there are several cases where convergence is achieved even though the theoretical convergence condition is not met.

The use of BLR with a good choice of compression threshold \(\tau _b\) generally results in substantial reductions of the LU-IR and GMRES-IR execution times. As the BLR threshold increases, the operational complexity and, consequently, the execution time of the factorization and solve operations decreases; conversely, the number of iterations increases up to the point where convergence may not be achieved anymore. The optimal BLR threshold value that delivers the lowest execution time obviously depends on the problem. It must be noted that even though the GMRES-IR variants achieve convergence for larger \(\tau _b\) values, this leads to an excessive number of iterations whose cost exceeds the improvement provided by BLR; as a result, these variants are slower than LU-IR ones (\(u_f=\) s but also \(u_f=\) d) in all cases. Consequently, the single-precision factorization LU-IR variant generally achieves the best execution time on this set of problems, similarly to what was observed in Section 5.3, with a few exceptions. On perf009ar the double-precision factorization LU-IR variant is the best due to the fact that similarly to the full rank case (see Section 5.3.1) the BLR factorization is less than twice as fast when single precision is used instead of double for the same \(\tau _b\) value; additionally, a substantial number of iterations is needed to achieve convergence. It is worth mentioning that on this matrix the GMRES-IR variant with \(u_g=u_p=\) d is faster than the single-precision factorization LU-IR variant (36.9 s versus 40.3 s) and consumes less memory than the double-precision factorization LU-IR variant (20.0 GB versus 37.1 GB). On CarBody25M, DMUMPS is the fastest variant as in the full rank case; this is due to the fact that, on this problem, BLR does not achieve a good reduction of the operational complexity and, therefore, of the execution time.

As for the storage, the use of BLR leads to a different outcome with respect to the case where a full-rank factorization is used (see Section 5.3) where the single-precision factorization LU-IR variant is the best. This is due to the combination of two factors. First, when BLR is used, the relative weight of the active memory is higher, because it corresponds to data that is not compressed due to the choice of parameters we have made; consequently, the memory consumption peak is often reached during the factorization rather than during the iterative refinement. Second, the memory consumption of the factorization decreases monotonically when the BLR threshold is increased. As a result of these two effects, the GMRES-IR variants achieve the lowest memory consumption on this set of problems, because they can preserve convergence for larger values of \(\tau _b\) than the LU-IR variants can. For example, on tminlet3M the GMRES-IR variant with \(u_g=u_p=\) d consumes almost \(15\%\) less memory than the LU-IR one with \(u_f=\) s (70.9 GB versus 82.4 GB), on thmgas the GMRES-IR variant with \(u_g=u_p=\) d consumes almost \(30\%\) less memory than variant LU-IR with \(u_f=\) s (43.7 GB versus 61.4 GB), and on matrix 24 the GMRES-IR variant with \(u_g=u_p=\) d consumes more than \(35\%\) less memory than variant LU-IR with \(u_f=\) d (41.8 GB versus 64.8 GB). It is worth pointing out that the value of \(\tau _b\) for which GMRES-IR achieves the lowest possible memory consumption is not always the largest value for which convergence is still possible. This is because for a large number of iterations the memory needed to store the Krylov basis may exceed the savings obtained with BLR. This problem can be overcome or mitigated by choosing an appropriate value for the \(\tau _g\) threshold or, similarly, using a restarted GMRES method; we leave this analysis for future work.

We finally compare the two GMRES-IR variants \(u_g=u_p=\) d and \(u_g=u_p=\) s. When \(u_g=u_p=\) s, GMRES-IR avoids the cast of the LU factors from single to double precision, and thus reduces memory consumption compared with \(u_g=u_p=\) d. However, as explained above, the relative weight of the factors with respect to the active memory is smaller as \(\tau _b\) increases, and so the reduction achieved by GMRES-IR with \(u_g=u_p=\) s grows smaller until the point where both variants achieve a similar memory consumption. On our matrix set, for the values of \(\tau _b\) where the LU-IR with \(u_f=\) s does not converge, GMRES-IR with \(u_g=u_p=\) s does not achieve significant memory reductions compared with GMRES-IR with \(u_g=u_p=\) d (at best 7% on perf009ar, 18.6 GB versus 20.0 GB).

We now turn our attention to the use of static pivoting. We report in Table 4 the execution time and memory consumption of three iterative refinement variants for different values of the static pivoting threshold \(\tau _s\). All variants are stopped when they reach a forward error on the solution equivalent to the one of DMUMPS. If \(\tau _s=\) “partial,” then the factorization is run in standard MUMPS threshold partial pivoting. It should be noted that in this case the double-precision factorization LU-IR variant is equivalent to DMUMPS.

Once again the observed convergence behaviors are in good agreement with Theorem 4.1 as explained below. In the case of static pivoting, the execution time of the factorization does not depend on \(\tau _s\); to minimize the overall solution time, the goal is therefore to achieve the fastest possible convergence. This is a complex issue: A smaller perturbation \(\tau _s\) does not always mean a faster convergence, because the value of \(\tau _s\) also directly impacts the growth factor \(\rho _n\). Thus, there is an optimal value of \(\tau _s\), which is clearly problem dependent, that leads to the fastest convergence by balancing the \(u_f\rho _n\) and \(\tau _s\) terms in the convergence condition. To confirm this, Table 4 reports \(\underline{\rho _n}\), a lower bound on the true growth factor that can be used as a cheap but rough indicator of the behavior of \(\rho _n\) (the true \(\rho _n\) would be extremely expensive to compute for such large matrices). There is a clear trend of \(\underline{\rho _n}\) decreasing as \(\tau _s\) increases, which explains, for example, why on lfm_aug5M convergence is achieved for large \(\tau _s\). For many matrices in our set, such as tminlet3M in Table 4, static pivoting slightly accelerates the factorization without excessively deteriorating the convergence, and so allows for modest time gains overall. However, for some matrices such as lfm_aug5M, static pivoting requires many iterations and can lead to significant slowdowns compared with partial pivoting. It is, however, interesting to note that, on lfm_aug5M, if the use of partial pivoting is impossible (for instance, because the available solver does not support it), then the GMRES-IR variant provides the best overall execution time.

Finally, in Table 5, we present the execution time and memory consumption of three iterative refinement variants (the same as in Section 5.4.2) for different values of the BLR compression threshold \(\tau _b\) and a fixed value of the static pivoting perturbation \(\tau _s\). All variants are stopped when they reach a forward error on the solution equivalent to the one of DMUMPS. If \(\tau _s=\) partial, then the factorization is run in standard MUMPS threshold partial pivoting and the results are then equivalent to the BLR results of Section 5.4.1.

Theorem 4.1 applied to the case where BLR and static pivoting are used together states that the convergence conditions should be affected by the largest perturbations \(\max (\tau _s, \tau _b)\) and the term \(\rho _nu_f\), which depends on the growth factor. Our experiments confirm this: Values of \(\tau _s\) or \(\tau _b\) for which a given variant was not converging with BLR or static pivoting alone still do not converge when they are combined, and, conversely, variants that were converging for BLR and static pivoting alone still converge when these two approximations are used together. lfm_aug5M with \(\tau _s=10^{-6}\) illustrates an interesting point of the error bound \(\max (\tau _s,\tau _b)+\rho _nu_f\); convergence is only achieved for the variant that uses a double-precision factorization (\(u_f={\rm\small D}\)), even for values of \(\tau _b\) that are much larger than the unit roundoff of single precision. This shows that the rule of thumb that the factorizaton precision should be chosen as low as possible as long as \(u_f \le \tau _b\) is not true in presence of large element growth, since a smaller value of \(u_f\) can be beneficial to absorb a particularly large \(\rho _n\).

While the reductions in execution time obtained by using static pivoting instead of partial pivoting were modest for the full-rank factorization, they are larger for the BLR factorization. Taking tminlet3M as an example, in full-rank the single-precision factorization LU-IR variant is only 1.12\(\times\) (136 s/121 s) faster after the activation of the static pivoting (see Table 4), whereas in BLR it is 1.26\(\times\) (88 s/70 s) faster (see Table 3). These better reductions are explained by the fact that in the BLR factorization, static pivoting also allows the panel reduction operation to be processed with low-rank operations [7], which leads to a reduction of flops and thus a faster execution time.