[R] row-echelon form (was no subject)

[Douglas Bates]

"Simon Fear" <Simon.Fear at synequanon.com> writes:

> I think one needs an LU decomposition rather than QR. 
> However, I couldn't find anything off the shelf to do 
> an LU, other than learning that determinant() now 
> uses LU instead of QR or SVD, so the code to do it must
> be in there for those that want it.

I was about to write that the version of the Matrix package for
r-devel (to be R-1.9.0) has an LU decomposition but then I checked and
found that, although the underlying C code is there, I haven't written
an accessor function.  In any case the LU decomposition is computed
and stored whenever it is needed, say for calculating a condition
number or for solving a linear system.

> library(Matrix)   # version 0.7-3 for R-1.9.0
> mm = Matrix(rnorm(9), nc = 3)
> mm
           [,1]       [,2]       [,3]
[1,] -0.4186151 -0.3959071 -0.1717103
[2,] -0.5334526 -2.2817253 -1.5398915
[3,] -0.5993042 -1.3396313  0.2062784
> rcond(mm)
[1] 0.04505896
> str(mm)
 list()
 - attr(*, "x")= num [1:9] -0.419 -0.533 -0.599 -0.396 -2.282 ...
 - attr(*, "Dim")= int [1:2] 3 3
 - attr(*, "rcond")= Named num 0.0451
  ..- attr(*, "names")= chr "O"
 - attr(*, "factorization")=List of 1
  ..$ LU: list()
  .. ..- attr(*, "x")= num [1:9] -0.599  0.890  0.699 -1.340 -1.089 ...
  .. ..- attr(*, "Dim")= int [1:2] 3 3
  .. ..- attr(*, "pivot")= int [1:3] 3 2 3
  .. ..- attr(*, "class")= atomic [1:1] LU
  .. .. ..- attr(*, "package")= chr "Matrix"
 - attr(*, "class")= atomic [1:1] geMatrix
  ..- attr(*, "package")= chr "Matrix"
> expand(mm at factorization$LU)
$L
          [,1]       [,2] [,3]
[1,] 1.0000000  0.0000000    0
[2,] 0.8901200  1.0000000    0
[3,] 0.6985019 -0.4955766    1

$U
           [,1]      [,2]       [,3]
[1,] -0.5993042 -1.339631  0.2062784
[2,]  0.0000000 -1.089293 -1.7235040
[3,]  0.0000000  0.000000 -1.1699244

As I understand it the LU decomposition is much more widely used in
numerical linear algebra than is the row-reduced echelon form.