Should one use both Eigen's Tensor.contract Tensor.reduce to achieve similar functionality to np.einsum?

Say I want to multiply a=shape(3, 2, 2) with b=shape(3, 2, 2) matrices together. With np.einsum this can easily be expressed as np.einsum("sij,sjk->sil", a, b).

I'm trying to do something in similar in Eigen using Tensor.contract but at most I've been able to produce a (3,2,3,2) output shape (see output c below).

#include <iostream>

#include "third_party/eigen3/unsupported/Eigen/CXX11/Tensor"


int main(void)
{
    // np.einsum("sij,sjk->sik")
    Eigen::DefaultDevice my_device;
    Eigen::Tensor<float, 3> a(3, 2, 2);
    Eigen::Tensor<float, 3> b(3, 2, 2);

    Eigen::array<Eigen::IndexPair<int>,1> dims = {Eigen::IndexPair<int>(2, 1)};

    a.setConstant(1);
    b.setConstant(2);

    Eigen::Tensor<float, 4> c(3, 2, 3, 2);
    Eigen::Tensor<float, 3> d(3, 2, 2);
    Eigen::array<int, 1> reduce_dims({2});

    c.device(my_device) = a.contract(b, dims);
    d.device(my_device) = a.contract(b, dims).sum(reduce_dims);
    std::cout << c << std::endl;
    std::cout << d << std::endl;
}

Questions

1. It seems that contract and reduce must be combined to produce the einsum behaviour. Is this correct?

2. This is a simple case, but I want to end up performing more complex einstein-summation-like behaviour. Something to the effect of np.einsum("stacij,tacjk, tckl->tacil"). Is it reasonable to to use Eigen's contract and reduce to achieve this?

3. The ijkl dimensions will be at most 4, while the other dimensions will be much larger. Will Eigen perform well with these dimension sizes.