Broadcast operator over shared index

Question

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.

Answer 1

You could write the loops explicitly as you show.

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

So for your application, it would be something like:

auto del = delta(idx_samp,prime(idx_samp),prime(idx_samp,2));
auto V = (P*U*prime(Q,idx_samp)*del).noPrime();

As is mentioned in the other post, this creates some unnecessary intermediate tensors so may not be the most efficient solution. Please let us know if either approach is not fast enough, it would be nice to have a more direct and efficient way of doing this in ITensor (but that kind of operation has not been very common so far in physics applications, so hasn't been a high priority to make special functionality for it).