[R] Solutions of equation systems

[Ravi Varadhan]

Since you have an under-determined system of linear equations, you have infinitely many solutions.  One reasonable solution is the "minimum-norm" solution.  You can obtain this from singular-value decomposition as follows:

# Minimum norm solution `x':  A %*% x = b, such that ||x|| is minumum

d <- svd(A)

x.min <- c(d$v %*% diag(1/d$d, length(d$d)) %*% t(d$u) %*% b)  # min-norm solution

Hope this helps,
Ravi.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Berend Hasselman <bhh at xs4all.nl>
Date: Saturday, August 15, 2009 10:03 am
Subject: Re: [R] Solutions of equation systems
To: r-help at r-project.org


>  
>  
>  Moreno Mancosu wrote:
>  > 
>  > Dear Greg,
>  > 
>  > I tried to use function "solve", but I have some problems with it.
>  > Let's try with a simple equation:
>  > 
>  > x + 3y -2z = 5
>  > 3x + 5y + 6z =7
>  > 
>  >  > a<-matrix(c(1,3,-2,3,5,6),nrow=2,ncol=3)
>  >  > a
>  >      [,1] [,2] [,3]
>  > [1,]    1   -2    5
>  > [2,]    3    3    6
>  >  > b<-matrix(c(5,7),nrow=2,ncol=1)
>  >  > b
>  >      [,1]
>  > [1,]    5
>  > [2,]    7
>  > 
>  > When I try to solve this system, I find this error:
>  > 
>  >  > solve(a,b)
>  > Error in solve.default(a, b) : 'b' must be compatible with 'a'
>  > 
>  > Also, when i try to solve a system with n equation and n variables, 
> 
>  > solve works perfectly.
>  > 
>  
>  If you do ?solve then you will see that 
>  the matrix a must be non square (it says so in the description of a in
>  section Arguments)/
>  and in section Details right at the end it is written that "qr.solve 
> can
>  handle non-square systems".
>  
>  Your example can be solved in several ways.
>  See this
>  
>  <example>
>  A <- matrix(c(1,2,3,5,6,7),nrow=2,byrow=T)
>  b <- c(2,8) 
>  
>  A <- matrix(c(1,3,-2,3,5,6),nrow=2,ncol=3) 
>  b <- matrix(c(5,7),nrow=2,ncol=1)
>  
>  # a particular solution
>  xsolve <- qr.solve(A,b)
>  xsolve
>  
>  # another particular solution
>  
>  xany <- solve(qr(A,LAPACK=T),b)   
>  xany
>  
>  # QR decomposition of A-transposed
>  
>  ATQR <- qr(t(A),LAPACK=TRUE)
>  ATQR
>  
>  qr.Q(ATQR)
>  qr.R(ATQR)
>  
>  # minimum norm solution comes from QR decomposition of t(A)
>  # and solving t(R) %*% t(Q) %*% x = b in two steps
>  # ! Lapack does a pivoted QR
>  
>  z <- solve(t(qr.R(ATQR)),b[ATQR$pivot])
>  xmin <- qr.Q(ATQR) %*% z
>  xmin
>  
>  # null space of A  
>  library(MASS)
>  A.null <- Null(t(A))
>  A.null
>  
>  # test
>  A %*% (xmin+5*A.null)  
>  </example>
>  
>  ! Use the manual
>  
>  Berend
>  -- 
>  View this message in context: 
>  Sent from the R help mailing list archive at Nabble.com.
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.