title: Matrix multiplication on CPU (numpy) and GPU (gnumpy) give different results
body: <p>I'm using <a href="http://www.cs.toronto.edu/~tijmen/gnumpy.html" rel="noreferrer">gnumpy</a> to speed up some computations in training a neural network by doing them on GPU.
  I'm getting the desired speedup but am a little bit worried about the differences in the results of numpy (cpu) vs gnumpy (gpu).</p>
  
  <p>I have the following test script to illustrate the problem:</p>
  
  <pre><code>import gnumpy as gpu
  import numpy as np
  
  n = 400
  
  a = np.random.uniform(low=0., high=1., size=(n, n)).astype(np.float32)
  b = np.random.uniform(low=0., high=1., size=(n, n)).astype(np.float32)
  
  ga = gpu.garray(a)
  gb = gpu.garray(b)
  
  ga = ga.dot(gb)
  a  = a.dot(b)
  
  print ga.as_numpy_array(dtype=np.float32) - a
  </code></pre>
  
  <p>which provides the output:</p>
  
  <pre><code>[[  1.52587891e-05  -2.28881836e-05   2.28881836e-05 ...,  -1.52587891e-05
      3.81469727e-05   1.52587891e-05]
   [ -5.34057617e-05  -1.52587891e-05   0.00000000e+00 ...,   1.52587891e-05
      0.00000000e+00   1.52587891e-05]
   [ -1.52587891e-05  -2.28881836e-05   5.34057617e-05 ...,   2.28881836e-05
      0.00000000e+00  -7.62939453e-06]
   ..., 
   [  0.00000000e+00   1.52587891e-05   3.81469727e-05 ...,   3.05175781e-05
      0.00000000e+00  -2.28881836e-05]
   [  7.62939453e-06  -7.62939453e-06  -2.28881836e-05 ...,   1.52587891e-05
      7.62939453e-06   1.52587891e-05]
   [  1.52587891e-05   7.62939453e-06   2.28881836e-05 ...,  -1.52587891e-05
      7.62939453e-06   3.05175781e-05]]
  </code></pre>
  
  <p>As you can see, the differences are around the magnitude of 10^-5.</p>
  
  <p>So the question is: <strong>should I be worried about these differences or is this the expected behaviour?</strong></p>
  
  <p>Additional information:</p>
  
  <ul>
  <li>GPU: GeForce GTX 770;</li>
  <li>numpy version: 1.6.1</li>
  </ul>
  
  <p>I noticed the problem when I used gradient checking (with finite difference approximation) to verify that the small modifications I made to switch from numpy to gnumpy didn't break anything. As one may expect the gradient checking did not work with 32 bit precision (gnumpy does not support float64), but to my surprise the errors differed between CPU and GPU when using the same precision.</p>
  
  <p>The errors on CPU and GPU on a small test neural network are given below:
  <img src="https://i.imgur.com/IBj3eKt.png?1" alt="gradient checking errors"></p>
  
  <p>Since the error magnitudes are similar, I guess that these differences are OK?</p>
  
  <p>After reading the <a href="http://developer.download.nvidia.com/assets/cuda/files/NVIDIA-CUDA-Floating-Point.pdf" rel="noreferrer">article</a>, referenced in the comment by BenC, I'm quite sure that the differences can be mostly explained by one of the devices using the fused multiply-add (FMA) instruction and the other not.</p>
  
  <p>I implemented the example from the paper:</p>
  
  <pre><code>import gnumpy as gpu
  import numpy as np
  
  a=np.array([1.907607,-.7862027, 1.147311, .9604002], dtype=np.float32)
  b=np.array([-.9355000, -.6915108, 1.724470, -.7097529], dtype=np.float32)
  
  ga = gpu.garray(a)
  gb = gpu.garray(b)
  
  ga = ga.dot(gb)
  a  = a.dot(b)
  
  print "CPU", a
  print "GPU", ga
  print "DIFF", ga - a
  
  &gt;&gt;&gt;CPU 0.0559577
  &gt;&gt;&gt;GPU 0.0559577569366
  &gt;&gt;&gt;DIFF 8.19563865662e-08
  </code></pre>
  
  <p>...and the difference is similar to FMA vs serial algorithm (though for some reason both results differ from the exact result more than in the paper).</p>
  
  <p>The GPU I'm using (GeForce GTX 770) supports FMA instruction while the CPU does not (I have an Ivy Bridge Intel® Xeon(R) CPU E3-1225 V2, but intel introduced the FMA3 instruction in their products with Haswell).</p>
  
  <p>Other possible explanations include the different math libraries used in the background or differences in the sequence of operations caused by, for example, the different level of parallelization on CPU vs GPU.</p>