[R] LU bug in Matrix package

[Thomas Lumley]

On Thu, 28 Dec 2006, yangguoyi.ou wrote:

> There is a bug in Matrix package, please check it, thanks!

You didn't say what R code you ran to get that output or why you think it 
is wrong.

Let us experiment to see if we can guess what the problem is from the 
limited information provided
> x<-t(Matrix(1:25,5))
> x
5 x 5 Matrix of class "dgeMatrix"
      [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
[5,]   21   22   23   24   25

> lux<-lu(x)
Warning message:
Exact singularity detected during LU decomposition.

Now check that the decomposition is correct
> with(expand(a), P%*%L%*%U)
5 x 5 Matrix of class "dgeMatrix"
      [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
[5,]   21   22   23   24   25

Ok, it is correct.
Now let's Look at the matrices

> expand(a)
$L
5 x 5 Matrix of class "dtrMatrix"
      [,1]       [,2]       [,3]       [,4]       [,5]
[1,] 1.00000000          .          .          .          .
[2,] 0.04761905 1.00000000          .          .          .
[3,] 0.52380952 0.50000000 1.00000000          .          .
[4,] 0.28571429 0.75000000 0.35400130 1.00000000          .
[5,] 0.76190476 0.25000000 0.50000000 0.00000000 1.00000000

$U
5 x 5 Matrix of class "dtrMatrix"
      [,1]          [,2]          [,3]          [,4]          [,5]
[1,]  2.100000e+01  2.200000e+01  2.300000e+01  2.400000e+01  2.500000e+01
[2,]             .  9.523810e-01  1.904762e+00  2.857143e+00  3.809524e+00
[3,]             .             . -1.919092e-15 -3.711332e-15 -5.614541e-15
[4,]             .             .             .  3.781500e-16  6.288326e-16
[5,]             .             .             .             .  0.000000e+00

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

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


Perhaps the output was a trimmed version of this? The L and U matrices 
agree, but you didn't show the permutation matrix.

Did you just miss the fact that the decomposition is P%*%L%*%U? If you 
aren't used to S4 classes it is easy to miss the fact that class?LU 
gives help on the structure of the "LU" class, and if you try to guess 
what the components are without the documentation it is easy to be 
confused.


 	-thomas




>
>
>
>
> Matlab result:
>
>
>
> x =
>
>     1     2     3     4     5
>
>     6     7     8     9    10
>
>    11    12    13    14    15
>
>    16    17    18    19    20
>
>    21    22    23    24    25
>
>
>
>>> lu(x)
>
>
>
> ans =
>
>
>
>   21.0000   22.0000   23.0000   24.0000   25.0000
>
>    0.0476    0.9524    1.9048    2.8571    3.8095
>
>    0.7619    0.2500   -0.0000   -0.0000   -0.0000
>
>    0.5238    0.5000    0.4839         0    0.0000
>
>    0.2857    0.7500   -0.0645         0    0.0000
>
>
>
>
>
> Gsl result:
>
>
>
> 21 22 23 24 25
>
> 0.047619 0.952381 1.90476 2.85714 3.80952
>
> 0.761905 0.25 -3.44169e-015 -6.88338e-015 -1.03251e-014
>
> 0.52381 0.5 -0.0322581 -1.77636e-015 -1.66533e-015
>
> 0.285714 0.75 -0.0645161 -0 2.22045e-016
>
>
>
>
>
> R result:
>
> $L
>
> 5 x 5 Matrix of class "dtrMatrix"
>
>     [,1]       [,2]       [,3]       [,4]       [,5]
>
> [1,] 1.00000000          .          .          .          .
>
> [2,] 0.04761905 1.00000000          .          .          .
>
> [3,] 0.52380952 0.50000000 1.00000000          .          .
>
> [4,] 0.28571429 0.75000000 0.35400130 1.00000000          .
>
> [5,] 0.76190476 0.25000000 0.50000000 0.00000000 1.00000000
>
>
>
> $U
>
> 5 x 5 Matrix of class "dtrMatrix"
>
>     [,1]          [,2]          [,3]          [,4]          [,5]
>
> [1,]  2.100000e+01  2.200000e+01  2.300000e+01  2.400000e+01  2.500000e+01
>
> [2,]             .  9.523810e-01  1.904762e+00  2.857143e+00  3.809524e+00
>
> [3,]             .             . -1.919092e-15 -3.711332e-15 -5.614541e-15
>
> [4,]             .             .             .  3.781500e-16  6.288326e-16
>
> [5,]             .             .             .             .  0.000000e+00
>
>
>
>
>
>
>
>
>
>
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle