I am not using numpy but Eigen::Tensor C++ API, which only has contraction operations, this is just to enable me think through implementation from python.
So 'ij, ijk -> ik' is basically like doing a for loop for each of the first dimensions.
a = np.random.uniform(size=[10, 4])
b = np.random.uniform(size=[10, 4, 4])
vec = []
for i in range(10):
vec.append(a[i].dot(b[i]))
print(np.stack(vec, axis=0))
## or with einsum
print(np.einsum('ij,ijk->ik', a, b))
This can not seem to be done easily with tensordot. Any suggestions?