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

[John Fox]

Dear Doug and Simon,

The U matrix in the LU decomposition is upper triangular, but (as I
understand it) it's not in row-echelon because the leading entries aren't 1
and the columns with leading entries aren't fully swept.

I'm certainly not making any claims, other than pedagogical, for the REF. In
fact, the function that I sent originated as a demonstration for a class. By
the way, I think that setting the tolerance in the function to
Machine$double.eps (which is what I did) is too small: Gabor Grothendieck
sent me an example in which my function failed with the default tolerance.

Regards,
 John

> -----Original Message-----
> From: Douglas Bates [mailto:bates at bates4.stat.wisc.edu] On 
> Behalf Of Douglas Bates
> Sent: Friday, March 05, 2004 8:26 AM
> To: Simon Fear
> Cc: John Fox; Spencer Graves; r-help at stat.math.ethz.ch
> Subject: Re: [R] row-echelon form (was no subject)
> 
> "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.