Origin: https://github.com/dask/dask/issues/4361
Summary: Here's a temporary work-around that resolve the issue with NumPy 1.16.
 See also: https://github.com/numpy/numpy/issues/12712 for discussion
 about the bug.

diff --git a/dask/array/core.py b/dask/array/core.py
index 5cf0f082a..bdb31a833 100644
--- a/dask/array/core.py
+++ b/dask/array/core.py
@@ -952,6 +952,10 @@ def __array_ufunc__(self, numpy_ufunc, method, *inputs, **kwargs):
                 return NotImplemented
 
         if method == '__call__':
+            if numpy_ufunc is np.matmul:
+                from .routines import matmul
+                # special case until apply_gufunc handles optional dimensions
+                return matmul(*inputs, **kwargs)
             if numpy_ufunc.signature is not None:
                 from .gufunc import apply_gufunc
                 return apply_gufunc(numpy_ufunc,
diff --git a/dask/array/tests/test_array_core.py b/dask/array/tests/test_array_core.py
index 97e6b0b72..2f73688e7 100644
--- a/dask/array/tests/test_array_core.py
+++ b/dask/array/tests/test_array_core.py
@@ -857,6 +857,16 @@ def test_matmul():
     assert_eq(operator.matmul(z, a), operator.matmul(c, x))
 
 
+def test_matmul_array_ufunc():
+    # regression test for https://github.com/dask/dask/issues/4353
+    x = np.random.random((5, 5))
+    y = np.random.random((5, 2))
+    a = from_array(x, chunks=(1, 5))
+    b = from_array(y, chunks=(5, 1))
+    result = b.__array_ufunc__(np.matmul, '__call__', a, b)
+    assert_eq(result, x.dot(y))
+
+
 def test_T():
     x = np.arange(400).reshape((20, 20))
     a = from_array(x, chunks=(5, 5))