Is this a bug?
In [21]: a = ones((2,2))
In [22]: b = ones((2,2,2,2))
In [23]: tensordot(a, tensordot(b, a, 2), 2)
Out[23]: array(16.0)
It seems to me that consistency with the dot product would require a scalar result, not a 0-dim array.
Also, do we have a plain old tensor product? Outer flattens the indices.
Chuck