[R] Solutions of equation systems
Ravi Varadhan
rvaradhan at jhmi.edu
Sat Aug 15 18:29:21 CEST 2009
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.
More information about the R-help
mailing list