# Re: R-alpha: svd

Douglas Bates (bates@stat.wisc.edu)
Sat, 14 Dec 1996 17:48:26 -0600 (CST)

Message-Id: <m0vZ3og-000hi2C@franz.stat.wisc.edu>
Date: Sat, 14 Dec 1996 17:48:26 -0600 (CST)
From: Douglas Bates <bates@stat.wisc.edu>
To: Ellis / Gilbert <la-jassine@aix.pacwan.net>
Subject: Re: R-alpha: svd

>>>>> "Ellis" == Ellis / Gilbert <la-jassine@aix.pacwan.net> writes:

Ellis> In S, as I recall, the inverse of an arbitrary square matrix
Ellis> X can be found by
v <- svd(X)
inv.X <- v\$u %*% diag(v\$d) %*% t(v\$v)

Actually that is X, not the inverse of X.

Ellis> but in R it appears to be
inv.X <- v\$v %*% diag(v\$d) %*% t(v\$u)

Not quite.  To get the inverse you would need to take the inverse of
the singular values as in
inv.X <- v\$v %*% diag(1/v\$d) %*% t(v\$u)

As pointed out on this list, the simplest way of getting an inverse of
a square matrix is as solve(X) which uses a QR decomposition, not a
singular value decomposition.  A numerical analyst will tell you that
you should look very carefully at the computation before you decide
that you need to take an inverse of a matrix in the first place.  Most
times you don't need it; you are just used to writing formulas in
terms of inverses and matrix multiplications.  I cringe every time I
see code for least squares computations that features something like
beta <- solve(t(X) %*% X) %*% t(X) %*% y
Besides doing much more work than necessary, it is a less stable way
of doing the calculation than
beta <- qr.beta(qr(X), y)

A week or two ago there were several messages over S-news on the more
complicated question of computing the generalized inverse of a matrix
of arbitrary dimension.  That computation used a formula like the one
above but compressed a bit.  For example it is a waste to create a
diagonal matrix to be used in a matrix multiplication.  You can write
the same result as
inv.X <- v\$v %*% (t(v\$u) * 1/v\$d)

The next issue is what do you do about singular values that are near
zero.  For that I would suggest that you consult the archives of the
S-news list to see why the final form of the function is
ginv <-
function (X, tol = sqrt(.Machine\$double.eps))
## Generalized Inverse of a Matrix
{
svdX <- svd(X)
nonZero <- svdX\$d > tol * svdX\$d[1]
svdX\$v[ , nonZero] %*% (t(svdX\$u[, nonZero]) * 1/svdX\$d[nonZero])
}
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-testers mailing list -- For info or help, send "info" or "help",
To [un]subscribe, send "[un]subscribe"
(in the "body", not the subject !)  To: r-testers-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-