[R] Basic Linear Algebra

[John Fox]

Dear Tom,

I'm not sure that I follow the question correctly, but I think that you want
to solve the over determined system tdata %*% x == sumtd.

Some time ago, I believe that I posted to the list a function for Gaussian
elimination that I use in teaching basic linear algebra. Here's an updated
version of that function, and its application to your problem:

GaussianElimination <- function(A, b, tol=sqrt(.Machine$double.eps),
verbose=FALSE, fractions=FALSE){
    # A: coefficient matrix
    # b: right-hand size vector
    # tol: tolerance for checking for 0 pivot
    # verbose: if TRUE, print intermediate steps
    # fractions: try to express nonintegers as rational numbers
    # If b is absent returns the reduced row-echelon form of A.
    # If b is present, reduces A to RREF carrying b along.
    if (fractions) {
        mass <- require(MASS)
        if (!mass) stop("fractions=TRUE needs MASS package")
        }
    if ((!is.matrix(A)) || (!is.numeric(A))) stop("argument must be a
numeric matrix")
    if ((!missing(b)) && ((!(length(b) == nrow(A))) || (!is.numeric(b))))
        stop("argument must be a numeric vector of length equal to the
number of columns in A")
    n <- nrow(A)
    m <- ncol(A)
    if (!missing(b)) A <- cbind(A, b)
    for (i in 1:min(c(m, n))){
        currentColumn <- A[,i]
        currentColumn[1:n < i] <- 0
        which <- which.max(abs(currentColumn))  # find maximum pivot in
current column at or below current row
        pivot <- A[which, i]
        if (abs(pivot) <= tol) next     # check for 0 pivot
        if (which > i) A[c(i, which),] <- A[c(which, i),]  # exchange rows
        A[i,] <- A[i,]/pivot            # pivot
        row <- A[i,]                    
        A <- A - outer(A[,i], row)      # sweep
        A[i,] <- row                    # restore current row
        if (verbose) if (fractions) print(fractions(A)) else print(round(A,
round(abs(log(tol,10)))))
        }
    for (i in 1:n) if (max(abs(A[i,1:m])) <= tol) A[c(i,n),] <- A[c(n,i),]
# 0 rows to bottom
    if (fractions) fractions (A) else round(A, round(abs(log(tol,10))))
    }
        
> GaussianElimination(tdata, sumtd)
                 b
[1,] 1 0 0 0 -3 30
[2,] 0 1 0 0 -1 12
[3,] 0 0 1 0 -2 20
[4,] 0 0 0 1  0 12

A caveat, however: Unlike some other members of this list, I'm not an expert
in numerical linear algebra and therefore can't vouch for the quality of the
algorithm -- it's meant basically for demonstration.

I hope this helps,
 John

--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
-------------------------------- 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Mulholland, Tom
> Sent: Friday, January 07, 2005 2:00 AM
> To: R-Help (E-mail)
> Subject: [R] Basic Linear Algebra
> 
> I don't normally have to go anywhere near this stuff , but it 
> seems to me that this should be a straight-forward process in R.
> 
> For the purposes of this enquiry I thought I would use 
> something I can work out on my own. 
> 
> So I have my matrix and the right hand results from that matrix
> 
> tdata <- 
> matrix(c(0,1,0,-1,-1,2,0,0,-5,-6,0,0,3,-5,-6,1,-1,-1,0,0),byro
> w = T,ncol = 5) sumtd <- c(0,0,0,-2)
> 
> > tdata
>      [,1] [,2] [,3] [,4] [,5]
> [1,]    0    1    0   -1   -1
> [2,]    2    0    0   -5   -6
> [3,]    0    0    3   -5   -6
> [4,]    1   -1   -1    0    0
> > sumtd
> [1]  0  0  0 -2
> > 
> 
> which I can calculate to give me 3x+30, x+12, 2x+20,  12x, x
> 
> Would someone be kind enough to point me in the right 
> direction as to which tools I should be using and any sage 
> words of advice for the barely informed.
> 
> Tom Mulholland
> 
>          _              
> platform i386-pc-mingw32
> arch     i386           
> os       mingw32        
> system   i386, mingw32  
> status                  
> major    2              
> minor    0.0            
> year     2004           
> month    10             
> day      04             
> language R
> 
> ______________________________________________
> 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