LU Decomposition using R

user_8384135 0
Viewed 4391

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)
linear-algebra math r matrix
1 Answers

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
user_47061710
Related