Linear equation solvers by LU decomposition

The LU decomposition is a factoring of an arbitrary square matrix into an upper-triangular and lower-triangular pair of factors. This factoring is very useful for solving linear equations, since the inversion of such upper/lower matrices is trivial. We select a standard form for the lower triangular matrix; ones on its diagonal. For example

$\left(\begin{array}{ccc}2&-1&-2\\ -4&6&3\\ -4&-2&8\end{array}\right)=\left(\begin{array}{ccc}1&0&0\\ -2&1&0\\ -2&-1&1\end{array}\right)\left(\begin{array}{ccc}2&-1&-2\\ 0&4&-1\\ 0&0&3\end{array}\right)$

LU decomposition can be done very efficiently, it is significantly faster and easier to program than QR decomposition, and it is technically much simpler.
$\bullet$ Start with

$LU=\left(\begin{array}{ccc}1&0&0\\ L_{21}&1&0\\ L_{31}&L_{32}&1\end{array}\right)\left(\begin{array}{ccc}U_{11}&U_{12}&U_{13}\\ 0&U_{22}&U_{23}\\ 0&0&U_{33}\end{array}\right)$

$=\left(\begin{array}{ccc}m_{11}&m_{12}&m_{13}\\ m_{21}&m_{22}&m_{23}\\ m_{31}&m_{32}&m_{33}\end{array}\right)=M$

Dotting the first row of $L$ into each column of $U$ we see

$U_{11}=m_{11}, \quad U_{12}=m_{12}\quad\cdots\quad U_{1n}=m_{1n}$

$\bullet$ Now dot the second row of $L$ into each column of $U$ and solve

$L_{21}U_{11}=m_{11}, \quad L_{21}U_{12}+U_{22}=m_{22}, \quad L_{21}U_{13}+U_{23}=m_{23}$

Solve the first equation, so $L_{21}=m_{11}/U_{11}$ is known, now solve the remaining equations for $U_{2 n}$ :

$U_{22}=m_{22}- L_{21}U_{12}, \qquad U_{23}=m_{23}-L_{21}U_{13}$

so the second rows of both $L$ and $U$ are gotten from $M$ and the first row of $U$
You can see that if $m_{11}=0$ , the process fails at the very first stage, and so we usually perform a permutation or pivot $P$ that rearranges the rows of of $M$ so that $PM$ has the property that $|(PM)_{kk}|>|(PM)_{k+i,k}|$ for all $i$ . In other words the diagonal entries of $PM$ are larger in absolute value than any entries below them in the same column. We apply the LU algorithm to the permuted matrix

$LU=PM$

# in R
Pivot <- function(M){
  N <- ncol(M)
  P <- matrix(rep(0,N^2),nrow=N,ncol=N)
  for (k in 1:N){P[k,k]<- 1.0}
    Mc <- M
    for (k in 1:(N-1)){
      largest <- which.max(abs(Mc[k:N,k]))
        if(largest>1){
          Mc[c(k,largest+k-1),] <- Mc[c(largest+k-1,k),]
          P[c(k,largest+k-1),] <- P[c(largest+k-1,k),]}     
    }
    return (list("perm"=P,"matrix"=Mc))
}

$\bullet$ Finally dot the last row of $L$ into $U$ ;

$L_{31}U_{11}=m_{31}, \quad L_{31}U_{12}+L_{32}U_{22}=m_{32}, \quad L_{31}U_{13}+L_{32}U_{23}+U_{33}=m_{33}$

Solve the first, now $L_{31}=m_{31}/U_{11}$ is known, use this to solve the second

$L_{32}={m_{32}-L_{31}U_{12}\over U_{22}}$

so now $L_{32}$ known, and the third equation contains one unknown

$U_{33}=m_{33}-\Big(L_{31}U_{13}+L_{32}U_{23}\Big)$

Some computations (for example the first row of $U$ ) can be saved by copying $M$ into $U$ immediately and creating the rest of $U$ right on top of $M$ , destroying $M$ in the process (construction in situ)

# this R code returns L,U,P
U <- function(M){
  N <- ncol(M)
  L <- matrix(rep(0,N^2),nrow=N,ncol=N)
  U <- matrix(rep(0,N^2),nrow=N,ncol=N)
  Q <- Pivot(M)
  for (i in 1:N){
    for (k in i:N){
      sm <- -sum(L[i,1:i]*U[1:i,k])
      U[i,k] <- Q$matrix[i,k]+sm}
    for (k in i:N){
      if(i==k){L[i,i] <- 1.0}
      else{sm <- -sum(L[k,1:i]*U[1:i,i])
        L[k,i]<- (Q$matrix[k,i]+sm)/U[i,i]}
    }
  }
     return (list("l"=L,"u"=U,"p"=t(Q$perm)))
}

Inverting triangular linear systems
Back substitution

$\left(\begin{array}{cccc} u_{11}& u_{12} & \cdots & u_{1n}\\ 0 & u_{22} & \cdots & u_{2n}\\ \vdots & \vdots & \ddots &\vdots\\ 0 & \cdots & 0 & u_{nn}\end{array}\right)\left(\begin{array}{c} x_1\\ x_2\\ \vdots\\ x_n\end{array}\right)= \left(\begin{array}{c} b_1\\ b_2\\ \vdots\\ b_n\end{array}\right)$

(1) $\begin{eqnarray*} u_{nn}x_n & = & b_n, \quad x_n={b_n\over u_{nn}}\nonumber\\ u_{n-1,n-1}x_{n-1} +u_{n-1,n}x_n& = & b_{n-1}, \quad x_n={b_{n-1}-u_{n-1,n}b_n\over u_{n-1,n-1}}\nonumber\end{eqnarray*}$

so for the $j^{th}$ row element

$x_j={b_j-\sum_{k=j+1}^n u_{j,k}x_k\over u_{jj}}$

Forward substitution

$\left(\begin{array}{cccc} \ell_{11}&02} & \cdots & 0\\ \ell_{21} & \ell_{22} & \cdots & 0\\ \vdots & \vdots & \ddots &\vdots\\ \ell_{n1} & \ell_{n2} & \cdots & \ell_{nn}\end{array}\right)\left(\begin{array}{c} x_1\\ x_2\\ \vdots\\ x_n\end{array}\right)= \left(\begin{array}{c} b_1\\ b_2\\ \vdots\\ b_n\end{array}\right)$

(2) $\begin{eqnarray*} \ell_{11}x_1 & = & b_1, \quad x_1={b_1\over \ell_{11}}\nonumber\\ \ell_{21}x_1 +\ell_{22}x_2& = & b_2, \quad x_2={b_2-\ell_{21}x_1\over \ell_{22}}\nonumber\end{eqnarray*}$

$x_j={b_j-\sum_{k=1}^{j-1}\ell_{jk}x_k\over \ell_{jj}}$

Fsub <- function(L,b){
    N <- ncol(L)
    y <-  matrix(rep(0,N),nrow=N,ncol=1)
    for (j in 1:N){
      y[j]<- (b[j]-sum(L[j,1:(j-1)]*y[1:j-1]))/L[j,j]}
    return(y)
    }

To write a solver for $M\cdot \bm{x}=\bm{b}$ ,
$\bullet$ apply the LU decomposition

$PM=LU, \qquad PM\cdot\bm{x}=P\cdot\bm{b}=LU\cdot\bm{x}$

$\bullet$ solve $L\cdot\bm{y}=P\cdot\bm{b}$ for $\bm{y}$ , then
$\bullet$ solve $U\cdot\bm{x}=\bm{y}$ for $\bm{x}$ .

You probably can’t write any linear algebra code that can compete with GEMM or the solvers in LAPACK and BLAS when it comes to speed, but the point of “rolling your own” is to learn the algorithms, and to create self-contained codes that don’t require numerous calls to external libraries.

Home 2.0

Linear equation solvers by LU decomposition

Related Posts

Linear equation solvers by LU decomposition

Published by Jeff on August 7, 2021August 7, 2021

Related Posts

Protected: AI/ML-1. Neural nets in Julia

Algorithms: Tridiagonal Hessenberg form using Krylov subspaces

Algorithms: Devising quadrature (integration) formulas