Problem 2: Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.

Here is my attempt for a 3x3 matrix…

LU decomposition from lecture notes - my LU

x<-matrix(c(1,-2,3,4,8,4,-3,5,7),nrow=3, ncol=3)
x
##      [,1] [,2] [,3]
## [1,]    1    4   -3
## [2,]   -2    8    5
## [3,]    3    4    7

Find the LU decomposition of the above matrix x

MyLU<-function(x) {
E21=matrix(c(1,-x[2,1]/x[1,1],0,0,1,0,0,0,1),nrow=3)
b<-E21 %*% x
E31=matrix(c(1,0,-b[3,1]/b[1,1],0,1,0,0,0,1),nrow=3)
c<-E31 %*% E21 %*% x
E32=matrix(c(1,0,0,0,1,-c[3,2]/c[2,2],0,0,1),nrow=3)
U=E32%*%E31 %*% E21 %*% x
U
L<-solve(E21) %*% solve(E31) %*% solve(E32)
L%*%U
(L%*%U==x)
result <- list( L=L, U=U )
    return( result )
}

MyLU(x)
## $L
##      [,1] [,2] [,3]
## [1,]    1  0.0    0
## [2,]   -2  1.0    0
## [3,]    3 -0.5    1
## 
## $U
##      [,1] [,2] [,3]
## [1,]    1    4 -3.0
## [2,]    0   16 -1.0
## [3,]    0    0 15.5

This code was found online -https://rdrr.io/cran/matrixcalc/src/R/lu.decomposition.R

with minor modifications.

This is more advanced and can accommodate any square matrix. It is more sophisticated and better than my code above. It is more versatile in accepting any size matrix.

Find the LU decomposition of above matrix X

x<-matrix(c(1,-2,3,4,8,4,-3,5,7),nrow=3, ncol=3)
LU<-function(x) {

if ( nrow(x)!=ncol(x))
       stop( "No LU decomp - not a square matrix" )
  
    n <- nrow( x )
    L <- matrix( 0, nrow=n, ncol=n )
    U <- matrix( 0, nrow=n, ncol=n )
    diag( L ) <- rep( 1, n )
    for ( i in 1:n ) {
        ip1 <- i + 1
        im1 <- i - 1
        for ( j in 1:n ) {
            U[i,j] <- x[i,j]
            if ( im1 > 0 ) {
                for ( k in 1:im1 ) {
                    U[i,j] <- U[i,j] - L[i,k] * U[k,j]
                }
            }
        }
        if ( ip1 <= n ) {
            for ( j in ip1:n ) {
                L[j,i] <- x[j,i]
                if ( im1 > 0 ) {
                    for ( k in 1:im1 ) {
                        L[j,i] <- L[j,i] - L[j,k] * U[k,i]
                    }
                }
                if ( U[i,i] == 0 )
                    stop( "argument x is a singular matrix" )
                L[j,i] <- L[j,i] / U[i,i]
            }    
        }
    }
    result <- list( L=L, U=U )
    return( result )
}
LU(x)
## $L
##      [,1] [,2] [,3]
## [1,]    1  0.0    0
## [2,]   -2  1.0    0
## [3,]    3 -0.5    1
## 
## $U
##      [,1] [,2] [,3]
## [1,]    1    4 -3.0
## [2,]    0   16 -1.0
## [3,]    0    0 15.5