[Rd] svd() (Linpack) problems/bug for ill-conditioned matrices (PR#594)
maechler@stat.math.ethz.ch
maechler@stat.math.ethz.ch
Wed, 5 Jul 2000 18:31:45 +0200 (MET DST)
After fixing princomp(), recently,
{tiny negative eigen-values are possible for non-negative
definite matrices}
Fritz Leisch drew my attention to the fact the not only eigen() can be
funny, but also svd().
Adrian Trappleti found that the singular values returned
can be "-0" instead of "0". This will be a problem in something like
sd <- svd(Mat) $ d
max(sd) / min(sd)
which gives -Inf instead of Inf in the above case.
I did some more searching and testing with svd(), could easily
reproduce the "-0" problem,
hilbert <- function(n, k=n) { i <- 1:n; 1 / outer(i - 1, i[1:k], "+") }
M1 <- cbind(hilbert(12,8), pi*1e-15, 0)
(sd <- svd(M1, nu=0, nv=0) $ d)
sd[1] / sd[10] #--> -Inf
but even worse: On both Solaris (2.5.1, gcc 2.95.1,
f77 "WorkShop Compilers 4.2"
and Linux (Redhat 6.2, enhanced with gcc 2.95.2;
g77 version egcs-2.91.66 19990314/Linux
(egcs-1.1.2 release)
(from FSF-g77 version 0.5.24-19981002)
I have managed to make svd() stop with
something like "error 170 in dsvdc" which basically comes from the
underlying LINPACK routine DSVDC() which did not converge in 30 QR iterations.
example code appended
--- The examples I've found are machine dependent and disjoint for Linux &
Solaris.
HOWEVER: S-PLUS 5.1 on both machines does NOT fail for the identical
matrices!
S-PLUS has basically .Fortran("dsvdcs" ,
and looks like it is calling Linpack as well, but note that the subroutine
has an extra "s" appended (``S version of Linpack'' ??).
---
Then I had a short look at GSL (GNU scientific library)
http://sourceware.cygnus.com/gsl/
which calls a LAPACK (instead of LINPACK) routine for svd which seems to
use Householder- instead of QR- iterations...
The whole makes me think that we should really shift away from Linpack and
towards Lapack for R's numerical linear algebra....
Here is my code which gives both an error (but not the same!) for R
on Linux and Solaris and all runs through S-PLUS 5.1 on both machines.
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
## "better":
hilbert <- function(n, k=n) { i <- 1:n; 1 / outer(i - 1, i[1:k], "+") }
M1 <- cbind(hilbert(12,8), 3.14e-15, 0)
M2 <- cbind(hilbert(15,8), .2,0)
H2c <- hilbert(200,180)
## Gives SVD error "error 170 in dsvdc" on (MM's) Solaris :
ii.sol <- c(87, 40, 120, 164, 93, 98, 108, 156, 77, 52, 127, 50, 58, 123,
17, 147, 10, 152, 103, 136, 20, 34, 6, 30, 44, 130, 47, 35, 168,
4, 115, 3, 29, 71, 56, 128, 90, 18, 60, 25, 125, 74, 91, 16,
9, 111, 46, 137, 53, 110, 21, 94, 78, 104, 100, 99, 76, 126,
62, 61, 117, 54, 140, 66, 150, 19, 149, 114, 63, 139, 2, 166,
135, 23, 107, 8, 158, 39, 69, 75, 92, 112, 86, 124, 13, 118,
151, 80, 146, 68, 43, 159, 36, 31, 144, 85, 122, 141, 11, 167,
157, 1, 132, 97, 133, 79, 138, 49, 45, 72, 119, 84, 59, 131,
5, 22, 116, 154, 142, 48, 73, 67, 161, 7, 155, 42, 41, 143, 57,
33, 162, 102, 134, 101, 37, 148, 64, 160, 24, 27, 55, 83, 121,
38, 106, 113, 28, 109, 70, 14, 32, 12, 65, 96, 129, 89, 163,
88, 105, 15, 95, 145, 81, 165, 51, 153, 82, 169, 26)
H2c.sol <- cbind(H2c[,ii.sol], 1)
## NEW rec: -9.413e+14 ; nc = 164
## Gives SVD error 165 on Linux:
ii <- c(34, 101, 120, 53, 140, 108, 112, 9, 10, 96, 154, 89, 17, 131,
130, 123, 126, 51, 95, 85, 147, 87, 132, 60, 58, 64, 150, 139,
116, 36, 78, 107, 118, 104, 127, 61, 82, 84, 63, 23, 114, 124,
83, 162, 55, 99, 50, 141, 8, 24, 25, 19, 65, 15, 2, 143, 145,
148, 138, 136, 57, 22, 26, 1, 152, 129, 21, 69, 59, 20, 163,
98, 92, 106, 88, 56, 70, 117, 37, 48, 115, 79, 68, 142, 4, 122,
76, 28, 155, 29, 45, 13, 30, 47, 62, 31, 121, 119, 125, 74, 144,
27, 93, 100, 71, 40, 149, 16, 157, 41, 33, 113, 133, 54, 14,
156, 11, 97, 72, 75, 77, 49, 111, 128, 18, 35, 42, 6, 160, 32,
3, 91, 66, 73, 44, 39, 81, 90, 153, 135, 43, 151, 159, 102, 161,
12, 110, 5, 164, 67, 80, 46, 94, 7, 105, 134, 38, 146, 86, 103,
137, 52, 158, 109)
H2c.lin <- cbind(H2c[,ii], 1)
(sd <- svd(M1, nu=0, nv=0) $ d) ; sd[1] / sd[length(sd)] #--> -Inf
(sd <- svd(M2, nu=0, nv=0) $ d) ; sd[1] / sd[length(sd)] #--> -Inf
summary(svd(H2c.sol)$d)## Error on Solaris only
summary(svd(H2c.lin)$d)## Error on Linux only (RH 6.2, gcc 2.95.2)
## BTW: quite a strong negative ratio
range( ev <- eigen(var(H2c.lin), only=T)$val )
ev[1] / ev[length(ev)] # -5.868e+15 (solaris); -1.2115e+16 (linux)
### Try different ii and always do the following line :
summary(svd(cbind(H2c[,ii], 1))$d)
## BTW: Negative Eigenvalues
range( eigen(var(M2), only=T)$val )
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._