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