[Rd] Numerical accuracy of matrix multiplication

Alexis Sarda alexis.sarda at gmail.com
Tue Sep 20 17:33:49 CEST 2016


I just realized that I was actually using a different random number
generator, could that be a valid reason for the discrepancy?

The code should be:

RNGkind("L'Ecuyer")
set.seed(883)
x <- rnorm(100)
x %*% x - sum(x^2) # equal to 1.421085e-14


Regards,
Alexis Sarda.



On Tue, Sep 20, 2016 at 5:27 PM, Martin Maechler <maechler at stat.math.ethz.ch
> wrote:

> >>>>> peter dalgaard <pdalgd at gmail.com>
> >>>>>     on Fri, 16 Sep 2016 13:33:11 +0200 writes:
>
>     > On 16 Sep 2016, at 12:41 , Alexis Sarda <alexis.sarda at gmail.com>
> wrote:
>
>     >> Hello,
>     >>
>     >> while testing the crossprod() function under Linux, I noticed the
> following:
>     >>
>     >> set.seed(883)
>     >> x <- rnorm(100)
>     >> x %*% x - sum(x^2) # equal to 1.421085e-14
>     >>
>     >> Is this difference normal? It seems to be rather large for double
> precision.
>     >>
>
>     > It's less than .Machine$double.eps, relative (!) to x  %*% x ~= 100.
>
> indeed!
>
> Still, it gives exactly 0 on my platform(s), where I'm using R's
> own version of BLAS / Lapack.
>
> Are you perhaps using an "optimized" BLAS / LAPACK , i.e, one
> that is fast but slightly less so accurate ?
>
> Martin Maechler,
> ETH Zurich
>
>
>     > -pd
>
>     >> Regards,
>     >> Alexis.
>     >>
>     >> [[alternative HTML version deleted]]
>     >>
>     >> ______________________________________________
>     >> R-devel at r-project.org mailing list
>     >> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>     > --
>     > Peter Dalgaard, Professor,
>     > Center for Statistics, Copenhagen Business School
>     > Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>     > Phone: (+45)38153501
>     > Office: A 4.23
>     > Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
>
>     > ______________________________________________
>     > R-devel at r-project.org mailing list
>     > https://stat.ethz.ch/mailman/listinfo/r-devel
>

	[[alternative HTML version deleted]]



More information about the R-devel mailing list