Balanced incomplete factorization preconditioner with pivoting

In this work we study pivoting strategies for the preconditioner presented in Bru (SIAM J Sci Comput 30(5):2302–2318, 2008) which computes the LU factorization of a matrix A. This preconditioner is based on the Inverse Sherman Morrison (ISM) decomposition [Preconditioning sparse nonsymmetric linear systems with the Sherman–Morrison formula. Bru (SIAM J Sci Comput 25(2):701–715, 2003), that using recursion formulas derived from the Sherman-Morrison formula, obtains the direct and inverse LU factors of a matrix. We present a modification of the ISM decomposition that allows for pivoting, and so the computation of preconditioners for any nonsingular matrix. While the ISM algorithm at a given step computes only a new pair of vectors, the new pivoting algorithm in the k-th step also modifies all the remaining vectors from k+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k+1$$\end{document} to n. Thus, it can be seen as a right looking version of the ISM decomposition. The results of numerical experiments with ill-conditioned and highly indefinite matrices arising from different applications show the robustness of the new algorithm, since it is able to solve problems that are not possible to solve otherwise.


Introduction
This paper is concerned with the computation of robust preconditioners by using standard pivoting techniques for solving ill-conditioned sparse nonsingular linear systems of equations of the form  (1) such that it can be solved efficiently by means of iterative methods still remains one of the most active research areas in numerical linear algebra. Preconditioned Krylov methods have been traditionally linked to the solution of large and sparse linear systems due to the relative small amount of memory and computation time needed to obtain an approximate solution compared with direct methods. There exist different techniques that can be used successfully to compute preconditioners, as incomplete LU factorizations, approximate inverses, algebraic methods, etc. (see [2] and the references therein). But in recent years they have also been employed in the context of mixed precision techniques for solving dense linear systems. The accuracy of an initial solution obtained with an LU factorization computed in single precision is improved by iterative refinement using the LU factorization as preconditioner [10,15]. Ill-conditioned nonsingular linear systems arise in many areas of scientific and engineering applications, and computing numerically stable LU factorizations for these linear systems becomes a challenge [1,13]. Pivoting has been originally used to compute good LU factorizations for such problems [14], but there have been also some work for incomplete factorizations [6,18,21]. Also MATLAB has incorporated this possibility into his function ilu, that computes the Incomplete LU factorization of a matrix, [17].
There are different pivoting techniques being partial, complete and rook pivoting the more important ones [14,19]. Basically, at a given step of Gaussian elimination pivoting looks for an element sufficiently large in magnitude in the remaining submatrix, the Schur complement, to use it as the next pivot. These techniques involve row and possibly column permutations of the matrix that supposes a computational overhead. In this sense, partial pivoting is the cheapest pivoting technique, since it looks only in the first column of the Schur complement. Close behind is rook pivoting [20], it selects a pivot with maximum absolute value in its row and column, moving first to the biggest entry in magnitude in the first column, then it moves in the corresponding row, and then again in the column, and so on until the requirement is fulfilled. Finally, complete pivoting is the most expensive one, but guarantees the largest pivot at any stage however because the pivot is the entry of biggest magnitude in all the Schur complement.
In this work we study pivoting techniques for the balanced incomplete factorization preconditioner, BIF. BIF preconditioning is based on the incomplete Sherman-Morrison decomposition, ISM. The ISM decomposition uses recursion formulas derived from the Sherman-Morrison formula and was introduced in [7] as a method for computing approximate inverse preconditioners. In [8] the authors show that, applying the ISM algorithm to a symmetric and positive definite matrix A it is possible to compute an incomplete Cholesky factorization, later in [9] they showed that applying the ISM algorithm to A and A T , it is possible to compute an incomplete LDU factorization. Moreover, in both cases the inverse factors are also available and they influence the computation of the Cholesky or LDU factorization, and vice versa. In addition, the availability of the direct and inverse factors is exploited to implement norm based dropping rules [5]. The numerical results show that BIF is a robust algorithm comparable to other techniques as ILU(τ ) [3], ILUT (Threshold Incomplete LU) [22] and RIF (Robust Incomplete Factorization) [4]. Nevertheless, as mentioned above, computing stable (incomplete) factorizations for general ill-conditioned problems still require the application of pivoting techniques, except for some kind of matrices that can be solved with high accuracy, [16]. In this paper we show that with a slight modification of the ISM recursion formulas it is possible to incorporate pivoting to BIF, and obtain an algorithm which can be efficiently implemented and with similar efficiency that the well known incomplete LU with pivoting (ILUP). So, this study completes our previous work.
The paper is organized as follows. In Sect. 2 an overview of the ISM decomposition is presented. In Sect. 3 we introduce and analyze the right looking ISM decomposition. It is shown that the Schur complement computed with Gaussian elimination is available at each step of the modified algorithm and therefore, it is possible to incorporate any standard pivoting technique. In Sect. 4 the BIF algorithm with pivoting is presented and the results for several ill-conditioned matrices are reported. Our experiments show that BIF with pivoting is able to solve such a challenging problems and it is comparable to ILUT with partial pivoting. Finally, the main conclusions are presented in Sect. 5.

The ISM decomposition
The ISM decomposition was introduced in [7] as an algorithm to compute approximate inverse preconditioners since it obtains a factorization of the (shifted) inverse matrix of A, as where s > 0 is a given scalar and the columns of the matrices Z and V s are computed using the recursion formulas for k = 1, 2, . . . , n. In (3) the vector e k (e k ) denotes the k−th column (row) of the identity matrix, y k = (a k − se k ) T where a k denotes the k-th row of A, and are the entries of the diagonal matrix D s . It was proved in [8] that for symmetric matrices the factorization A = L DL T and the decomposition (2) satisfy The algorithm to get the decomposition of A uses explicitly the computed factors of A −1 , that is, A −1 is implicitly factorized at the same time. In [9] it is proved the following result Theorem 1 (Theorem 2.1 of [9]) Let A = L DU be the LDU decomposition of A, and let (2). Then Observe that L does not appear in (5). Therefore, to get the LU factorization for general matrices it is necessary to compute also the ISM decomposition of A T that gives as result where we have denoted with tilde the factors of the ISM decomposition of A T .
It is well known that a nonsingular matrix A has an LU factorization if there exists a lower unit triangular matrix L and an upper triangular matrix U , such that A = LU . The LDU factorization is obtained from the LU factorization by taking D as the diagonal matrix whose entries are the diagonal entries of U , and applying its inverse to U as D −1 U . Both factorizations are closely related with Gaussian elimination. Note that not all the nonsingular matrices have LU factorization since a zero pivot can be found during the Gaussian elimination process. However it is always possible to permute some rows, and maybe some columns of the matrix in such a way that the permuted matrix P AQ has LU factorization. Here P and Q are permutation matrices acting on rows and columns of A, respectively.
The idea is that it is possible to find permutation matrices P and Q such that at the k-th step of the Gaussian elimination process one obtains the matrix where the Schur complement A 12 is nonsingular and its first diagonal element is nonzero. Then, the permuted matrix P AQ is factorized as Here, A 11 , A 12 , A 21 and A 22 represent the submatrices of the reordered matrix P AQ, and the size of A 11 is k × k.
Note that in practice, the permutation matrices P and Q are not known in advance and therefore LU factorization algorithms determine which rows and columns must be interchanged during the elimination process. In the next section we show that it is possible to obtain the factorization (6) from the ISM decomposition by modifying its recursion formulas.

Right looking ISM algorithm
To implement pivoting in the ISM decomposition it is necessary to know the Schur complement of the LU factorization. To accomplish that, the vectors z k and v k must be computed in a different way. Instead of computing only one pair of vectors in the k-th step of the algorithm according to equations (3), the modification consist in updating also the remaining vectors, from k + 1 to n. That is, the right part of the matrices Z and V are updated in each step. The following MATLAB code implements the new right looking version of ISM. Next, we will show that the Schur complement S (k) is available from the matrix V . As usual, we denote by u k j the entry in the row k and column j of the matrix U , and by l ik the entry in the row i and column k of the matrix L. Observe that after the k-th step the first k columns of Z and V are computed. From A = L DU and Theorem 1, AZ = L D, and after the k-th step and the j- From (7) and considering that z k has zero entries bellow the row k, an important equality for the proof of our main result is We denote by V

Theorem 2 If A is a nonsingular matrix, then at the k-th step of the right looking ISM
Proof We are going to prove equation (10) by induction on the steps.
For k = 1 let us consider the element a (1) i j of the Schur complement S (1) , that is entries with i, j > 1. It is well known In the right looking ISM the ( j, i) entry of the matrix V (1) where we have used Eqs. (8) and (9). Working in the same way when i = j, we have v (1) Then V Let us prove the equality for the k-th step. Consider the entries of the Schur complement S (k) , that is a where we have used Eqs. (8) and (9). Working in the same way, Since we need to compute also the factorization of A T , observe that. Also note that the pivoting strategy should be decided by looking into the Schur complement contained in V s , or that inṼ s , but not both. In contrast, for complete pivoting it is clear that V s orṼ s produce the same pivot in exact arithmetic so any of them or both may be used.

Numerical experiments
In this section we report the results of some numerical experiments with a set of matrices from The SuiteSparse Matrix Collection [11] and the Harwell-Boeing collection [12]. The matrices are listed in Table 1 where their size, number of nonzeros, condition number and application are indicated. They correspond to very ill-conditioned and highly indefinite problems for which Gaussian elimination without pivoting fails to compute good quality L and U factors, so the same is expected to be the case for incomplete LU factorizations (see [5,10]). Partial, rook and complete pivoting techniques have been tested. The experiments have been implemented and run in MATLAB R2022a. As iterative solvers the MATLAB implementation of full GMRES [23] and BiCGStab [24] were used. The right hand side vector was computed such that the solution was the vector of all ones, and the initial guess was the zero vector. The iterations were stopped when the initial residual was reduced by 8 orders of magnitude with a maximum number of 1, 000 iterations. To compare the results obtained with BIF the problems were also solved with the MATLAB's incomplete LU preconditioner with partial pivoting, ILUTP. The implementation of the BIF preconditioner is based on the algorithm described in [9] but with the right looking modification described in Sect. 3. In [9] the authors show that in the ISM factorization the computation of the direct and inverse LU factors is interleaved and they mutually influence each other. These characteristics allows for the use of advanced dropping rules, see [5]. We will not discuss in detail these rules but we recall that their application requires the estimation of the norm of the columns of the LU factors and their inverses. Since approximations of these factors are explicitely available the application of this kind of dropping rules is straightforward.
For simplicity, all the experiments have been done with the the parameter s of the ISM decomposition equal to one. The algorithm is implemented such that the ISM decompositions of A and A T are computed at the same time. Therefore, accessing A and A T simultaneously is needed. The pivot is chosen from the Schur complement contained in V s rather thanṼ s . We note that for complete pivoting the same pivot could be obtained working either with V s orṼ s but we choose working with V s for simplicity. Thus, if in the k-th step of the algorithm the pivot strategy determines that the new pivot is in column p and row q, since V s stores the transpose of the Schur complement, then the columns p and k in V s and A T must be interchanged. Observe that inṼ s ,Z and A the rows p and k must be pivoted instead. By the same reason the rows q and k must be interchanged in V s , Z and A T , and the corresponding columns inṼ s and A. Also, other vectors whose elements depend on the column and row ordering, for instance vectors storing the norms of the columns needed for the dropping rule, must be reordered accordingly. Algorithm 2 sketches the pivoted version of the BIF algorithm described. The algorithm first determines the pivot position according to the pivoting strategy applied, complete, rook or partial. Then, permutation of the matrices and vectors are performed. After that, norms of the columns of the LU factors and their inverses are updated and the dropping rule is applied. The recursion formula of the factorization is applied. Finally, after n steps the LU factors are extracted from matrices V andṼ .
In Tables 2 and 3 the pivoting strategy is indicated with C, P and R for the complete, partial and rook pivoting strategies, respectively. Densit y is the ratio between the number of nonzeros of the preconditioner and the number of nonzeros of the matrix, that is nnz(L)+nnz(U ) .
Column iter shows the number of iterations of the solver and droptol is the tolerance used to drop elements in BIFP and ILUTP. The other columns are self explanatory. To reduce the numbers in the tables, a blank space means that the value is the same appearing in previous  rows. For instance, in Table 2 the droptol value for BIFP was always 10 −6 and therefore it appears only in the first row. The same holds for the preconditioner densities which are the same for GMRES and BiCSTAB and therefore only indicated once. Next, we will comment on the results. We note that the matrices tested can not be solved without pivoting with both BIFP and ILUTP preconditioners. Thus, pivoting is an essential tool to gain robustness for these factorizations. Starting with the University of Florida test matrices, Table 2, we observe for the adder group that there are not big differences between the different pivoting strategies for BIFP. Density is small, except for adder_dcop_06. The same can be said for the number of iterations spent by both iterative solvers. For the rest of matrices one can see that BIFP with complete pivoting computes sparser preconditioners than partial and rook pivoting. The iteration count does not present remarkable differences except for the oscil_dcop_01 matrix for which GMRES with partial pivoting, although with larger nonzero density, doubles the number of iterations.
For the Harwell-Boeing matrices reported in Table 3 the first thing to observe is that the preconditioners are quite dense, with the exception of partial pivoting for the matrix orani678. Complete pivoting still produces less fill-in in the preconditioner and the number of iterations is similar to the other pivoting strategies. Note however that partial pivoting performs extraordinary well for the matrix orani678 since it is able to converge in the same number of iterations with a very sparse preconditioner.
Finally, comparing the performance of BIFP with ILUTP we did not observed significant differences, specially with the preconditioned BiCGStab method. We recall that ILUTP uses partial pivoting and we observe that BIFP with this pivoting strategy performed closely in most cases.

Conclusions
In this paper we have presented an improved version of the BIF preconditioner that incorporates pivoting. The algorithm relies on a modification of the recursion formulas such that the Schur complement of standard Gaussian elimination is available at each step of the factorization. Thus, the application of different pivoting techniques, as for instance partial, rook and complete pivoting, can be done in a straightforward manner. Incorporating pivoting turns out to be an important step in order to achieve our initial goal of obtaining a more robust preconditioner since it is able to solve very ill-conditioned and indefinite problems that it may not be possible to solve in other way. The results of the numerical experiments with several matrices arising in different applications confirm that BIF with pivoting is a robust algorithm. Partial, rook and complete pivoting has been tested. Although complete pivoting very often produces sparser preconditioners with a competitive iteration count, rook and partial pivoting perform also quite well. Taking into account that partial and rook are less expensive from a computational point of view since they need less comparisons in order to determine the pivot, these two techniques may be preferable as default. Also, a comparison with ILU with partial pivoting (ILUP) has been done and one can see that the results between them are fairly close. As a final note on future work, since with the ISM decomposition one can also compute incomplete approximate inverse preconditioners, it can be worth to explore its application for ill-conditioned problems, as it is done for instance in [15].