LU Decomposition using R

Question

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)

Not reproducible. What package(s) are you using? `library` is missing. Edit your request. — Bhas, Sep 14 '18 at 05:54
There are not many LU factorization R questions on Stack Overflow BTW. As a pointer to all readers, the following two are rather educational: 1> [LU decomposition with row pivot](https://stackoverflow.com/q/3642188/4891738); 2> [Write a trackable R function that mimics LAPACK's dgetrf for LU factorization](https://stackoverflow.com/q/51687808/4891738). — Zheyuan Li, Sep 14 '18 at 18:07

score 7 · Answer 1 · answered Sep 14 '18 at 06:35

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

LU Decomposition using R

1 Answers1