I've been messing around a bit with ITensor for some ML/signal processing applications and I think the library is awesome! However, I'm not sure how to go about the following: I have two "matrix" tensors - say P and Q - whose rows represent samples from a signal and whose columns represent features of the signal. In code:

uint num_samp = 100;
auto idx_samp = Index(num_samp "samples");
auto idx_f1 = Index(2, "features1");
auto idx_f2 = Index(2, "features2");
auto P = randomITensor(idx_samp, idx_f1);
auto Q = randomITensor(idx_samp, idx_f2);


Now I have some "operator" U indexed by idx_f1 and idx_f2 and some other index idx_u. I want to contract each row of P and Q with U producing a new tensor indexed by idx_samp and idx_u. In pseudocode:

U = randomITensor(idx_f1, idx_f2, idx_u);
V = ITensor(idx_samp, idx_u);
for n in 1:num_samp
V[idx_samp=n,:] = P[idx_samp=n,:]*U*Q[idx_samp=n,:]


But I'm not sure how to go about doing this in ITensor. Of course, I could just treat everything as a matrix and not worry about it, but there are some other steps where it's nice to have ITensor keep track of bookeeping on indices. I'm new both to ITensor and thinking about things as tensors so apologies if the question has a straightforward answer.

An alternative is to write it in terms of delta tensors to constrain idx_samp of P and Q to be the same value, similar to this question here: http://www.itensor.org/support/1407/computing-the-partial-kronecker-product-of-two-itensors?show=1407#q1407
auto del = delta(idx_samp,prime(idx_samp),prime(idx_samp,2));