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