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
In-Reply-To: <1361602670-885614@pacwan.mm-soft.fr>

>>>>> "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
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-