DOI: 10.1515/auom-2015-0044 An. Şt. Univ. Ovidius Constanţa Vol. 23(3),2015, 17–27 Preconditioners based on the ISM factorization∗ Rafael Bru and Juana Cerdán and José Marı́n and José Mas and Miroslav Tůma Abstract In this paper we survey our work on preconditioners based on the Inverse Sherman-Morrison factorization. The most important theoretical results are also summarized and some numerical conclusions are provided. 1 Introduction The solution of linear systems Ax = b, (1) where A is a nonsingular large and sparse matrix, is usually computed applying an iterative method, typically a Krylov method, see [17]. To get convergence, or to improve it, a preconditioner can be used. A preconditioner is a matrix M that approximates the inverse of A and can be constructed in several ways, even implicitly. When a preconditioner is used, instead of system (1) one of the systems M Ax = M b AM y = b, M 1/2 AM 1/2 y=M (left preconditioning) x = My 1/2 b x=M (right preconditioning) 1/2 y (2) (symmetric preconditioning) Key Words: Preconditioned iterative methods, sparse matrices, incomplete decompositions, approximate inverses, Sherman-Morrison formula. 2010 Mathematics Subject Classification: Primary 65F08, 65F10; Secondary 65F50. Received: November, 2014. Revised: November, 2014. Accepted: February, 2015. ∗ This work was supported by the Spanish Ministerio de Economı́a y Competitividad under grant MTM2014-58159-P and by the project 13-06684S of the Grant Agency of the Czech Republic. 17 PRECONDITIONERS BASED ON THE ISM FACTORIZATION 18 is solved by an iterative method. There are many ways to compute a preconditioner. Two of them are of interest for us, the one that computes an approximation of A−1 , the so called approximate inverse preconditioners; the second is based on the computation of an incomplete factorization of A, typically the well known Incomplete LU (ILU) factorization or the Incomplete Cholesky (IC) factorization for the symmetric positive definite case. Let us focus now on approximate inverse preconditioners, (see [3] for a detailed discussion). These preconditioners are applied by matrix-vector products in each iteration of the Krylov method. Then they are highly parallel. Some methods to obtain approximate inverse preconditioners are based on minimizing kIn − AM kF , subject to some sparsity constraints on M . The use of the Frobenius norm allows to reduce the problem to solving independent linear least square problems, which can be done also in parallel. Examples of this class are the SPAI [14] and MR [13] preconditioners. Other approaches compute the preconditioner M as a product of triangular factors by an Abiorthogonalization process, as the AINV method, see [1] and [2]. On the contrary, other preconditioners try to approximate A computing an incomplete LU [16] or Cholesky factorization. They are implicit because instead of computing the matrix M in (2), what it is computed is a factorization of its inverse, resulting in a factorization of a matrix that is close to the coefficient matrix A. Then, in each iteration instead of a matrix vector product with the matrix M two linear systems with triangular matrices must be solved. This implies that these kind of preconditioners are not as parallel as the approximate inverse ones. However they have some advantages. As A is usually much more sparse than its inverse, this implicit preconditioners can also be much more sparse than the explicit ones to get an equivalent improvement on the convergence of the iterative method. To get sparse preconditioners it is necessary to annihilate entries during computation. A threshold is selected and entries are zeroed if their absolute values are less than this quantity, usually multiplied by some quantity as, for example, the norm of the corresponding row. That is, an entry is zeroed if it is small relative to some other relevant quantity. In addition a maximum number of entries can be allowed by row or column. This technique is the most usual when computing preconditioners of both types summarized previously. The Inverse Sherman-Morrison (ISM) decomposition introduced in [8] that can be used to derive preconditioners of both types: approximate inverse preconditioners and incomplete LU or Cholesky preconditioners. In the next section we will review the general framework of the factorization and some theoretical results about the decomposition. In section 3 we will review preconditioners and we will also mention some other results. Our goal is to present PRECONDITIONERS BASED ON THE ISM FACTORIZATION 19 them in the order they have been known. 2 The ISM factorization In this section we present the ISM factorization. The factorization was introduced in [8] and later developed in several works. During the process some changes in the original notation have been introduced to improve readability. In the sequel we will use the new notation. The ISM factorization gives a shifted factorization of the inverse of a matrix and, as stated later, it also contains some of the LU factors of the matrix and its inverse. The starting point is the next result. Theorem 1 ([8, Theorem 2.1 and Corollary 2.2]). Let A and A0 be two nonsingular matrices, and let {xk }nk=1 and {yk }nk=1 be two sets of vectors such that n X xk ykT . (3) A = A0 + k=1 ykT A−1 k−1 xk Suppose that rk = 1 + 6= 0 for k = 1, . . . , n, where Ak−1 = A0 + Pk−1 −1 −1 −1 1 T −1 T i=1 xi yi . Then Ak is nonsingular and Ak = Ak−1 − rk Ak−1 xk yk Ak−1 . Moreover the vectors k−1 X v T A−1 xk i 0 zk := xk − zi (4) r i i=1 and vk := yk − k−1 X i=1 ykT A−1 0 zi vi ri (5) are well defined and satisfy the relations A−1 k−1 xk T −1 yk Ak−1 = A−1 0 zk T −1 = v k A0 (6) and T −1 rk = 1 + ykT A−1 0 zk = 1 + vk A0 xk (7) for k = 1, . . . , n. In addition −1 T −1 A−1 = A−1 0 −A 0 ZDV A0 , where Z = [z1 , . . . , zn ], D = diag(r1 , . . . , rn ) and V = [v1 , . . . , vn ]. (8) PRECONDITIONERS BASED ON THE ISM FACTORIZATION 20 Note that conditions rk 6= 0 , k = 1, . . . , n allow the use of the ShermanMorrison theorem to prove nonsingularity of matrices Ak . But taking into account (6) and (7) the process can be rewritten without these matrices. Different choices of A0 , {xk } and {yk } vectors can lead to different fac−1 torizations of A−1 . The choice that produces what we call the Inverse 0 −A Sherman-Morrison factorization is A0 = sI, x k = ek , yk = (ak − sek )T , (9) where s is a positive parameter, ek and ek are the kth column and row respectively of the identity matrix, and ak is the kth row of A. In the sequel we will denote the rows and columns of a matrix using superindexes and subindexes respectively. With the choices (9), equations (4), (5) and (7) become zk = ek − k−1 X i=1 vk = yk − k−1 X i=1 vki zi , sri ak zi vi , sri (10) ak zk vkk y T zk = =1+ rk = 1 + k s s s Then, if Zs = [z1 , . . . , zn ], Ds = diag(r1 , . . . , rn ) and Vs = [v1 , . . . , vn ], the ISM factorization is given by s−1 I − A−1 = s−2 Zs Ds−1 VsT . (11) In [8, Lemma 3.1] the relation of the factors for different values of the parameter s are established, in particular it is proved that Zs actually does not depend on s, so the subindex will be removed in the sequel. Algorithm 1 computes the ISM factorization of a given matrix A. It is worth to note the relation between vectors zk and vk . The ISM factorization is strongly related with the LU factorization as was proved in [9]. Theorem 2. A square matrix A has LDU factorization, A = LDU , if and only if it has ISM factorization s−1 I − A−1 = s−2 ZDs−1 VsT . Then D = s−1 Ds , U = Z −1 , Vs = U T D − sL−T . PRECONDITIONERS BASED ON THE ISM FACTORIZATION 21 Algorithm 1 Computation of the exact ISM factorization (11). (1) for k = 1, . . . , n zk = ek vk = (ak − sek )T for i = 1, . . . , k − 1 zk = zk − vsrkii zi k vk = vk − asrzii vi end for rk = 1 + vkk /s end for (2) Return Z = [z1 , z2 , . . . , zn ], V = [v1 , v2 , . . . , vn ] and D = diag(r1 , r2 , . . . , rn ). 3 Preconditioners based on the ISM factorization In this section we review all preconditioners derived up to now from the ISM factorization. They are approximate inverse preconditioners and also ILU type preconditioners. We also mention the application of them to least squares problems. 3.1 Approximate inverse preconditioners The exact process described in Algorithm 1 can be used to compute a preconditioner if a dropping strategy of small entries is implemented. In [8] the dropping strategy was based on the size of entries and several results about the incomplete factorization are proved. The most important one is that the incomplete process can be applied without breakdown to M-matrices. The approximate inverse preconditioner (AISM) was compared with AINV, showing similar robustness and performance. A block extension of AISM (bAISM) was introduced in [12]. If A is block partitioned in such a way that diagonal blocks are square, the block ISM reads −1 −1 T −1 A−1 = A−1 V A0 , 0 − A0 U T (12) where all matrices have the same block partition than A, U is block upper triangular and T is block diagonal. Two interesting results were proved. If A0 is taken as the identity matrix, then by [12, Lemma 3.2.] the kth diagonal block of T is exactly the kth pivot in the block LU factorization of A. Then the incomplete process does not breakdown if A is an H-matrix, [12, Theorem 3.4]. A similar result, [12, Theorem 3.7], is proved if A0 is the diagonal block matrix formed with the PRECONDITIONERS BASED ON THE ISM FACTORIZATION 22 diagonal blocks of A. It is worth to mention that this second case reduces to the first one if a right block Jacobi scaling is previously done to A [12, Theorem 3.8]. 3.2 ILU type preconditioners based on ISM In [9] the relation between the ISM and LDU factorizations, that was previously stated for pivots, was extended in a fruitful way as stated in Theorem 2. From this result it follows that the ISM factorization contains the D (maybe scaled) and U factors of the LDU factorization of A, and the factors L−1 , U −1 of the LDU factorization of A−1 , the other factor D−1 is easily computed as sDs−1 . This situation is specially interesting when A is a symmetric positive definite matrix. Then all factors of the LDLT factorizations of A are available and it is possible to get a preconditioner computing an incomplete Cholesky factorization. Moreover, since all factors of the Cholesky factorization of A−1 are also computed, special dropping strategies can be used, namely those based on the results obtained by M. Bollhöfer and Y. Saad for LU factorizations, [6, 7]. This robust strategy drops entries by value, but the value is compared with the norm of a row of the inverse. That is, given a drop tolerance τ an entry of the incomplete factor L̂ is zeroed if |ˆljk |kek L̂−1 k ≤ τ. Also the dual dropping criteria can be used. If `ˆjk denotes the entry of jth row, kth column of the incomplete factor L̂−1 , it is zeroed if |`ˆjk |kek L̂k ≤ τ. There is a balance between entries of the matrices L̂ and L̂−1 computed by the preconditioner obtained implementing this dual dropping strategies to sparsify matrices in the ISM process given in Algorithm 1. So it is called Balanced Incomplete Factorization (BIF). This balance gives raise to robust preconditioners as numerical experiments presented in [9] confirm. In addition they also show that the BIF preconditioner has a very good behavior from the point of view of the selection of the dropping parameter. As it is well known as the value of the dropping parameter decreases the density of the preconditioner increases. Although theoretically the number of iterations should decrease, for some methods and matrices this is not true. Even, in some cases, the preconditioner stops being able to get convergence for some interval of the dropping parameter, and then the method converges again. This is not the case of BIF, the method is robust and the number of iterations varies smoothly with a decreasing trend, see Figures 5.1–5.3 in [9]. On the other hand we get 23 PRECONDITIONERS BASED ON THE ISM FACTORIZATION −sU −1 −sL−T LD UTD Vs Ṽs Figure 1: Factors of the LDU factorizations of A and A−1 included in Vs and Ṽs . Green parts are factors of A, and blue parts are factors of A−1 . good results with very sparse preconditioners. Let us mention at this point that all codes we use to compute preconditioners implement also a restriction in the number of nonzero entries of the matrix, which probably is the reason that explain why we can not get denser preconditioners. 3.3 ISM preconditioners for general nonsingular matrices: nBIF As mentioned above, when A is a nonsymmetric matrix from the ISM factorization one can retrieve all factors in A−1 = U −1 D−1 L−1 , and factors D and U in A = LDU . To get the missed factor U it is proposed in [10] to compute the ISM factorization of AT . To fix notation let s−1 I − A−T = s−2 Z̃ D̃s−1 ṼsT , the ISM factorization of the inverse of AT . Then Vs = U T D − sL−T and Ṽs = LD − sU −1 , as depicted in Figure 1. As a consequence Z and Z̃ are no longer needed. PRECONDITIONERS BASED ON THE ISM FACTORIZATION 24 Moreover, as stated in Lemma 3.1 of [10] vpk = s k−1 X ṽki ṽkp − vpi . dk d i=p+1 i for p < k, A dual expression allows computation of ṽpk , for p < k, and then the strict upper triangular parts of Vs and Ṽs are computed in a new way. Entries of factors of the inverse are used to compute entries of factors of the matrix and vice versa. In BIF they were used only to drop entries. Numerical results in [10] show the same that was already true for BIF. The method is robust and the dependence of the dropping tolerance is not critical. 3.4 BIF for least squares problems It is clear that BIF can be applied to the iterative solution of least squares problems min kb − Axk2 x when the coefficient matrix has full column rank (and hence has at least as many rows as columns). If this is the case the matrix AT A is symmetric and positive definite and a BIF preconditioner can be computed. Anyway some particular aspects of the least squares problems should be taken into account. Since AT A is much denser than A, to get a sparse preconditioner a lot of entries should be dropped, and some instabilities may appear. To prevent this possibility two additional strategies are used in [11] to improve robustness. The first one is to replace the way that the pivots are computed. Instead of to compute them as in (10), that is as sri , they are computed as ziT AT Azi . The second one follows the ideas of Tismenetsky in [18], by storing additional entries that are used only in computations, but finally discarded a better preconditioner is computed, see Algorithm 2.2 in [11]. We refer to this preconditioner as lsBIF. In [11] lsBIF was compared with IC, an incomplete Cholesky factorization preconditioner see [4] for implementation details; AINV; and Cholesky incomplete modified Gram Schmidt, CIMGS [19], see [15] for details of the implementation. All preconditioners were used to accelerate the conjugate gradient for least squares (CGLS) [5] iterative method. Results show the robustness of the method and also its stability with respect the choice of the threshold. PRECONDITIONERS BASED ON THE ISM FACTORIZATION 4 25 Conclusions The ISM factorization can be used to compute robust, sparse and efficient preconditioners, taking advantage of the fact that it computes factors of the LDU factorization of the coefficient matrix and also of its inverse. Numerical experiments exhibits good performance for all preconditioners developed: direct ones, as BIF for symmetric matrices, nBIF for general matrices and lsBIF for least squares problems with full column rank; as well as those of approximate inverse type (AISM and bAISM). The relation between the factorization and the LDU decomposition also means that the incomplete factorization can be computed breakdown free for H-matrices, and in the case of full column rank matrices for the normal equations. References [1] M. Benzi, C. D. Meyer, and M. Tůma. A sparse approximate inverse preconditioner for the conjugate gradient method. SIAM J. Sci. Comput., 17(5):1135–1149, 1996. [2] M. Benzi and M. Tůma. A sparse approximate inverse preconditioner for nonsymmetric linear systems. SIAM J. Sci. Comput., 19(3):968–994, 1998. [3] M. Benzi and M. Tůma. A comparative study of sparse approximate inverse preconditioners. Appl. Numer. Math., 30(2-3):305–340, 1999. [4] M. Benzi and M. Tůma. A robust preconditioner with low memory requirements for large sparse least squares problems. SIAM J. Sci. Comput., 25(2):499–512, 2003. [5] Å. Björck. Numerical methods for Least Squares Problems. Philadelphia, 1996. SIAM, [6] M. Bollhöfer. A robust and efficient ILU that incorporates the growth of the inverse triangular factors. SIAM J. Sci. Comput., 25(1):86–103, 2003. [7] M. Bollhöfer and Y. Saad. On the relations between ILU s and factored approximate inverses. SIAM J. Matrix Anal. Appl., 24(1):219–237, 2002. [8] R. Bru, J. Cerdán, J. Marı́n, and J. Mas. Preconditioning sparse nonsymmetric linear systems with the Sherman-Morrison formula. SIAM J. Sci. Comput., 25(2):701–715, 2003. PRECONDITIONERS BASED ON THE ISM FACTORIZATION 26 [9] R. Bru, J. Marı́n, J. Mas, and M. Tůma. Balanced incomplete factorization. SIAM J. Sci. Comput., 30(5):2302–2318, 2008. [10] R. Bru, J. Marı́n, J. Mas, and M. Tůma. Improved balanced incomplete factorization. SIAM J. Matrix Anal. Appl., 31(5):2431–2452, 2010. [11] R. Bru, J. Marı́n, J. Mas, and M. Tůma. Preconditioned iterative methods for solving linear least squares problems. SIAM J. Sci. Comput., 36(4):A2002–A2022, 2014. [12] J. Cerdán, T. Faraj, N. Malla, J. Marı́n, and J Mas. Block approximate inverse preconditioners for sparse nonsymmetric linear systems. ETNA, 37:23–40, 2010. [13] E. Chow and Y. Saad. Approximate inverse preconditioners via sparsesparse iterations. SIAM J. Sci. Comput., 19(3):995–1023, 1998. [14] M. J. Grote and T. Huckle. Parallel preconditioning with sparse approximate inverses. SIAM J. Sci. Comput., 18(3):838–853, 1997. [15] I. E. Kaporin. High quality preconditioning of a general symmetric positive definite matrix based on its U T U + U T R + RT U decomposition. Numer. Linear Algebra Appl., 5:483–509, 1998. [16] Y. Saad. ILUT: a dual threshold incomplete LU factorization. Numer. Linear Algebra Appl., 1(4):387–402, 1994. [17] Y. Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Co., Boston, 1996. [18] M. Tismenetsky. A new preconditioning technique for solving large sparse linear systems. Linear Algebra Appl., 154–156:331–353, 1991. [19] X. Wang, K. A. Gallivan, and R. Bramley. CIMGS: an incomplete orthogonal factorization preconditioner. SIAM J. Sci. Comput., 18(2):516–536, 1997. Rafael Bru, Instituto de Matemática Multidisciplinar, Universitat Politécnica de Valencia, Camı́ de Vera s/n, 46022 Valencia, Spain. Email: rbru@imm.upv.es Juana Cerdán, Instituto de Matemática Multidisciplinar, Universitat Politécnica de Valencia, Camı́ de Vera s/n, 46022 Valencia, Spain. Email: jcerdan@imm.upv.es PRECONDITIONERS BASED ON THE ISM FACTORIZATION José Marı́n, Instituto de Matemática Multidisciplinar, Universitat Politécnica de Valencia, Camı́ de Vera s/n, 46022 Valencia, Spain. Email: jmarinma@imm.upv.es José Mas, Instituto de Matemática Multidisciplinar, Universitat Politécnica de Valencia, Camı́ de Vera s/n, 46022 Valencia, Spain. Email: jmasm@imm.upv.es Miroslav Tůma, Institute of Computer Science, Academy of Sciences of the Czech Republic, Pod vodárenskouvěžı́ 2, 182 07 Prague 8, Czech Republic. Email: tuma@cs.cas.cz 27