title: Why is a.dot(b) faster than a@b although Numpy recommends a@b
body: <p>According to the answers from this <a href="https://stackoverflow.com/questions/3890621/how-does-multiplication-differ-for-numpy-matrix-vs-array-classes#">question</a> and also according to <a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html?highlight=matrix%20multiplication" rel="noreferrer">numpy</a>, matrix multiplication of 2-D arrays is best done via <code>a @ b</code>, or <code>numpy.matmul(a,b)</code> as compared to <code>a.dot(b)</code>.</p>
  
  <blockquote>
    <p>If both a and b are 2-D arrays, it is matrix multiplication, but using
    matmul or a @ b is preferred.</p>
  </blockquote>
  
  <p>I did the following benchmark and found contrary results. </p>
  
  <p><strong>Questions:</strong> Is there an issue with my benchmark? If not, why does Numpy not recommend <code>a.dot(b)</code> when it is faster than <code>a@b</code> or <code>numpy.matmul(a,b)</code>?</p>
  
  <p>Benchmark used python 3.5 numpy 1.15.0.</p>
  
  <pre><code>$ pip3 list | grep numpy
  numpy                         1.15.0
  $ python3 --version
  Python 3.5.2
  </code></pre>
  
  <p><strong>Benchmark code:</strong></p>
  
  <pre><code>import timeit
  
  setup = '''
  import numpy as np
  a = np.arange(16).reshape(4,4)
  b = np.arange(16).reshape(4,4)
  ''' 
  test = '''
  for i in range(1000):
      a @ b
  '''
  test1 = '''
  for i in range(1000):
      np.matmul(a,b)
  '''
  test2 = '''
  for i in range(1000):
      a.dot(b)
  '''
  
  print( timeit.timeit(test, setup, number=100) )
  print( timeit.timeit(test1, setup, number=100) )
  print( timeit.timeit(test2, setup, number=100) )
  </code></pre>
  
  <p><strong>Results:</strong></p>
  
  <pre><code>test  : 0.11132473500038031
  test1 : 0.10812476599676302
  test2 : 0.06115105600474635
  </code></pre>
  
  <p><strong>Add on results:</strong></p>
  
  <pre><code>&gt;&gt;&gt; a = np.arange(16).reshape(4,4)
  &gt;&gt;&gt; b = np.arange(16).reshape(4,4)
  &gt;&gt;&gt; a@b
  array([[ 56,  62,  68,  74],
         [152, 174, 196, 218],
         [248, 286, 324, 362],
         [344, 398, 452, 506]])
  &gt;&gt;&gt; np.matmul(a,b)
  array([[ 56,  62,  68,  74],
         [152, 174, 196, 218],
         [248, 286, 324, 362],
         [344, 398, 452, 506]])
  &gt;&gt;&gt; a.dot(b)
  array([[ 56,  62,  68,  74],
         [152, 174, 196, 218],
         [248, 286, 324, 362],
         [344, 398, 452, 506]])
  </code></pre>