R devel - Jul 2000 - svd() (Linpack) problems/bug for ill-conditioned matrices (PR#594)

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._