A generalization of the optimal diagonal approximate inverse preconditioner

A generalization of the optimal diagonal approximate inverse preconditionerI Luis Gonz´alez∗, Antonio Su´arez, Eduardo Rodr´ıguez Department of Mathematics, University of Las Palmas de Gran Canaria, 35017 Las Palmas de Gran Canaria, Spain Abstract The classical optimal (in the Frobenius sense) diagonal preconditioner for large sparse linear systems Ax =bis generalized and improved. The new proposed approximate inverse preconditioner Nis based on the minimization of the Frobenius norm of the residual matrix AM −I, where Mruns over a certain linear subspace of n×nreal matrices, defined by a prescribed sparsity pattern. The number of nonzero entries of the n×npreconditioning matrix Nis less than or equal to 2n, and nof them are selected as the optimal positions in each of the ncolumns of matrix N. All theoretical results are justified in detail. In particular, the comparison between the proposed preconditioner Nand the optimal diagonal one is theoretically analyzed. Finally, numerical experiments reported confirm the theory and illustrate that our generalization of the optimal diagonal preconditioner improves (in general) its efficiency, when they do not coincide. Keywords: Approximate inverse preconditioner, Frobenius norm minimization, Diagonal preconditioner 1. Introduction The discretization of many different PDEs (modeling physical problems) by any adequate numerical method (finite differences, finite elements, finite volumes, meshless, etc.), generally leads to a large linear system IShort running title: A generalization of the optimal diagonal preconditioner ∗Corresponding author. Email address: [email protected] (Luis Gonz´alez) Preprint submitted to Computers &Mathematics with Applications September 10, 2013 Ax =b, A ∈Rn×n, x, b ∈Rn×1(1.1) in which the matrix Ais nonsingular and sparse. The solution of these linear systems is usually performed by iterative methods based on Krylov subspaces (see, e.g., [1, 2, 3, 4]). To improve the convergence of these Krylov methods, system (1.1) can be preconditioned with an adequate preconditioning matrix N, transforming it into any of the equivalent problems NAx =Nb, (1.2) ANy =b, x =Ny, (1.3) that is, the left and right preconditioned systems, respectively. In this paper, we address only the case of right-sided preconditioners (1.3), but analogous results can be obtained for the left-sided preconditioners (1.2). The study of preconditioning strategies for large linear systems is at present one of the most relevant research areas in Numerical Linear Algebra. In [5], we can find a very complete survey about this question. The preconditioning of system (1.1) is performed in order to obtain a preconditioned matrix AN as close as possible to the identity in some sense, and the preconditioner Nis called an approximate inverse of A. The different strategies to construct approximate inverse preconditioners can be grouped into three categories [6]: approximate inverse methods based on Frobenius norm minimization, factorized sparse approximate inverses (see, e.g., [7, 8, 9] and the references therein), and preconditioning methods consisting of an incomplete factorization followed by an approximate inversion of the incomplete factors. The idea of using Frobenius norm minimization for preconditioning purposes was first described in [10], and other early works can be found in [11, 12, 13]. Some posterior approaches in this sense can be found, for instance, in [14, 15, 16, 17, 18, 19] and in the references therein. In some cases, the Frobenius norm based preconditioners are parametrized by prescribed sparsity patterns. Otherwise, among the Frobenius norm minimization preconditioners not extracted from sparse matrix subspaces, let us mention here the preconditioners for structured matrices obtained by orthogonal projections onto unitary matrix algebras (like, for instance, circulant preconditioners for Toeplitz matrices); see, e.g., [20, 21, 22] and the references therein. 2 In [23, 24], the search of Frobenius norm based approximate inverses with a prescribed sparsity pattern is generalized by considering a more general case of linear parametrization where preconditioners belong to an arbitrary matrix subspace Sof Rn×n. This procedure leads to a natural generalization of the classical Moore-Penrose inverse, the so-called S-Moore-Penrose inverse introduced in [25]. The closeness of the preconditioned matrix AN to the identity may be measured by using a suitable matrix norm like, for instance, the Frobenius norm k·kF. In this way, the problem of obtaining the best preconditioner N (with respect to the Frobenius norm) of system (1.1) in the subspace Sof Rn×nis reduced to the minimization problem min M∈S kAM −IkF=kAN −IkF(1.4) and the solution Nto problem (1.4) will be referred to as the “optimal” preconditioner of system (1.1) over the subspace S. It is important to highlight that, throughout this paper, the term “optimal” means that the approximate inverse Nis the matrix that minimizes the Frobenius norm on AN −Iover a certain subspace Sof Rn×n, but the preconditioner Nis not necessarily optimal in any other sense of the word. Let us briefly describe the basic idea of this work. Our starting point is the well-known optimal diagonal preconditioner; see, e.g., [4]. This is exactly the solution Dto problem (1.4) for the subspace of all n×ndiagonal matrices, and it is often used as a simple preconditioner for sparse linear systems. Sometimes, the preconditioner Dis efficient and it leads to fast convergence. For instance, this is usually the case when matrix Ais symmetric positive definite [2]. However, in other cases, the diagonal preconditioner D is not effective enough for convergence. Then, we improve Din the following natural way. Since, obviously, the diagonal matrix Dhas one and only one nonzero element per column, this suggests the idea of considering the best approximate inverse (in the Frobenius sense) of matrix Aamong all the n×nmatrices that have exactly one nonzero element per column. Call each of such nonzero elements the optimal row position or entry for its corresponding column. Then, our proposed preconditioner Nwill contain the ndiagonal entries and those optimal entries per column which do not coincide with the diagonal ones. Finally, Nwill be exactly the solution to problem (1.4) for the subspace S ⊂ Rn×n, defined by the above described sparsity pattern. 3 Obviously, the so defined approximate inverse Nof matrix Ageneralizes D, and it has at least nnonzero entries (the diagonal ones), and at most 2n nonzero entries. Moreover, the preconditioning matrix Nhas another advantage, compared with the classical diagonal approximate inverse D. Namely, the reiteration of the preconditioning technique with the optimal diagonal approximate inverse makes no sense. On the contrary, when using our new preconditioner N, the well-known multistep preconditioning strategy (see, e.g., [15, 26]) not only makes sense, but as we shall prove, each step of this reiterated preconditioning strategy strictly reduces the Frobenius norm of the residual matrix, whenever the preconditioner Nobtained in the previous step is not diagonal. We propose a simple, natural generalization of the optimal diagonal preconditioner, which improves it (in the sense of Eq. (1.4)). Theoretical results will be justified in detail and illustrated with some numerical experiments. In addition, the proposed preconditioner is also compared with the AINV approximate inverse preconditioner [6]. This paper has been organized as follows. In Section 2, we recall explicit expressions for both the solution Nto problem (1.4) and its corresponding minimum Frobenius norm kAN −IkF, valid for any matrix subspace S ⊂ Rn×n. Next, in Section 3, we derive explicit expressions for the proposed preconditioner Nand for kAN −IkF. Numerical experiments are presented in Section 4. Finally, Section 5 closes the paper with some concluding remarks. 2. A preliminary lemma Now, we present a preliminary lemma required to make this paper selfcontained. Taking advantage of the prehilbertian character of the matrix Frobenius norm, the solution Nto problem (1.4) can be directly obtained using the orthogonal projection theorem. Here and in the following, orthogonality is with respect to the Frobenius inner product h·,·iF. More precisely, the matrix product AN is the orthogonal projection of the identity onto the subspace AS. Consequently, an explicit formula for matrix Ncan be obtained by expressing the orthogonal projection AN of the identity matrix onto the subspace ASby its expansion with respect to an orthonormal basis of AS [23]. This is the idea of the following lemma. 4 Lemma 2.1. Let A∈Rn×nbe nonsingular. Let Sbe a linear subspace of Rn×nof dimension d, and {M1, ..., Md}a basis of Ssuch that {AM1, ..., AMd} is an orthogonal basis of AS. Then, the solution to problem (1.4) is N= d X i=1 tr (AMi) kAMik2 F Mi,(2.1) and the minimum Frobenius norm is kAN −Ik2 F=n− d X i=1 [tr (AMi)]2 kAMik2 F .(2.2) Remark 2.1. If we have a basis {Mi}d i=1 of subspace Ssuch that the corresponding basis {AMi}d i=1 of subspace ASis not orthogonal, then we only need to use the Gram-Schmidt orthogonalization procedure to obtain an orthogonal basis of AS, in order to apply Lemma 2.1. This procedure has been formalized in [23], obtaining several explicit expressions for both the optimal preconditioner Ndefined by (1.4) and kAN −IkF, that have been applied to the sparse preconditioning of large linear systems arising from real-world cases. For different spectral properties of matrix AN, and for the theoretical effectiveness analysis of the optimal approximate inverse preconditioners N defined by Eq. (1.4), we refer the reader to [24, 25]. 3. The proposed approximate inverse preconditioner In this section, the proposed preconditioner Nof system (1.1) is introduced. First, we need to give a definition and to set some notations. Comparing two different approximate inverses for the same matrix A, as stated by the following definition, is an essential point for preconditioning purposes. Definition 3.1. Let A, N, N0∈Rn×nand suppose that Ais nonsingular. Then, we say that Nis better approximate inverse of Athan N0, or that N improves N0as approximate inverse of Aif and only if kAN −IkF<kAN0−IkF. 5 Throughout this paper, the subspace of all n×ndiagonal matrices is denoted by Dn. From now on, Mi,j denotes the n×nmatrix whose only nonzero term is mij = 1, eidenotes the ith column of the identity matrix (i.e., Aeiis the ith column of A), and the symbols k·k2and h·,·i2stand for the usual Euclidean vector norm and inner product, respectively. Remark 3.1. Note that since Mi,j =eieT j, then the only non-null column of matrix AMi,j is its jth one, which coincides with the ith column Aeiof matrix A. Consequently, we have tr (AMi,j) = tr AeieT j=aji,kAMi,jk2 F= AeieT j 2 F=kAeik2 2.(3.1) Moreover, hAMi,j, AMi0,jiF=AeieT j, Aei0eT jF=hAei, Aei0i2,(3.2) hAMi,j, AMi0,j0iF=AeieT j, Aei0eT j0F= 0 for all j6=j0,(3.3) so that any system of matrices {AMi,j}n j=1 is orthogonal with respect to the Frobenius inner product. When all diagonal entries of matrix Aare not null, the preconditioner diag a−1 11 , a−1 22 , . . . , a−1 nn is often used as an approximate inverse preconditioner of system (1.1); see, e.g., [2]. However, in general, this is not the optimal choice (in the sense of Eq. (1.4)) among the diagonal approximate inverses. Indeed, as it is well-known, the best diagonal preconditioner Dof system (1.1), that is, the solution to problem (1.4) for the subspace Dnis given by (see, e.g., [27]) D= n X j=1 ajj kAejk2 2 Mjj =diag a11 kAe1k2 2 ,··· ,ann kAenk2 2!,(3.4) while the corresponding minimum Frobenius norm is given by  AD −I 2 F=n− n X j=1 a2 jj kAejk2 2 .(3.5) Obviously, the optimal diagonal approximate inverse (3.4) of matrix A, has exactly one nonzero element per column. As mentioned in Section 1, 6 this suggests the idea of considering the best approximate inverse of matrix Aamong all the n×nmatrices that have exactly one nonzero element per column, the so-called optimal row position or entry per column. Suppose that the nnonzero optimal entries are placed at positions (i1,1) ,(i2,2) ,...,(in, n), i.e., for each column j= 1,2, . . . , n, the optimal entry for preconditioning the linear system (1.1), by using Eq. (1.4), is placed at the ijth row. Then, our new preconditioner Nis defined as follows. (i) If ij=j, the best entry in the jth column is the diagonal one. We select this entry, and no other entries are added to (j, j) in column j. (ii) If ij6=j, the best entry in the jth column is not the diagonal one. We select the diagonal entry (j, j) and, besides, the optimal entry (ij, j) is added to column j. (iii) Finally, our preconditioner Nis defined as the solution to problem (1.4) for the subspace Sof Rn×nwhose only nonzero entries are the ones defined by steps (i) and (ii), i.e., S=Sn:= span {Mj,j}n j=1 ∪Mij,j |ij6=jn j=1. In this way, the number of nonzero entries of each column j= 1,2, . . . , n of matrix Nis either 1 (if ij=j) or 2 (if ij6=j). Hence, the total number of nonzero entries of the n×npreconditioning matrix Nis at least nand at most 2n. For instance, let n= 4 and suppose that for a certain coefficient matrix A∈R4×4, the optimal positions per column are (i1,1) = (1,1) ,(i2,2) = (4,2) ,(i3,3) = (3,3) ,(i4,4) = (1,4) . Then, the sparsity patterns of the preconditioning matrices Dand Nwill be D=    n11 000 0n22 0 0 0 0 n33 0 000n44    , N =    n11 0 0 n14 0n22 0 0 0 0 n33 0 0n42 0n44    , where n42 and n14 are the new (optimal) entries in N, not appearing in D. Hence S4=span {M1,1, M2,2, M3,3, M4,4, M4,2, M1,4}. 7 Remark 3.2. Note that the preconditioning matrix Ngeneralizes the optimal diagonal approximate inverse D. Indeed, in the special case that the optimal entry for each column j= 1,2, . . . , n is the diagonal one, we have ij=jfor all j= 1,2, . . . , n ⇒N=D. Moreover, the preconditioner Nimproves, in general, the optimal diagonal preconditioner D. Indeed, since Dand Nare the solutions to problem (1.4) for the subspaces Dn=span {Mj,j}n j=1 and Sn=span {Mj,j}n j=1 ∪Mij,j |ij6=jn j=1, respectively, then we have Sn⊇ Dn⇒ kAN −IkF≤ AD −I F. The following is the main result of this paper. It provides us with explicit expressions for both matrix Nand the minimum Frobenius norm kAN −IkF. Theorem 3.1. Let A∈Rn×nbe nonsingular. Let Nbe the solution to problem (1.4) for the subspace Sn=span {Mj,j}n j=1 ∪Mij,j |ij6=jn j=1.(3.6) Then, for each j= 1,2, . . . , n, its corresponding index ijis defined by the condition ajij  Aeij 2 = max |aj1| kAe1k2 ,|aj2| kAe2k2 ,··· ,|ajn| kAenk2.(3.7) Moreover, N= n X j= 1 ij=j ajj kAejk2 2 Mj,j + n X j= 1 ij6=j ajj  Aeij 2 2−ajijAej, Aeij2 kAejk2 2 Aeij 2 2−Aej, Aeij2 2 Mj,j + n X j= 1 ij6=j ajijkAejk2 2−ajj Aej, Aeij2 kAejk2 2 Aeij 2 2−Aej, Aeij2 2 Mij,j (3.8) 8 and the corresponding minimum Frobenius norm is given by kAN −Ik2 F=n− n X j=1 a2 jj kAejk2 2 − n X j= 1 ij6=j ajijkAejk2 2−ajj Aej, Aeij22 kAejk2 2kAejk2 2 Aeij 2 2−Aej, Aeij2 2.(3.9) Proof. First, we determine the (optimal) positions {(ij, j)}n j=1 of the nonzero entries in the best approximate inverse of matrix Aamong all n×nmatrices that have exactly one nonzero element per column. Let j∈ {1,2, . . . , n}be arbitrary, but fixed. The optimal approximate inverse Ni,j, among all the n×nmatrices whose only nonzero term is placed at the ith row, jth column, can be obtained as the solution to problem (1.4) for the one-dimensional subspace S=span {Mi,j}. That is, using Eqs. (2.1) and (3.1), we obtain Ni,j =tr (AMi,j) kAMi,jk2 F Mi,j =aji kAeik2 2 Mi,j, for which, using Eqs. (2.2) and (3.1), we have kANi,j −Ik2 F=n−[tr (AMi,j)]2 kAMi,jk2 F =n−a2 ji kAeik2 2 . Consequently, the index i∈ {1,2, . . . , n}that minimizes kANi,j −Ik2 Ffor each fixed column j, is the one that maximizes the quotient a2 ji kAeik2 2 , that is, the index ij, defined by Eq. (3.7). Now, consider the set T={j∈ {1,2, . . . , n} | ij6=j}. There are two possible cases. Case 1. If ij=jfor all j= 1,2, . . . , n then T=∅. In this case, Sn=span {Mj,j}n j=1 ∪Mij,j |ij6=jn j=1=span {Mj,j}n j=1 =Dn. 9 Proof. Using the obvious fact that kA−Ik2 F=n−2tr (A)−kAk2 F and Eqs. (3.9) and (3.15), we obtain kA−Ik2 F−kAN −Ik2 F=kAk2 F−2tr (A) + n X j=1 a2 jj kAejk2 2 + n X j= 1 ij6=j ajijkAejk2 2−ajj Aej, Aeij22 kAejk2 2kAejk2 2 Aeij 2 2−Aej, Aeij2 2 ≥ n X j=1 kAejk2 2− n X j=1 2ajj + n X j=1 a2 jj kAejk2 2 + n X j= 1 ij6=j ajij  Aeij 2−|ajj| kAejk2!2 1 sin2θj = n X j=1 kAejk2 2−ajj2 kAejk2 2 + n X j= 1 ij6=j ajij  Aeij 2−|ajj| kAejk2!2 1 sin2θj . Finally, if kAejk2 26=ajj for at least one index j= 1,2, . . . , n then the left sum in Eq. (3.17) contains at least one positive summand, and thus we conclude that kA−IkF>kAN −IkF. Moreover, if matrix Nis not diagonal then ij6=jfor at least one column j∈ {1,2, . . . , n}(case 2 in the proof of Theorem 3.1). Hence, the right sum in Eq. (3.17) contains at least one summand, which is necessarily positive due to Eq. (3.7), and thus kA−IkF>kAN −IkF. Remark 3.6. Corollary 3.1 has established the comparison between the optimal diagonal preconditioner Dand the proposed preconditioner N, in the 16 following terms. First, note that from Eqs. (3.7) and (3.16), we conclude that Nimproves Din the sense of Definition 3.1, i.e., If N6=D⇒ ∃ j∈ {1,2, . . . , n}s.t. ij6=j⇒ AD −I F>kAN −IkF. Second, Eq. (3.16) provides us with the following analysis of the improvement achieved when using preconditioner Ninstead of preconditioner D. For the second item, we use the obvious fact that the function f(θ) = 1 sin2θis strictly decreasing in the interval 0,π 2, and strictly increasing in the interval π 2, π. (i) For each j= 1,2, . . . , n such that ij6=j, the more the maximum quotient |ajij| kAeijk2 (given by Eq. (3.7)) exceeds the quotient |ajj | kAejk2 , the larger the difference  AD −I 2 F−kAN −Ik2 Fwill be, and thus the more the preconditioner Nimproves the diagonal preconditioner D(in the sense of Definition 3.1). (ii) For each j= 1,2, . . . , n such that ij6=j, the closer the angle θjbetween the jth and the ijth columns of the coefficient matrix Ais either to 0 or to π(i.e., the larger the difference between θjand π 2is), the larger f(θj), and then the larger the difference  AD −I 2 F−kAN −Ik2 Fwill be, and thus the more the preconditioner Nimproves the diagonal preconditioner D(in the sense of Definition 3.1). Remark 3.7. Note that the right-most sum in Eq. (3.17) coincides with the sum in Eq. (3.16). Thus, the above two comments (i) and (ii) in Remark 3.6, concerning the comparison between  AD −I Fand kAN −IkF(analyzed in Corollary 3.1), remain true for the comparison between kA−IkFand kAN −IkF(analyzed in Corollary 3.2). Remark 3.8. Call N1=Nthe approximate inverse of matrix A, constructed in Theorem 3.1. According to the well-known multistep preconditioning strategy (see, e.g., [15, 26]), we can obtain a sequence N1,N1N2, N1N2N3, . . . of approximate inverses of Awhere, for every k≥2, matrix Nk is the best sparse approximate inverse of matrix AN1N2···Nk−1, among all matrices defined by the sparsity pattern (3.6). Note that since subspace Dnof all n×ndiagonal matrices is closed for the matrix product, then we have min M∈Dnk(AN1)M−IkF=k(AN1)N2−IkF=kA(N1N2)−IkF min M∈DnkA(N1M)−IkF= min M∈DnkAM −IkF=kAN1−IkF, 17 and then, due to the uniqueness of solution of problem (1.4), we conclude that N1N2=N1. This means that the multistep strategy does not make sense for the optimal diagonal preconditioner. However, this does not happen with the optimal preconditioners Nbelonging to the subspaces Sndefined by the prescribed sparsity patterns (3.6) and thus, in our case the multistep preconditioning strategy makes sense. In particular, Eq. (3.18) implies that each step kof our multistep preconditioning strategy strictly reduces the Frobenius norm, whenever the preconditioner Nk−1, obtained in the previous step, is not diagonal. 4. Numerical experiments We present some numerical experiments to illustrate the behavior of the proposed preconditioner N. We compare the preconditioned linear system using our preconditioner Nwith both the unpreconditioned linear system and the preconditioned system using the optimal diagonal preconditioner D. At the end of this section, the preconditioner Nis also compared with the approximate inverse preconditioner AINV. We have studied a number of linear systems Ax =b, where the test coefficient matrices are taken from the University of Florida Sparse Matrix Collection [28]. We carried out our numerical problems with the Krylov solvers GMRES [29] and BiCGStab [30]. Both solvers led to similar results for most test matrices, but a small advantage was observed when using the latter for solving the systems preconditioned with matrix N. For this reason, we only present here the results obtained with the (right-preconditioned) BiCGStab. In any case, our purpose in this paper is to analyze the effectiveness of the proposed preconditioner (especially in comparison with the optimal diagonal approximate inverse), rather than to compare different Krylov subspace methods. The initial guess was always x0= 0, and the right-hand side vector was b= [1,...,1]T. The stopping criterion was either kb−Axkk2 kbk2 <10−8, or when this condition about the relative residual was not satisfied, within 2niterations (nbeing the order of the coefficient matrix A). We run all numerical experiments in double precision arithmetic, on Intel(R)Xeon(R) E5620 with 2.40 GHz clock frequency and 24GB of main memory using GNU Octave 3.2.4. 18 In Table 1, nand nnz (A) stand for the order and the number of nonzero entries of matrix A, respectively. This table also provides the Frobenius norms of matrices A−I,AD −Iand AN −I. In Table 2, nnz (N) denotes the number of nonzero entries of our preconditioner N, which is compared, in its second column, with the number nof entries of the optimal diagonal preconditioner D(obviously, n≤nnz (N)≤ 2n). In the third and fourth columns, D-time and N-time denote the CPU time (in seconds) for constructing the preconditioners Dand N, respectively. In the three right-most columns, Unprec-iter, D-iter, and N-iter stand for the number of iterations of the BiCGStab method for the unpreconditioned system, and the preconditioned systems with Dand N, respectively. When convergence is not attained, within the maximum number 2nof allowed iterations, we indicate it by writing “†”, in any of the corresponding columns Unprec-iter, D-iter and N-iter. Numerical tests reported confirm the theoretical results and illustrate the effectiveness of the proposed preconditioner in comparison with the optimal diagonal one. Test problems have been grouped together into five classes, according to the behavior of the preconditioner Nin comparison with D. Looking at the two right-most columns (D-iter and N-iter) in Table 2, one can easily identify each of these five groups of test matrices. For each of these classes, problems are arranged in increasing order of the size nof the test coefficient matrices. The first five test problems correspond to matrices for which the optimal diagonal preconditioner and our preconditioner coincide, i.e., N=D. Of course, this is in accordance with the theory, since Nhas been defined as a generalization of D, and they do coincide when ij=jfor all j= 1,2, . . . , n. Obviously, in such cases, Table 2 shows that nnz (N) = n(each column of Nconsists only of its diagonal entry), and the number of iterations needed for convergence coincide for both preconditioners. The rest of test matrices corresponds to the case N6=D, so that the number nnz (N) of nonzero entries of Nwill be greater than n. When nnz (N) = 2n, this means that ij6=jfor all j= 1,2, . . . , n, and each column of Nconsists of two nonzero entries. 19 Table 1 The test matrices and the Frobenius norms of A−I,AD −Iand AN −I. Matrix n nnz (A)kA−IkF AD −I FkAN −IkF orsirr 2 886 5970 1.55 ×10618 18 sherman1 1000 3750 1.00 ×10315.6 15.6 sherman4 1104 3786 1.21 ×10310.4 10.4 bcsstk09 1083 18437 8.57 ×10821.2 21.2 sherman3 5005 20033 1.36 ×10727.2 27.2 hor 131 434 4182 4.34 ×10214.7 14.7 rdb450l 450 2580 5.01 ×10216.2 13.7 pores 3 532 3474 6.63 ×10515.1 14.9 steam2 600 5660 5.27 ×1010 16.8 12.3 young3c 841 3988 6.40 ×10314.8 14.7 bcsstk10 1086 22070 2.97 ×10822.6 22.6 olm500 500 1996 2.24 ×10518.3 15.6 olm1000 1000 3996 1.26 ×10625.8 22.1 tols1090 1090 3546 1.23 ×10716.8 14.6 fpga trans 01 1220 7382 1.22 ×10324.6 22.4 adder dcop 14 1813 11246 1.81 ×10324.3 23.9 adder dcop 15 1813 11246 1.81 ×10324.4 24 adder dcop 16 1813 11246 1.81 ×10324.1 23.7 adder dcop 17 1813 11246 1.81 ×10324 23.6 adder dcop 20 1813 11246 1.81 ×10324 23.5 adder trans 02 1814 14579 1.81 ×10314.6 14.1 tols2000 2000 5184 5.40 ×10721.6 18.7 psmigr 1 3140 543160 3.54 ×10615.4 15.4 tols4000 4000 8784 2.98 ×10829.5 25.5 meg4 5860 25258 9.87 ×10517.4 15.8 steam1 240 2248 8.41 ×1078.95 8.92 sherman5 3312 20793 1.44 ×10432.4 32.3 sherman2 1080 23094 7.00 ×10930.7 28.9 adder dcop 10 1813 11232 1.81 ×10325.2 24.7 20 Table 2 Convergence results for Dand N. Matrix nnz (N)/n D-time N-time Unprec-iter D-iter N-iter orsirr 2 886/886 0.072 0.360 1232 448 448 sherman1 1000/1000 0.084 0.368 407 192 192 sherman4 1104/1104 0.096 0.408 92 64 64 bcsstk09 1083/1083 0.124 0.832 188 154 154 sherman3 5005/5005 1.56 6 †453 453 hor 131 435/434 0.028 0.140 †369 255 rdb450l 900/450 0.024 0.380 †254 51 pores 3 560/532 0.032 0.188 †717 637 steam2 750/600 0.044 0.324 448 10 7 young3c 849/841 0.072 0.304 1080 974 804 bcsstk10 1088/1086 0.136 0.936 †736 657 olm500 1000/500 0.032 0.428 † † 307 olm1000 2000/1000 0.084 0.996 † † 657 tols1090 1568/1090 0.088 0.688 † † 1560 fpga trans 01 1352/1220 0.128 0.684 † † 205 adder dcop 14 1834/1813 0.248 1.224 † † 437 adder dcop 15 1836/1813 0.248 1.220 † † 310 adder dcop 16 1837/1813 0.248 1.224 † † 418 adder dcop 17 1843/1813 0.248 1.228 † † 329 adder dcop 20 1857/1813 0.248 1.248 † † 518 adder trans 02 1832/1814 0.260 1.384 229 †9 tols2000 2842/2000 0.256 1.560 † † 1534 psmigr 1 3143/3140 4.348 46.951 1812 †58 tols4000 5642/4000 0.896 4.492 † † 1328 meg4 5968/5860 2.104 8.584 † † 5 steam1 320/240 0.012 0.108 †236 388 sherman5 3327/3312 0.756 3.53 2251 849 892 sherman2 1628/1080 0.144 1.38 † † † adder dcop 10 1872/1813 0.248 1.256 † † † In some of these cases (e.g., for the test matrices displayed in rows 6-11), the number of iterations needed to reach convergence is reduced when using our preconditioner Ninstead of D. Moreover, in many cases, like for instance for the test matrices displayed in rows 12-25 in Table 2, convergence is not attained with the diagonal pre21 conditioner D, but however it is reached when we use our preconditioner N. We think that this is the main advantage of our preconditioner pointed out by the numerical tests. In particular, we highlight the case of the test matrix psmigr 1. In this case, although the preconditioner Nhas only 3 nonzero entries more than the diagonal preconditioner D(nnz (N)/n = 3143/3140), the system transits from non-convergence to convergence if we use N(with a small number of iterations) instead of Das preconditioning matrix. Moreover, for this test problem, both preconditioners only differ in the numerical values of 6 entries: the nnz (N)−n= 3 diagonal entries (j, j) for which ij6=j, and the corresponding 3 nondiagonal entries (ij, j) that do not appear in the diagonal matrix D. The numerical values of the remaining n−3 = 3137 nonzero diagonal entries of Dand Ncoincide (see Eqs. (3.4) and (3.8)). So, Dand Nare very similar but, as mentioned, we obtain convergence with the latter, but not with the former. This happens not only for the test problem psmigr 1, but also for other test matrices reported. Next, we present two problems, namely the third and fourth ones from the bottom, for which the number of iterations required for reaching convergence is greater when we use the preconditioner Ninstead of D. Finally, for the last two test matrices reported, convergence is reached neither with Nnor with D. This can be explained by the small number of nonzero entries in both preconditioners, which makes them inefficient for some very ill-conditioned systems. For all test problems reported, the last three columns in Table 1 confirm our theoretical results, in the sense that kA−IkF≥ AD −I F≥ kAN −IkF. On one hand, the first inequality is due to the fact that the optimal diagonal preconditioner D(by definition) minimizes kAM −IkFover the subspace of all n×ndiagonal matrices (and the identity is a diagonal matrix). On the other hand, the second inequality is an obvious consequence of the set inclusion Sn⊇ Dn(as commented in Remark 3.2). The large difference (observed for all test matrices reported) between kA−IkFand any of the values  AD −I Fand kAN −IkFis mainly due to the following fact. While, obviously, kA−IkFcan be arbitrarily large, both norms  AD −I Fand kAN −IkFare never greater than √n. The reason is that, from Eq. (2.2) we immediately derive that the optimal approximate 22 inverse M(in the Frobenius sense) over any matrix subspace S ⊂ Rn×n (and, in particular, M=Dand M=N) always satisfies kAM −IkF≤√n. Also, we can observe that the test matrices for which the difference between  AD −I Fand kAN −IkFis larger usually correspond to those cases for which the ratio between the numbers of nonzero entries in matrices Nand D (denoted by nnz (N)/n in Table 2) is larger. In addition to this ratio, other parameters that also determine the difference between (the squares of) the Frobenius norms of AD −Iand AN −Ihave been analyzed in Remark 3.6. Summing up, in almost all the cases, the iterations required by the BiCGStab method when using preconditioner Nare fewer than those needed by this solver when using the diagonal preconditioner D. The differences between the CPU times for constructing the preconditioners Dand N, as well as the differences between the execution times of the BiCGStab method for both preconditioners, are not significant, and in most cases the total CPU times for solving the system with Dand Nare of the same order of magnitude. In any case, the small increment in computational cost when using the preconditioner Ninstead of D, is compensated by the fact that, in many cases, convergence is reached with the preconditioner Nbut not with the optimal diagonal preconditioner D. For both preconditioners Dand N, their respective sparsity patterns consist of small numbers of nonzero entries (nentries for D, and a number of entries between nand 2nfor N). On one hand, this implies low computational costs for constructing them. On the other hand, for both of them, the number of required iterations is large in comparison with other more dense and expensive optimal approximate inverse preconditioners based on the same idea (Frobenius norm minimization). In any case, as the numerical experiments have shown, the proposed preconditioner improves the classical diagonal one, with a small increment in the computation cost required for constructing the former instead of the latter. To finish this section, we compare the proposed preconditioner Nwith a more expensive approximate inverse preconditioner, namely the well-known AINV preconditioner [6]. We have implemented the AINV preconditioner, with a drop tolerance Tol = 0.25, for the same set of test matrices used for comparing our preconditioner Nwith the optimal diagonal one D. For most of these test problems, AINV was found to be more efficient (in terms of the overall solution time and using the Krylov solver BiCGStab) than N. However, for some test matrices, namely sherman4,steam2,adder trans 02, 23 psmigr 1 and meg4, the proposed preconditioner was found to be more efficient than AINV. In conclusion, our proposed preconditioner was more effective than the optimal diagonal Dfor most of the numerical problems considered in this paper, and it was more effective than the AINV preconditioner in a few cases. 5. Summary and conclusions In this paper, a new approximate inverse preconditioner Nfor large sparse linear systems has been constructed and theoretically analyzed. Nhas been defined as the optimal preconditioner (in the Frobenius sense) among all the n×nmatrices whose only nonzero entries, for each column j= 1,2, . . . , n, are the diagonal one (j, j) and, in addition, the optimal entry (ij, j) in column j, whenever it does not coincide with the diagonal one. In this way, our preconditioning matrix Ngeneralizes the optimal diagonal preconditioner D. Explicit expressions for both matrix Nand the minimum Frobenius norm kAN −IkFhave been presented. We have proved that, whenever N6=D, the preconditioner N(which has at least nand at most 2nnonzero entries) improves D, in the sense of the Frobenius norm. We have also analyzed the difference between the Frobenius norms of A−Iand AN −I. Numerical experiments have confirmed the theoretical results, presenting a number of test matrices for which the proposed preconditioner improves the convergence of the optimal diagonal one, when they do not coincide. In particular, Table 1 shows that the Frobenius norm kAN −IkFis always smaller (as we have theoretically shown) and, in fact, much smaller in most cases than the Frobenius norm kA−IkF. Table 2 illustrates the reduction, in most cases, of the number of iterations when we use the preconditioning matrix Ninstead of D. The main advantage of our preconditioner pointed out by the numerical tests is the following. For many test matrices, the small additional CPU time required for constructing Ninstead of Dis compensated by the fact that system transits from non-convergence to convergence when using N instead of Das preconditioning matrix. Moreover, this also occurs for some test problems for which the number of nonzero entries of the preconditioner Nexceeds only by a very small quantity the number of nonzero entries of the diagonal preconditioner D, and the numerical values of most of the entries placed at the same diagonal positions coincide for both preconditioners (D and Nare very similar). 24 For future researches, it would be interesting to analyze in more detail the practical value of the preconditioner proposed in this paper. For this purpose, our preconditioner Ncould be compared, using new test matrices, with the AINV approximate inverse preconditioner and with other preconditioners not considered in this paper. In addition, it is worth trying to study some common characteristics/features of those test matrices for which the generalized preconditioner Nhas a better behavior for convergence purposes. Finally, regarding an additional line for future work, our method can be improved by considering the optimal preconditioner N(in the sense of the Frobenius norm) among all the n×nmatrices having exactly a (small) fixed number mj≥2 of nonzero entries for each column j= 1,2, . . . , n. The determination of the sparsity pattern of such preconditioner Nis not a simple problem because of the following fact. Assume that the jth column Nejof the preconditioner Nconsists of only one nonzero entry, i.e., mj= 1. Then, as shown in Section 3, the optimal position i1 j, jin the jth column Nejof Nfor minimizing kANej−ejk2can be easily determined simply by using Eq. (3.7). Similarly, we can easily determine the second, third,...,mjth best positions i2 j, j,i3 j, j,...,imj j, jin the jth column of Nfor minimizing kANej−ejk2. Unfortunately, the set of these best mjpositions (obtained separately, i.e., when Nejconsists of only one nonzero entry) for approximating the jth column ejof I, does not necessarily coincide with the optimal set of cardinality mjfor approximating ej(when the sparsity pattern of N is defined by the condition that Nejconsists of mj≥2 nonzero entries). Consequently, the determination of this optimal mj-set (in order to find the optimal sparsity pattern for an optimal preconditioner defined by Eq. (1.4)) is a very difficult problem (when mj≥2). Alternatively, one can consider the possibility of using an algorithm based on the LU factorization with partial pivoting to determine the optimal sparsity pattern for each column. Acknowledgments The authors thank the anonymous referees for their detailed revisions and valuable comments and suggestions, which have substantially improved the earlier version of this paper. This work was partially supported by the “Ministerio de Econom´ıa y Competitividad” of the Spanish Government, and FEDER, through Grant contract: CGL2011-29396-C03-01. 25