[SciPy-User] understanding machine precision

[josef.pktd at gmail.com]

On Tue, Dec 14, 2010 at 1:09 PM, Francesc Alted <faltet at pytables.org> wrote:
> A Tuesday 14 December 2010 18:42:40 josef.pktd at gmail.com escrigué:
>> I thought that we get deterministic results, with identical machine
>> precision errors, but I get (with some random a0, b0)
>>
>> >>> for i in range(5):
>>       x = scipy.linalg.lstsq(a0,b0)[0]
>>       x2 = scipy.linalg.lstsq(a0,b0)[0]
>>       print np.max(np.abs(x-x2))
>>
>>
>> 9.99200722163e-016
>> 9.99200722163e-016
>> 0.0
>> 0.0
>> 9.99200722163e-016
>>
>> >>> a0.shape
>>
>> (100, 10)
>>
>> >>> b0.shape
>>
>> (100, 3)
>>
>> Why is the result not always the same? just curious
>
> That's really funny!  Could you please come up with a self-contained
> example so as to see if others can reproduce that?

Essentially all I did was

    a0 = np.random.randn(100,10)
    b0 = a0.sum(1)[:,None] + np.random.randn(100,3)

I copied scipy.linalg pinv, pinv2 and lstsq to a local module, and the
results where not exactly the same as with the ones in scipy.
So, I did this check for scipy also.

the attached script produces different results in each run (on
WindowsXP 32), for example

lstsq
1.55431223448e-015
1.55431223448e-015
0.0
1.55431223448e-015
1.55431223448e-015

pinv
5.20417042793e-017
5.20417042793e-017
0.0
0.0
0.0

pinv2
0.0
0.0
5.76795555762e-017
5.76795555762e-017
4.51028103754e-017

Thanks,

Josef

>
> --
> Francesc Alted
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------

import numpy as np
import scipy.linalg

np.random.seed(12345)
a0 = np.random.randn(100,10)
b0 = a0.sum(1)[:,None] + np.random.randn(100,3)


print '\nlstsq'
for i in range(5):
    x = scipy.linalg.lstsq(a0,b0)[0]
    x2 = scipy.linalg.lstsq(a0,b0)[0]
    print np.max(np.abs(x-x2))

print '\npinv'
for i in range(5):
    x = scipy.linalg.pinv(a0)
    x2 = scipy.linalg.pinv(a0)
    print np.max(np.abs(x-x2))

print '\npinv2'
for i in range(5):
    x = scipy.linalg.pinv2(a0)
    x2 = scipy.linalg.pinv2(a0)
    print np.max(np.abs(x-x2))