On two laptops running 32-bit kubuntu, I have found that svd(), invoked
within R 2.9.1 as supplied with the current ubuntu package, returns very
incorrect results when presented with complex-valued input. One of the
laptops is a Dell D620, the other a MacBook Pro. I've also verified the
problem on a 32-bit desktop. On these same systems, R compiled from
source provides apparently correct results, so long as the --with-blas
option is not used or --with-blas is used and an atlas library compiled
from source is specified as the blas. The blas that gives incorrect
results is the ubuntu sse2 atlas package. On a 64-bit ubuntu desktop
system I find the svd() results are correct with the ubuntu packaged R
as well as with R compiled on the machine.
The error seems not to arise, at least sometimes, with simple 2x2
complex matrices. To check for the error, run the following:
a <- matrix(rnorm(16) + rnorm(16)*1i, 4)
b <- svd(a)
with(b, u %*% t(Conj(u)))
The result should be an identity matrix, by the definition of the
singular value decomposition.
Below are example runs, first on the 64-bit system where results are
correct, then on one of the 32-bit systems.
> version
_
platform x86_64-pc-linux-gnu
arch x86_64
os linux-gnu
system x86_64, linux-gnu
status
major 2
minor 9.1
year 2009
month 06
day 26
svn rev 48839
language R
version.string R version 2.9.1 (2009-06-26)
> a <- matrix(rnorm(16) + 1i * rnorm(16),4)
> a
[,1] [,2] [,3]
[1,] 1.064712-1.343633i -1.6787892-0.7356784i -1.3110026-1.7295312i
[2,] -0.829329-1.538217i 0.6514795+0.8242854i 1.3948013-0.4340075i
[3,] 2.241843+0.297037i 2.9147970+0.2397768i -0.5197081-0.5796579i
[4,] 1.567566-1.154438i -1.0900313+0.0121055i 0.0141203+0.6139178i
[,4]
[1,] 1.344998+1.231298i
[2,] -0.190861+1.817582i
[3,] 1.855617-1.244282i
[4,] -1.735007+0.725572i
> b <- svd(a)
> with(b, u %*% t(Conj(u)))
[,1] [,2]
[1,] 1.000000e+00+0.000000e+00i 9.971272e-18-1.092774e-16i
[2,] 9.971272e-18+1.062857e-16i 1.000000e+00-0.000000e+00i
[3,] -9.500311e-17+3.729126e-17i -1.319715e-16-3.354759e-16i
[4,] 1.585713e-16+2.581553e-16i 2.988332e-17+1.814683e-16i
[,3] [,4]
[1,] -9.500311e-17-2.546774e-17i 1.585713e-16-2.644267e-16i
[2,] -1.319715e-16+3.350337e-16i 2.988332e-17-1.820918e-16i
[3,] 1.000000e+00-0.000000e+00i -2.445638e-16+1.528022e-16i
[4,] -2.445638e-16-1.735266e-16i 1.000000e+00+0.000000e+00i
>
-----------------------------------------------------
> version
_
platform i486-pc-linux-gnu
arch i486
os linux-gnu
system i486, linux-gnu
status
major 2
minor 9.1
year 2009
month 06
day 26
svn rev 48839
language R
version.string R version 2.9.1 (2009-06-26)
> a <- matrix(rnorm(16) + 1i * rnorm(16), 4)
> a
[,1] [,2] [,3]
[1,] 1.4465327-1.4747010i 0.7291287-0.3789878i 1.1978147-0.2393877i
[2,] 0.8557766-0.9561998i -1.8390624-0.2263745i 0.1601790+0.8254724i
[3,] 1.1659289+1.3202786i 1.1182754-0.0327039i 0.9901992+0.5903004i
[4,] -0.4728973-0.2885660i -0.2541130+1.5875047i -0.9362061+0.1858397i
[,4]
[1,] -0.5949982+1.2969862i
[2,] 0.6272468+1.2662726i
[3,] 0.6722278-0.7768301i
[4,] 0.4064206+0.6788890i
> b <- svd(a)
> with(b, u %*% t(Conj(u)))
[,1] [,2] [,3]
[1,] 0.7429858-0.0000000i -0.0890804-0.1273023i -0.1145715-0.3955123i
[2,] -0.0890804+0.1273023i 0.8251623+0.0000000i -0.0217231-0.2699445i
[3,] -0.1145715+0.3955123i -0.0217231+0.2699445i 1.3072318-0.0000000i
[4,] -0.1236460+0.0628505i 0.1937076-0.0101354i -0.1682549+0.0544654i
[,4]
[1,] -0.1236460-0.0628505i
[2,] 0.1937076+0.0101354i
[3,] -0.1682549-0.0544654i
[4,] 1.0912424-0.0000000i
>
Dirk Eddelbuettel
2009-Aug-09 16:25 UTC
[R-sig-Debian] Inaccuracy in svd() with R ubuntu package
Hi Chris, Great bug report! On 9 August 2009 at 11:51, Chris Sims wrote: | On two laptops running 32-bit kubuntu, I have found that svd(), invoked | within R 2.9.1 as supplied with the current ubuntu package, returns very | incorrect results when presented with complex-valued input. One of the | laptops is a Dell D620, the other a MacBook Pro. I've also verified the | problem on a 32-bit desktop. On these same systems, R compiled from | source provides apparently correct results, so long as the --with-blas | option is not used or --with-blas is used and an atlas library compiled | from source is specified as the blas. The blas that gives incorrect | results is the ubuntu sse2 atlas package. On a 64-bit ubuntu desktop What happens when you remove the libatlas-sse2 package and keep the libatlas-base package? On my 32bit Debian system, I only see incorrect results (as per your test below) when libatlas3gf-sse2 is installed. The Atlas package is undergoing a transition. The previous maintainer stepped down, and Sylvestre Ledu (CC'ed) is almost done with new packages based on Atlas 3.8.3 -- they are in the Alioth SVN if you want to take a peek. Now, Sylvestre is currently travelling for a few weeks but this should get better 'real soon'. Hth, Dirk | system I find the svd() results are correct with the ubuntu packaged R | as well as with R compiled on the machine. | | The error seems not to arise, at least sometimes, with simple 2x2 | complex matrices. To check for the error, run the following: | | a <- matrix(rnorm(16) + rnorm(16)*1i, 4) | b <- svd(a) | with(b, u %*% t(Conj(u))) | | The result should be an identity matrix, by the definition of the | singular value decomposition. | | Below are example runs, first on the 64-bit system where results are | correct, then on one of the 32-bit systems. | | > version | _ | platform x86_64-pc-linux-gnu | arch x86_64 | os linux-gnu | system x86_64, linux-gnu | status | major 2 | minor 9.1 | year 2009 | month 06 | day 26 | svn rev 48839 | language R | version.string R version 2.9.1 (2009-06-26) | > a <- matrix(rnorm(16) + 1i * rnorm(16),4) | > a | [,1] [,2] [,3] | [1,] 1.064712-1.343633i -1.6787892-0.7356784i -1.3110026-1.7295312i | [2,] -0.829329-1.538217i 0.6514795+0.8242854i 1.3948013-0.4340075i | [3,] 2.241843+0.297037i 2.9147970+0.2397768i -0.5197081-0.5796579i | [4,] 1.567566-1.154438i -1.0900313+0.0121055i 0.0141203+0.6139178i | [,4] | [1,] 1.344998+1.231298i | [2,] -0.190861+1.817582i | [3,] 1.855617-1.244282i | [4,] -1.735007+0.725572i | > b <- svd(a) | > with(b, u %*% t(Conj(u))) | [,1] [,2] | [1,] 1.000000e+00+0.000000e+00i 9.971272e-18-1.092774e-16i | [2,] 9.971272e-18+1.062857e-16i 1.000000e+00-0.000000e+00i | [3,] -9.500311e-17+3.729126e-17i -1.319715e-16-3.354759e-16i | [4,] 1.585713e-16+2.581553e-16i 2.988332e-17+1.814683e-16i | [,3] [,4] | [1,] -9.500311e-17-2.546774e-17i 1.585713e-16-2.644267e-16i | [2,] -1.319715e-16+3.350337e-16i 2.988332e-17-1.820918e-16i | [3,] 1.000000e+00-0.000000e+00i -2.445638e-16+1.528022e-16i | [4,] -2.445638e-16-1.735266e-16i 1.000000e+00+0.000000e+00i | > | | ----------------------------------------------------- | > version | _ | platform i486-pc-linux-gnu | arch i486 | os linux-gnu | system i486, linux-gnu | status | major 2 | minor 9.1 | year 2009 | month 06 | day 26 | svn rev 48839 | language R | version.string R version 2.9.1 (2009-06-26) | > a <- matrix(rnorm(16) + 1i * rnorm(16), 4) | > a | [,1] [,2] [,3] | [1,] 1.4465327-1.4747010i 0.7291287-0.3789878i 1.1978147-0.2393877i | [2,] 0.8557766-0.9561998i -1.8390624-0.2263745i 0.1601790+0.8254724i | [3,] 1.1659289+1.3202786i 1.1182754-0.0327039i 0.9901992+0.5903004i | [4,] -0.4728973-0.2885660i -0.2541130+1.5875047i -0.9362061+0.1858397i | [,4] | [1,] -0.5949982+1.2969862i | [2,] 0.6272468+1.2662726i | [3,] 0.6722278-0.7768301i | [4,] 0.4064206+0.6788890i | > b <- svd(a) | > with(b, u %*% t(Conj(u))) | [,1] [,2] [,3] | [1,] 0.7429858-0.0000000i -0.0890804-0.1273023i -0.1145715-0.3955123i | [2,] -0.0890804+0.1273023i 0.8251623+0.0000000i -0.0217231-0.2699445i | [3,] -0.1145715+0.3955123i -0.0217231+0.2699445i 1.3072318-0.0000000i | [4,] -0.1236460+0.0628505i 0.1937076-0.0101354i -0.1682549+0.0544654i | [,4] | [1,] -0.1236460-0.0628505i | [2,] 0.1937076+0.0101354i | [3,] -0.1682549-0.0544654i | [4,] 1.0912424-0.0000000i | > | | | ---------------------------------------------------------------------- | _______________________________________________ | R-SIG-Debian mailing list | R-SIG-Debian at r-project.org | https://stat.ethz.ch/mailman/listinfo/r-sig-debian -- Three out of two people have difficulties with fractions.