6

The following function does not use row pivoting for LU decomposition. Is there an existing function in R that does LU decomposition with row pivot?

> require(Matrix)
> expand(lu(matrix(rnorm(16),4,4)))
$L
4 x 4 Matrix of class "dtrMatrix"
     [,1]        [,2]        [,3]        [,4]       
[1,]  1.00000000           .           .           .
[2,]  0.13812836  1.00000000           .           .
[3,]  0.27704442  0.39877260  1.00000000           .
[4,] -0.08512341 -0.24699820  0.04347201  1.00000000

$U
4 x 4 Matrix of class "dtrMatrix"
     [,1]       [,2]       [,3]       [,4]      
[1,]  1.5759031 -0.2074224 -1.5334082 -0.5959756
[2,]          . -1.3096874 -0.6301727  1.1953838
[3,]          .          .  1.6316292  0.6256619
[4,]          .          .          .  0.8078140

$P
4 x 4 sparse Matrix of class "pMatrix"

[1,] | . . .
[2,] . | . .
[3,] . . . |
[4,] . . | .
user236215
  • 7,278
  • 23
  • 59
  • 87
  • Share an educational Q & A on LU factorization: [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 17:39

2 Answers2

8

The lu function in R is using partial (row) pivoting. You did not give the original matrix with your example, so I will create a new example to demonstrate.

Function lu in R is computing A = PLU, which is equivalent to computing the LU decomposition of matrix A with its rows permuted by the permutation matrix P-1: P-1A = LU. See the Matrix package documentation for more information.

Example

> A <- matrix(c(1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1), 4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1   -1   -1
[3,]    1   -1   -1    1
[4,]    1   -1    1   -1

Here's the L factor:

> luDec <- lu(A)
> L <- expand(luDec)$L
> L
4 x 4 Matrix of class "dtrMatrix" (unitriangular)
     [,1] [,2] [,3] [,4]
[1,]    1    .    .    .
[2,]    1    1    .    .
[3,]    1    0    1    .
[4,]    1    1   -1    1

Here's the U factor:

> U <- expand(luDec)$U
> U
4 x 4 Matrix of class "dtrMatrix"
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    .   -2   -2    0
[3,]    .    .   -2   -2
[4,]    .    .    .   -4

Here's the permutation matrix:

> P <- expand(luDec)$P
> P
4 x 4 sparse Matrix of class "pMatrix"

[1,] | . . .
[2,] . . | .
[3,] . | . .
[4,] . . . |

We can see that LU is a row-permuted version of A:

> L %*% U
4 x 4 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1   -1   -1    1
[3,]    1    1   -1   -1
[4,]    1   -1    1   -1

Going back to the original identity A = PLU we can recover A (compare with A above):

> P %*% L %*% U
4 x 4 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1   -1   -1
[3,]    1   -1   -1    1
[4,]    1   -1    1   -1
989
  • 12,579
  • 5
  • 31
  • 53
David Alber
  • 17,624
  • 6
  • 65
  • 71
  • 1
    + 1 David. Excellent answer. I was wondering how we interpret the permutation matrix P. What do the bars and dots mean in the respective positions mean? – ColorStatistics Jun 19 '21 at 11:34
1

Perhaps this does the job. However, there is not a windows binary and I can't try it.

gd047
  • 29,749
  • 18
  • 107
  • 146
  • Magma requires the underlying magma library as well as an NVidia card for GPU computing. There are pivoting solutions in R -- I'd ask on r-help. – Dirk Eddelbuettel Sep 05 '10 at 13:23