# LU decomposition

(Difference between revisions)
 Revision as of 20:42, 15 December 2005 (view source)Tsaad (Talk | contribs) (towards a uniform notation for linear systems : A*Phi = B)← Older edit Latest revision as of 09:22, 21 November 2011 (view source)Yonah (Talk | contribs) (small typo fixed: first we solve intermediary vectory y and then x based on this y) (10 intermediate revisions not shown) Line 1: Line 1: - == LU Solvers == + == Description == - For the system of equations $A\cdot\phi=B$.
+ Consider the system of equations $A\vec{x}=\vec{b}$, where A[/itex] is an $n\times n$ nonsingular matrix.  $A$ may be decomposed into an lower triangular part $L$ and an upper triangular part $U$ that will lead us to a direct procedure for the solution of the original system.  This decomposition procedure is quite useful when more than one right-hand side (more than one $\vec{b}$) is to be used. - The solvers based on factorization are widely popular. The factorization of a non singular matrix '''A''' in two matrices '''L''' and '''U''', called lower and upper matrices, leads to a direct procedure for the evaluation of the inverse of the matrix. Thus making it possible to calculate the solution vector with given source vector. + - It now remains a two step process: a forward substitution followed by a backward substitution. + The algorithm is relatively straightforward - first, we determine the upper and lower triangular parts: + + :$A = LU.$ + + Then, + + :$A \vec{x} = (LU) \vec{x} = L(U\vec{x})=L\vec{y},$ + + where $y=Ux$.  Once we solve the system + + :$L\vec{y} = \vec{b},$ + + we will be able to find the solution to the original system by solving + + :$U\vec{x} = \vec{y}$ + + The first solution is a foward substitution, while the second solution is a backward substitution.  Both can be done efficiently once the factorization is available.  The forward substitution may be expressed as + + :$+ y_i = {1 \over {l_{ii} }}\left( {b_i - \sum\limits_{j = 1}^{i-1} {l_{ij} y_j } } \right), + i = 1,2,\ldots,n$ + + and the backward substitution process may be expressed as + + :$+ x_i = {1 \over {u_{ii} }}\left( {y_i - \sum\limits_{j = i + 1}^n {u_{ij} x_j } } \right), + i = n,n-1,\ldots,1.$ + + + LU decomposition essentially stores the operations of Gaussian elimination in "higher-level" form (see [[#References|Golub and Van Loan]]), so repeated solutions using the same left-hand side are computed without repetition of operations that are independent of the right-hand side. == Algorithm == == Algorithm == + ''Add factorization here''
+ Forward substitution
+ :    for k:=1 step until n do
+ ::    for i:=1 step until k-1
+ :::    $b_k=b_k-l_{ki}b_{i}$
+ ::    end loop (i)
+ ::    $b_{k}=b_{k}/l_{kk}$
+ :      end loop (k) + Backward substitution + :    for k:=n stepdown until 1 do
+ ::    for i:=k+1 step until n
+ :::    $b_k=b_k-u_{ki}b_{i}$
+ ::    end loop (i)
+ ::    $x_{k}=b_{k}/u_{kk}$
+ :      end loop (k) + + == Important Considerations == + As with [[Gaussian elimination]],  LU decomposition is probably best used for relatively small, relatively non-sparse systems of equations (with small and non-sparse open to some interpretation).  For larger and/or sparse problems, it would probably be best to either use an iterative method or use a direct solver package (e.g. [http://www.cse.psu.edu/~raghavan/software.html DSCPACK]) as opposed to writing one of your own. - :  Calculate '''LU''' factors of '''A''' + If one has a single lefthand-side matrix and many right-hand side vectors, then LU decomposition would be a good solution procedure to consider. At the very least, it should be faster than solving each system separately with Gaussian elimination. - ::$A\cdot \phi = (LU)\cdot \phi = L(U\cdot\phi)=LY=B$ + - : Solve $LY=B$ by forward substitution + - :: $Y = L^{-1} B$ + - :  Solve $U\phi =Y$ by backward substitution + - :: $\phi = U^{-1}Y$ + - {{stub}} + Here a Fortran code fragment for LU decomposition written by  D. Partenov - ---- + + do j = 1, n + do i = 1, j + if (i == 1) then + a(i,j) = a(i, j) + else + suma = 0.0 + do k = 1, i-1 + suma = suma + a(i,k)*a(k,j) + end do + a(i,j) = a(i,j) - suma + end if + end do + if ( j < n) then + do s = 1, n-j + i = j + s + if (j == 1) then + a(i,j) = a(i,j)/a(j,j) + else + suma = 0.0 + do k = 1, j-1 + suma = suma + a(i,k)*a(k,j) + end do + a(i,j) = (a(i,j) - suma)/a(j,j) + end if + end do + end if + end do + y(1) = b(1) + do i = 2, n + suma = 0.0 + do j = 1, i-1 + suma = suma + a(i,j)*y(j) + end do + y(i) = b(i) - suma + end do + i = n + j = n + x(i) = y(i)/a(i,j) + do s = 1, n -1 + i = n - s + suma = 0.0 + do j = i+1, n + suma = suma + a(i,j)*x(j) + end do + x(i) = (y(i) - suma)/a(i,i) + end do - ---- + == References == - Return to [[Numerical methods | Numerical Methods]] + *{{reference-book|author=Golub and Van Loan|year=1996|title=Matrix Computations|rest= 3rd edition, The Johns Hopkins University Press, Baltimore}} + *[http://en.wikipedia.org/wiki/LU_decomposition Wikipedia article ''LU decomposition]

## Description

Consider the system of equations $A\vec{x}=\vec{b}$, where $A$ is an $n\times n$ nonsingular matrix. $A$ may be decomposed into an lower triangular part $L$ and an upper triangular part $U$ that will lead us to a direct procedure for the solution of the original system. This decomposition procedure is quite useful when more than one right-hand side (more than one $\vec{b}$) is to be used.

The algorithm is relatively straightforward - first, we determine the upper and lower triangular parts:

$A = LU.$

Then,

$A \vec{x} = (LU) \vec{x} = L(U\vec{x})=L\vec{y},$

where $y=Ux$. Once we solve the system

$L\vec{y} = \vec{b},$

we will be able to find the solution to the original system by solving

$U\vec{x} = \vec{y}$

The first solution is a foward substitution, while the second solution is a backward substitution. Both can be done efficiently once the factorization is available. The forward substitution may be expressed as

$y_i = {1 \over {l_{ii} }}\left( {b_i - \sum\limits_{j = 1}^{i-1} {l_{ij} y_j } } \right), i = 1,2,\ldots,n$

and the backward substitution process may be expressed as

$x_i = {1 \over {u_{ii} }}\left( {y_i - \sum\limits_{j = i + 1}^n {u_{ij} x_j } } \right), i = n,n-1,\ldots,1.$

LU decomposition essentially stores the operations of Gaussian elimination in "higher-level" form (see Golub and Van Loan), so repeated solutions using the same left-hand side are computed without repetition of operations that are independent of the right-hand side.

## Algorithm

Forward substitution

for k:=1 step until n do
for i:=1 step until k-1
$b_k=b_k-l_{ki}b_{i}$
end loop (i)
$b_{k}=b_{k}/l_{kk}$
end loop (k)

Backward substitution

for k:=n stepdown until 1 do
for i:=k+1 step until n
$b_k=b_k-u_{ki}b_{i}$
end loop (i)
$x_{k}=b_{k}/u_{kk}$
end loop (k)

## Important Considerations

As with Gaussian elimination, LU decomposition is probably best used for relatively small, relatively non-sparse systems of equations (with small and non-sparse open to some interpretation). For larger and/or sparse problems, it would probably be best to either use an iterative method or use a direct solver package (e.g. DSCPACK) as opposed to writing one of your own.

If one has a single lefthand-side matrix and many right-hand side vectors, then LU decomposition would be a good solution procedure to consider. At the very least, it should be faster than solving each system separately with Gaussian elimination.

Here a Fortran code fragment for LU decomposition written by D. Partenov

do j = 1, n
do i = 1, j
if (i == 1) then
a(i,j) = a(i, j)
else
suma = 0.0
do k = 1, i-1
suma = suma + a(i,k)*a(k,j)
end do
a(i,j) = a(i,j) - suma
end if
end do
if ( j < n) then
do s = 1, n-j
i = j + s
if (j == 1) then
a(i,j) = a(i,j)/a(j,j)
else
suma = 0.0
do k = 1, j-1
suma = suma + a(i,k)*a(k,j)
end do
a(i,j) = (a(i,j) - suma)/a(j,j)
end if
end do
end if
end do
y(1) = b(1)
do i = 2, n
suma = 0.0
do j = 1, i-1
suma = suma + a(i,j)*y(j)
end do
y(i) = b(i) - suma
end do
i = n
j = n
x(i) = y(i)/a(i,j)
do s = 1, n -1
i = n - s
suma = 0.0
do j = i+1, n
suma = suma + a(i,j)*x(j)
end do
x(i) = (y(i) - suma)/a(i,i)
end do