- Email: [email protected]
On the ijk forms of the truncated LU decomposition
Journal of Computational and Applied Mathematics 233 (2009) 582–584
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Letter to the editor
On the ijk forms of the truncated LU decomposition A. Mohsen ∗ Engineering Mathematics & Physics Department, Engineering Faculty, Cairo University, Giza 12211, Egypt
article
info
Article history: Received 27 November 2008 Received in revised form 30 April 2009 Keywords: Matrix LU-decomposition The ijk forms Fredholm integral equations of the first kind Truncated LU-decomposition Vector and parallel processing
abstract Matrix LU decomposition has six ijk forms. Different forms have different computational complexities and storage requirements, particularly on vector and parallel computers. Other factors governing the choice of a particular form are considered. For treating Fredholm integral equations of the first kind, the truncated LU decomposition of the resulting system matrix is recommended. Required modifications to selected known ijk forms are presented. © 2009 Elsevier B.V. All rights reserved.
1. Introduction The solution of linear systems of equations is a basic step in most calculations in scientific computing. In general, matrix factorization is almost always the first step of much scientific computing and usually the one which places the heaviest demand in terms of computing resources. These include linear system solution, eigenvalue and least squares approximations and rank revealing problems. In solving linear systems, the sequential algorithm that performs Gauss elimination has a main loop with three nesting levels (indices i, j and k) that can be arranged according to six different organizations (called ijk forms) [1]. These variants correspond to different permutations of sequences of rows and columns computations of L and U. If all submatrices of the matrix A(n × n) are nonsingular, the LU factorization exists and is unique [2, Th.3.2.1]. While two implementations may be equivalent in terms of the number of operations, there may be substantial differences due to architectural considerations. A detailed exposition of the ijk forms of LU as well as Cholesky factorizations was considered by Ortega [3,4]. Several aspects of these forms and their properties on vector computers [3] as well as on parallel computers [4] were studied. Parallelization of the ijk forms has also been considered in detail in [5]. The efficiency of a particular ijk form depends on the sparsity of the system matrix. For the efficient iterative solution of sparse linear systems the choice of a proper preconditioner is crucial [6]. Incomplete LU (ILU) and sparse approximate inverse (SAI) represent two main preconditioners. Both set a threshold where one discards elements below this threshold to compute a sparse preconditioner. Different ijk forms will yield different preconditioners with different computational complexities and accuracies. This aspect should be carefully investigated. For example, the modified ILU (MILU) due to Gustafson [7] requires the row-sum in order to stabilize the ILU. In this case the forms containing row(U ) are more appropriate than the remaining forms.
∗
Tel.: +20 235678051. E-mail address: [email protected].
0377-0427/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2009.07.050
A. Mohsen / Journal of Computational and Applied Mathematics 233 (2009) 582–584
583
2. Forms suitable for truncated LU decomposition of a square matrix We consider the Fredholm integral equation of the first kind: f ( x) =
Z
K (x, t )u(t )dt .
This equation appears in many engineering applications including inverse problems and image restoration. It is widely used to estimate parameters and map geometries in areas as diverse as geophysical prospecting, medical imaging and nondestructive evaluation [8]. Replacing the integral equation by a quadrature often results in an ill-conditioned set of linear equations: Au = b. In such cases a reasonable result may be obtained when the approximating algebraic system is small. The solution oscillates violently when the size of the system is increased in an attempt to reduce the truncation error. A review of early methods used to treat the problem may be found in [9]. Of these methods, the truncated LU decomposition [10] is much simpler than the eigenvalue and svd methods. It needs fewer arithmetic operations and less computational time. Thus, we decompose A as A = L U. When L is unit lower triangular, any ill-conditioning in the matrix appears in d = diagonal of U [11]. With proper permutation, d has a decreasing value. Thus it plays a role similar to that of the eigenvalue of A. Let eps be the threshold for the truncation; we present appropriate ijk forms for finding the truncated LU decomposition when A is square. In this case the ijk forms containing row(U ) are more appropriate than the remaining forms. We first consider the ikj form which computes a row of L followed by a row of U. We use the MATLAB notation. Row(L) followed by row(U ):
for i =1: n for j =1: i-1,s = 0; for k =1: j-1,s = s + L(i,k)*U(k,j); end L(i,j) =(A(i,j)-s ) / U(j,j); end s=0 ; for k=1: i-1,s = s + L(i,k)*U(i,k); end U(i,i) = A(i,i) - s; if (abs (U(i,i)) < eps) , TRUNC=i , break ,end for j = i+1: n, s = 0; for k = 1: i-1, s = s +L(i,k)*U(k,j);end U(i,j) = A(i,j) - s; end L( i,i) = 1; end Similarly, we may compute the truncated LU decomposition using row(U ) followed by row(L). For the more popular kij form corresponding to row(U ) followed by column(L) we have: Row(U ) followed by column(L): for k=1:n-1 s= 0 ; for l =1:k-1,s = s + L(k,l)*U(l,k);end U(k,k) = A(k,k) - s ; if(abs(U(k,k))< ep),TRUNC=k , break , end for j = k+1:n , s=0 ; for l = 1:k-1, s = s + L(k,l)*U(l,j);end U(k,j)=A(k,j)-s; end L(k,k) = 1; for i = k+1: n,s = 0; for l = 1:k-1, s = s +L(i,l)*U(l,k);end L(i,k) = (A(i,k)-s) / U(k,k); end end Similarly, we may compute the truncated LU decomposition using column(L) followed by row(U ). 3. Conclusion For treating Fredholm integral equations of the first kind, the truncated LU decomposition of the resulting system matrix is recommended. Required modifications to selected ijk forms are presented. The four forms based on row(U ) computations are appropriate since the truncation is based on testing whether v = abs(U (i, i)) < eps. Thus, we calculate this quantity before we proceed further to compute U (i, j), j = i + 1 : n. Different forms have different computational complexities and storage requirements. Factors governing the choice of a particular form depend on the type of the system matrix, the computational environment and the application. For example, kij, r (row(U ) followed by column(L) with row interleaved storage) is recommended for parallel computation [3,4].
584
A. Mohsen / Journal of Computational and Applied Mathematics 233 (2009) 582–584
References [1] [2] [3] [4] [5] [6] [7] [8]
J.J. Dongarra, et al., Implementing linear algebra algorithms for dense matrices on a vector pipeline machine, SIAM Rev. 26 (1984) 91–112. G.H. Golub, C.F. Van Loan, Matrix Computations, The John Hopkins University Press, Baltimore, 1996 (Ch. 3). J.M. Ortega, The ijk forms of factorization methods: I-Vector computers, Parallel Comput. 7 (1988) 135–147. J. Ortega, C. Romine, The ijk forms of factorization methods: II Parallel computers, Parallel Comput. 7 (1988) 149–162. T.L. Freman, C. Phillips, Parallel Numerical Algorithms, Prentice Hall, NY, 1992. Y. Saad, H.A. van der Vorst, Iterative solution of linear systems in the 20th century, J. Comput. Appl. Math. 123 (2000) 1–33. I. Gustafson, A class of first order factorization methods, BIT Numer. Math. 18 (1978) 142–156. D.G. Dudley, T.M. Habashy, E. Wolf, Linear inverse problems in wave motion: Nonsymmetric first kind integral equations, IEEE Trans. Antennas and Propagation 48 (2000) 1607–1617. [9] A.V. Bitsadze, Integral Equations of First Kind, World Sci. Pub. Co., 1995. [10] N.N. Abdelmalek, T. Kasvand, Image restoration by Gauss LU decomposition, Appl. Opt. 18 (1979) 1684–1686. [11] G. Peters, J.H. Wilkinson, The least squares problem and pseudo-inverses, The Computer J. 13 (1970) 309–316.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Letter to the editor
On the ijk forms of the truncated LU decomposition A. Mohsen ∗ Engineering Mathematics & Physics Department, Engineering Faculty, Cairo University, Giza 12211, Egypt
article
info
Article history: Received 27 November 2008 Received in revised form 30 April 2009 Keywords: Matrix LU-decomposition The ijk forms Fredholm integral equations of the first kind Truncated LU-decomposition Vector and parallel processing
abstract Matrix LU decomposition has six ijk forms. Different forms have different computational complexities and storage requirements, particularly on vector and parallel computers. Other factors governing the choice of a particular form are considered. For treating Fredholm integral equations of the first kind, the truncated LU decomposition of the resulting system matrix is recommended. Required modifications to selected known ijk forms are presented. © 2009 Elsevier B.V. All rights reserved.
1. Introduction The solution of linear systems of equations is a basic step in most calculations in scientific computing. In general, matrix factorization is almost always the first step of much scientific computing and usually the one which places the heaviest demand in terms of computing resources. These include linear system solution, eigenvalue and least squares approximations and rank revealing problems. In solving linear systems, the sequential algorithm that performs Gauss elimination has a main loop with three nesting levels (indices i, j and k) that can be arranged according to six different organizations (called ijk forms) [1]. These variants correspond to different permutations of sequences of rows and columns computations of L and U. If all submatrices of the matrix A(n × n) are nonsingular, the LU factorization exists and is unique [2, Th.3.2.1]. While two implementations may be equivalent in terms of the number of operations, there may be substantial differences due to architectural considerations. A detailed exposition of the ijk forms of LU as well as Cholesky factorizations was considered by Ortega [3,4]. Several aspects of these forms and their properties on vector computers [3] as well as on parallel computers [4] were studied. Parallelization of the ijk forms has also been considered in detail in [5]. The efficiency of a particular ijk form depends on the sparsity of the system matrix. For the efficient iterative solution of sparse linear systems the choice of a proper preconditioner is crucial [6]. Incomplete LU (ILU) and sparse approximate inverse (SAI) represent two main preconditioners. Both set a threshold where one discards elements below this threshold to compute a sparse preconditioner. Different ijk forms will yield different preconditioners with different computational complexities and accuracies. This aspect should be carefully investigated. For example, the modified ILU (MILU) due to Gustafson [7] requires the row-sum in order to stabilize the ILU. In this case the forms containing row(U ) are more appropriate than the remaining forms.
∗
Tel.: +20 235678051. E-mail address: [email protected].
0377-0427/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2009.07.050
A. Mohsen / Journal of Computational and Applied Mathematics 233 (2009) 582–584
583
2. Forms suitable for truncated LU decomposition of a square matrix We consider the Fredholm integral equation of the first kind: f ( x) =
Z
K (x, t )u(t )dt .
This equation appears in many engineering applications including inverse problems and image restoration. It is widely used to estimate parameters and map geometries in areas as diverse as geophysical prospecting, medical imaging and nondestructive evaluation [8]. Replacing the integral equation by a quadrature often results in an ill-conditioned set of linear equations: Au = b. In such cases a reasonable result may be obtained when the approximating algebraic system is small. The solution oscillates violently when the size of the system is increased in an attempt to reduce the truncation error. A review of early methods used to treat the problem may be found in [9]. Of these methods, the truncated LU decomposition [10] is much simpler than the eigenvalue and svd methods. It needs fewer arithmetic operations and less computational time. Thus, we decompose A as A = L U. When L is unit lower triangular, any ill-conditioning in the matrix appears in d = diagonal of U [11]. With proper permutation, d has a decreasing value. Thus it plays a role similar to that of the eigenvalue of A. Let eps be the threshold for the truncation; we present appropriate ijk forms for finding the truncated LU decomposition when A is square. In this case the ijk forms containing row(U ) are more appropriate than the remaining forms. We first consider the ikj form which computes a row of L followed by a row of U. We use the MATLAB notation. Row(L) followed by row(U ):
for i =1: n for j =1: i-1,s = 0; for k =1: j-1,s = s + L(i,k)*U(k,j); end L(i,j) =(A(i,j)-s ) / U(j,j); end s=0 ; for k=1: i-1,s = s + L(i,k)*U(i,k); end U(i,i) = A(i,i) - s; if (abs (U(i,i)) < eps) , TRUNC=i , break ,end for j = i+1: n, s = 0; for k = 1: i-1, s = s +L(i,k)*U(k,j);end U(i,j) = A(i,j) - s; end L( i,i) = 1; end Similarly, we may compute the truncated LU decomposition using row(U ) followed by row(L). For the more popular kij form corresponding to row(U ) followed by column(L) we have: Row(U ) followed by column(L): for k=1:n-1 s= 0 ; for l =1:k-1,s = s + L(k,l)*U(l,k);end U(k,k) = A(k,k) - s ; if(abs(U(k,k))< ep),TRUNC=k , break , end for j = k+1:n , s=0 ; for l = 1:k-1, s = s + L(k,l)*U(l,j);end U(k,j)=A(k,j)-s; end L(k,k) = 1; for i = k+1: n,s = 0; for l = 1:k-1, s = s +L(i,l)*U(l,k);end L(i,k) = (A(i,k)-s) / U(k,k); end end Similarly, we may compute the truncated LU decomposition using column(L) followed by row(U ). 3. Conclusion For treating Fredholm integral equations of the first kind, the truncated LU decomposition of the resulting system matrix is recommended. Required modifications to selected ijk forms are presented. The four forms based on row(U ) computations are appropriate since the truncation is based on testing whether v = abs(U (i, i)) < eps. Thus, we calculate this quantity before we proceed further to compute U (i, j), j = i + 1 : n. Different forms have different computational complexities and storage requirements. Factors governing the choice of a particular form depend on the type of the system matrix, the computational environment and the application. For example, kij, r (row(U ) followed by column(L) with row interleaved storage) is recommended for parallel computation [3,4].
584
A. Mohsen / Journal of Computational and Applied Mathematics 233 (2009) 582–584
References [1] [2] [3] [4] [5] [6] [7] [8]
J.J. Dongarra, et al., Implementing linear algebra algorithms for dense matrices on a vector pipeline machine, SIAM Rev. 26 (1984) 91–112. G.H. Golub, C.F. Van Loan, Matrix Computations, The John Hopkins University Press, Baltimore, 1996 (Ch. 3). J.M. Ortega, The ijk forms of factorization methods: I-Vector computers, Parallel Comput. 7 (1988) 135–147. J. Ortega, C. Romine, The ijk forms of factorization methods: II Parallel computers, Parallel Comput. 7 (1988) 149–162. T.L. Freman, C. Phillips, Parallel Numerical Algorithms, Prentice Hall, NY, 1992. Y. Saad, H.A. van der Vorst, Iterative solution of linear systems in the 20th century, J. Comput. Appl. Math. 123 (2000) 1–33. I. Gustafson, A class of first order factorization methods, BIT Numer. Math. 18 (1978) 142–156. D.G. Dudley, T.M. Habashy, E. Wolf, Linear inverse problems in wave motion: Nonsymmetric first kind integral equations, IEEE Trans. Antennas and Propagation 48 (2000) 1607–1617. [9] A.V. Bitsadze, Integral Equations of First Kind, World Sci. Pub. Co., 1995. [10] N.N. Abdelmalek, T. Kasvand, Image restoration by Gauss LU decomposition, Appl. Opt. 18 (1979) 1684–1686. [11] G. Peters, J.H. Wilkinson, The least squares problem and pseudo-inverses, The Computer J. 13 (1970) 309–316.