CODEWITHSUNDEEP
rmathmatrixlinear-algebranumerical-methods
stat77
stat77

Reputation: 113

LU Decomposition using R

I am trying to run an LU decomposition using R. Here is the reproducible code. I am not understanding why my permutation matrix is different from the solution. The L and U matrices are correct. But for the permutation matrix, the 1st and the 2nd rows and 3rd and the 4th rows are interchanged. Hence, I am not getting the correct solution to the system of linear equations. Would appreciate your help.

install.packages("Matrix")
library(Matrix)
(A <- matrix(c(4, 3, -2, 5, 2, -4, 6, 1, -1, 2, -5, 6, 3, 5, -2, -3), nrow = 4))
(B <- matrix(c(16.9, -14, 25, 9.4), nrow = 4))

luA <- lu(A)
elu <- expand(luA)
(L <- elu$L)
(U <- elu$U)
(P <- elu$P)

(Y <- solve(L) %*% P %*% B)
(X <- solve(U) %*% Y)

Upvotes: 4

Views: 6903

Answers (1)

Sandipan Dey
Sandipan Dey

Reputation: 23099

With R implementation, it seems that we have A = PLU (instead of PA = LU). Hence, the following works:

all.equal(Matrix(A), with(elu, P %*% L %*% U))
# TRUE

(Y <- solve(L, solve(P) %*% B)) # solve LY = inv(P).B instead of LY = PB
(X <- solve(U, Y))
 
X
#4 x 1 Matrix of class "dgeMatrix"
#     [,1]
#[1,]  4.5
#[2,]  1.6
#[3,] -3.8
#[4,] -2.7 

all.equal(X, Matrix(solve(A, B)))
# TRUE

For double-checking, we can compute LU factorization with Gaussian elimination, for the given A and B:

# Gaussian elimination steps
E21 <- matrix(c(1,0,0,0,-3/4,1,0,0,0,0,1,0,0,0,0,1), nrow=4, byrow=T)
E31 <- matrix(c(1,0,0,0,0,1,0,0,1/2,0,1,0,0,0,0,1), nrow=4, byrow=T)
E41 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,-5/4,0,0,1), nrow=4, byrow=T)
E32 <- matrix(c(1,0,0,0,0,1,0,0,0,14/11,1,0,0,0,0,1), nrow=4, byrow=T)
E42 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,-3/11,0,1), nrow=4, byrow=T)
E43 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,13/4,1), nrow=4, byrow=T)

U <- E43 %*% E42 %*% E32 %*% E41 %*% E31 %*% E21 %*% A
L <- solve(E43 %*% E42 %*% E32 %*% E41 %*% E31 %*% E21)

all.equal(L %*% U, A)
# [1] TRUE

all.equal(solve(U, solve(L, B)), solve(A, B))
# [1] TRUE

Here the permutation matrix P is identity (no row swapping needed).

Upvotes: 8

Related Questions