title: Why is np.dot imprecise? (n-dim arrays)
body: <p>Suppose we take <code>np.dot</code> of two <code>'float32'</code> 2D arrays:  </p>
  
  <pre class="lang-py prettyprint-override"><code>res = np.dot(a, b)   # see CASE 1
  print(list(res[0]))  # list shows more digits
  </code></pre>
  
  <pre class="lang-py prettyprint-override"><code>[-0.90448684, -1.1708503, 0.907136, 3.5594249, 1.1374011, -1.3826287]
  </code></pre>
  
  <p>Numbers. Except, they can change:</p>
  
  <hr>
  
  <p><strong>CASE 1</strong>: slice <code>a</code></p>
  
  <pre class="lang-py prettyprint-override"><code>np.random.seed(1)
  a = np.random.randn(9, 6).astype('float32')
  b = np.random.randn(6, 6).astype('float32')
  
  for i in range(1, len(a)):
      print(list(np.dot(a[:i], b)[0])) # full shape: (i, 6)
  </code></pre>
  
  <pre><code>[-0.9044868,  -1.1708502, 0.90713596, 3.5594249, 1.1374012, -1.3826287]
  [-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
  [-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
  [-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
  [-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
  [-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
  [-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
  [-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
  </code></pre>
  
  <p>Results differ, even though the printed slice derives from the exact same numbers multiplied.
  <hr>
  <strong>CASE 2</strong>: flatten <code>a</code>, take a 1D version of <code>b</code>, <em>then</em> slice <code>a</code>:</p>
  
  <pre class="lang-py prettyprint-override"><code>np.random.seed(1)
  a = np.random.randn(9, 6).astype('float32')
  b = np.random.randn(1, 6).astype('float32')
  
  for i in range(1, len(a)):
      a_flat = np.expand_dims(a[:i].flatten(), -1) # keep 2D
      print(list(np.dot(a_flat, b)[0])) # full shape: (i*6, 6)
  </code></pre>
  
  <pre class="lang-py prettyprint-override"><code>[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  [-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  [-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  [-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  [-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  [-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  [-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  [-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
  </code></pre>
  
  <hr>
  
  <p><strong>CASE 3</strong>: stronger control; set all non-involved entires to <em>zero</em>: add <code>a[1:] = 0</code> to CASE 1 code. Result: discrepancies persist.</p>
  
  <hr>
  
  <p><strong>CASE 4</strong>: check indices other than <code>[0]</code>; like for <code>[0]</code>, results begin to stabilize a fixed # of array enlargements from their point of creation. <a href="https://pastebin.com/wwCYiKxp" rel="nofollow noreferrer">Output</a></p>
  
  <pre class="lang-py prettyprint-override"><code>np.random.seed(1)
  a = np.random.randn(9, 6).astype('float32')
  b = np.random.randn(6, 6).astype('float32')
  
  for j in range(len(a) - 2):
      for i in range(1, len(a)):
          res = np.dot(a[:i], b)
          try:    print(list(res[j]))
          except: pass
      print()
  </code></pre>
  
  <hr>
  
  <p>Hence, for the 2D * 2D case, results differ - but are consistent for 1D * 1D. From some of my readings, this appears to stem from 1D-1D using simple addition, whereas 2D-2D uses 'fancier', performance-boosting addition that may be less precise (e.g. pairwise addition does the opposite). Nonetheless, I'm unable to understand why discrepancies vanish in case 1 once <code>a</code> is sliced past a set 'threshold'; the larger <code>a</code> and <code>b</code>, the later this threshold seems to lie, but it always exists.</p>
  
  <p>All said: why is <code>np.dot</code> imprecise (and inconsistent) for ND-ND arrays? <a href="https://github.com/numpy/numpy/issues/14917" rel="nofollow noreferrer">Relevant Git</a></p>
  
  <hr>
  
  <p><strong>Additional info</strong>:</p>
  
  <ul>
  <li><b>Environment</b>: Win-10 OS, Python 3.7.4, Spyder 3.3.6 IDE, Anaconda 3.0 2019/10</li>
  <li><b>CPU</b>: <a href="https://ark.intel.com/content/www/us/en/ark/products/97185/intel-core-i7-7700hq-processor-6m-cache-up-to-3-80-ghz.html" rel="nofollow noreferrer">i7-7700HQ 2.8 GHz</a></li>
  <li><b>Numpy</b> v1.16.5</li>
  </ul>
  
  <p><strong>Possible culprit library</strong>: Numpy <a href="https://docs.anaconda.com/mkl-optimizations/" rel="nofollow noreferrer">MKL</a> - also BLASS libraries; thanks to <a href="https://stackoverflow.com/users/1004096/bi-rico">Bi Rico</a> for noting</p>
  
  <hr>
  
  <p><strong>Stress-test code</strong>: as noted, discrepancies exacerbate in frequency w/ larger arrays; if above isn't reproducible, below should be (if not, try larger dims). <a href="https://pastebin.com/kGFY3vwZ" rel="nofollow noreferrer">My output</a></p>
  
  <pre class="lang-py prettyprint-override"><code>np.random.seed(1)
  a = (0.01*np.random.randn(9, 9999)).astype('float32') # first multiply then type-cast
  b = (0.01*np.random.randn(9999, 6)).astype('float32') # *0.01 to bound mults to &lt; 1
  
  for i in range(1, len(a)):
      print(list(np.dot(a[:i], b)[0]))
  </code></pre>
  
  <hr>
  
  <p><strong>Problem severity</strong>: shown discrepancies are 'small', but no longer so when operating on a neural network with billions of numbers multiplied over a few seconds, and trillions over the entire runtime; reported model accuracy differs by entire 10's of percents, per <a href="https://github.com/keras-team/keras/issues/12400" rel="nofollow noreferrer">this thread</a>. </p>
  
  <p>Below is a gif of arrays resulting from feeding to a model what's basically <code>a[0]</code>, w/ <code>len(a)==1</code> vs. <code>len(a)==32</code>:</p>
  
  <p><img src="https://i.stack.imgur.com/nOqIG.gif" width="350"></p>
  
  <hr>
  
  <p><strong>OTHER PLATFORMS</strong> results, according and with thanks to <a href="https://stackoverflow.com/users/103081/paul">Paul</a>'s testing:</p>
  
  <p><strong>Case 1 reproduced (partly)</strong>:</p>
  
  <ul>
  <li>Google Colab VM -- Intel Xeon 2.3 G-Hz -- Jupyter -- Python 3.6.8</li>
  <li>Win-10 Pro Docker Desktop -- Intel i7-8700K -- jupyter/scipy-notebook -- Python 3.7.3</li>
  <li>Ubuntu 18.04.2 LTS + Docker -- AMD FX-8150 --  jupyter/scipy-notebook -- Python 3.7.3</li>
  </ul>
  
  <p><em>Note</em>: these yield much lower error than shown above; two entries on the first row are off by 1 in the least significant digit from corresponding entries in the other rows.</p>
  
  <p><strong>Case 1 not reproduced</strong>:</p>
  
  <ul>
  <li>Ubuntu 18.04.3 LTS -- Intel i7-8700K -- IPython 5.5.0 -- Python 2.7.15+ and 3.6.8 (2 tests)  </li>
  <li>Ubuntu 18.04.3 LTS -- Intel i5-3320M -- IPython 5.5.0 -- Python 2.7.15+</li>
  <li>Ubuntu 18.04.2 LTS -- AMD FX-8150 --  IPython 5.5.0 -- Python 2.7.15rc1 </li>
  </ul>
  
  <p><em>Notes</em>:</p>
  
  <ul>
  <li>The <a href="https://colab.research.google.com/drive/1yekKp3lsIzHPKMFcQ_Rs60E3ahaPQEd_" rel="nofollow noreferrer">linked</a> Colab notebook and jupyter environments show a far lesser discrepancy (and only for first two rows) than is observed on my system. Also, Case 2 never (yet) showed imprecision.</li>
  <li>Within this very limited sample, the current (Dockerized) Jupyter environment is more susceptible than the IPython environment. </li>
  <li><code>np.show_config()</code> too long to post, but in summary: IPython envs are BLAS/LAPACK-based; Colab is OpenBLAS-based. In IPython Linux envs, BLAS libraries are system-installed -- in Jupyter and Colab, they come from /opt/conda/lib</li>
  </ul>
  
  <hr>
  
  <p><strong>UPDATE</strong>: the accepted answer is accurate, but broad and incomplete. The question remains open for anyone who can explain the behavior at the <em>code level</em> - namely, an exact algorithm used by <code>np.dot</code>, and how it explains 'consistent inconsistencies' observed in above results (also see comments). Here are some direct implementations beyond my deciphering: <a href="https://www.netlib.org/clapack/cblas/sdot.c" rel="nofollow noreferrer">sdot.c</a> -- <a href="https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/arraytypes.c.src#L3534" rel="nofollow noreferrer">arraytypes.c.src</a></p>