Computing the partial Kronecker product of two ITensors

Question

I'm implementing several "fundamental" operations defined in arXiv:1405.7786 for matrix product states, such as the Hadamard product. Part of doing this is defining operations on the core tensors, such as the partial Kronecker product (see Definition 2.5 of the aforementioned paper, reproduced below). Define two tensors @@\mathbf{A}@@ and @@\mathbf{B}@@, where

$$
A \in \mathbb{R}^{R_1 \times R_2 \times \cdots \times R_M \times I_1 \times \cdots I_N}
$$

and

$$
B \in \mathbb{R}^{S_1 \times S_2 \times \cdots \times S_M \times I_1 \times \cdots I_N}
$$

with @@R_m@@ an auxiliary mode (type Link) and @@I_n@@ a physical mode (type Site). Then the partial Kronecker product is defined as

$$
\mathbf{C} = \mathbf{A} \boxtimes \mathbf{B} \in \mathbb{R}^{R_1 S_1 \times R_2 S_2 \times \cdots \times R_M S_M \times I_1 \times \cdots \times I_N}
$$

where the subtensors are given by

$$
\mathbf{C}(:, \dots, :, i_1, \dots, i_N) = \mathbf{A}(:, \dots, :, i_1, \dots, i_N) \otimes \mathbf{B}(:, \dots, :, i_1, \dots, i_N)
$$

If possible, I'd like to implement this in a general way, where it could be applied to any ITensor. I'm just getting hung up on how to do the indexing correctly and get the outer product that I need. Any help is greatly appreciated!

Answer 1

Hello,

If I understand the operation correctly, I think you could do it this way:

auto I1 = Index("I1",2);
auto I2 = Index("I2",2);
auto R1 = Index("R1",2);
auto R2 = Index("R2",2);
auto S1 = Index("S1",2);
auto S2 = Index("S2",2);

auto A = randomTensor(R1,R2,I1,I2);
auto B = randomTensor(S1,S2,prime(I1),prime(I2));
auto del1 = delta(I1,prime(I1),prime(I1,2));
auto del2 = delta(I2,prime(I2),prime(I2,2));

auto C = (A*del1*del2*B).noprime();

PrintData(C);

where here M=N=2, just for simplicity. For arbitrary M and N, it shouldn't be too hard to do it in terms of vectors of Indices and delta tensors.

Cheers,
Matt

Comments

Thanks so much, Matt! I think this is really close to what I need. Two quick things:

1. Can you explain what the double prime does in this case? I know it doesn't work without, but I'm not clear on why.

2. The dimensions of the tensor C in your example aren't quite consistent with the formula (should have product of R's and S's for each m). I fixed this with a combiner:

    auto T1 = combiner(R1, S1);
    auto T2 = combiner(R2, S2);
    auto C = ((A*del1*del2*B).noprime() * T1) * T2;

This is just a reshape as I understand it, so it shouldn't change the raw data of the tensor.

Thanks again! I never would have gotten this on my own; I'm still very new to tensors and ITensor.

---

Good questions!

1. The reason I made `delta(I1,prime(I1),prime(I1,2))` was just because we know for the sake of doing this contraction that we temporarily need three different versions of Index `I1`. `I1` is the version on ITensor `A`, `prime(I1)` is the version on ITensor `B`, and `prime(I1,2)` will be the version that will go on the output ITensor `C`. We need to make different versions to tell ITensor not to contract them with each other.

2. That's correct, I didn't include the combiner but that would be the right way to do it (note that you don't need the parentheses, the order of operations is left to right). In many cases, you may find that there is no need to combine the indices. It just depends on your use case.

---

Also, please note that in general contracting with a combiner may not be just a reshape, since the dimensions may have to be permuted first in order for them to be combined. But that is just an internal detail, you don't need to worry about that when using the combiner.

Cheers,
Matt

---

Thanks for the answer, Matt. One other comment I wanted to make is that, conceptually this may not be the fastest way to do this operation, since the delta tensors introduce some sparsity into the intermediate tensors. But if this part of your code is not a bottleneck in terms of taking up very much time, then Matt's solution ought to work well and help you to get the result you need. Please let us know, though, if this part is taking a very significant fraction of your code's overall running time.

Miles