[R] Decimal Accuracy Loss?

[peter dalgaard]

On Apr 6, 2011, at 17:58 , Petr Savicky wrote:

> On Wed, Apr 06, 2011 at 11:33:48AM -0400, Brigid Mooney wrote:
>> This is hopefully a quick question on decimal accuracy.  Is any
>> decimal accuracy lost when casting a numeric vector as a matrix?  And
>> then again casting the result back to a numeric?
>> 
>> I'm finding that my calculation values are different when I run for
>> loops that manually calculate matrix multiplication as compared to
>> when I cast the vectors as matrices and multiply them using "%*%".
>> (The errors are very small, but the process is run iteratively
>> thousands of times, at which point the error between the two
>> differences becomes noticeable.)
>> 
>> I've read FAQ # 7.31 "Why doesn't R think these numbers are equal?",
>> but just want to confirm that the differences in values are due to
>> differences in the matrix multiplication operator and manual
>> calculation via for loops, rather than information that is lost when
>> casting a numeric as a matrix and back again.
> 
> Others already confirmed that casting a numeric as a matrix and back
> again does not change the numbers. It is likely that the library
> operator "%*%" is more accurate than a straightforward for loop.
> For example, sum(x) uses a more accurate algorithm than iteration
> of s <- s + x[i] in double precision.

Even more likely, %*% is optimized for speed by reordering the additions and multiplications to take advantage of pipelining and CPU caches at several levels. This may or may not improve accuracy, but certainly does affect the last bits in the results. 

-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com