Partially contracting over an index

ITensor Julia Questions
1

Hi!

The short question statement:
When contracting two tensors over a certain index, is there a way of specifying that the sum should not go over all values of the index, but over a subset of them?

The longer version:
In general, for two matrices A and B, which are of dimensions d_1 \times d_2 and d_2 \times d_3 respectively, we would have the product (contraction):

C_{i,j} = [A \cdot B]_{i,j} = \sum_{k=1}^{d_2} A_{i,k} B_{k,j}.

Is there a way of specifying that the sum should not go over the full range of k from 1 to d_2? Like only from k_{\text{init}} to d_2?
I tried to look into the implementation of contract(A::ITensor, B::ITensor) but was not succesful.
What I want to do would be the equivalent of the partial contraction in the example below.

A = rand(3, 5)                              # 3×5 Matrix{Float64}
B = rand(5, 4)                              # 5×4 Matrix{Float64}
# full contraction
C = A * B                                   # 3×4 Matrix{Float64}
# partial contraction
C_red = A[:, 3:end] * B[3:end, :]           # 3×4 Matrix{Float64}

Of course there are some naïve ways of getting that:

  1. I could create a copy of A (or B) where the first k_{\text{init}}-1 columns (rows) are replaced with zeros, then the result would be the same.
  2. I could create to new tensors which are the “cropping” of A and B and contract those, but that would imply creating a new Index object for the smaller index, creating two new tensors, and copying all the values from A and B to these smaller tensors, to then contract them.

Here is some example code for both, but it feels overly complicated (and probably very inefficient). What would be the most “ITensor-like” way of doing that?

using ITensors

i = Index(3)
j = Index(4)
k = Index(5)
A = randomITensor(i, k)                # ITensor ord=2 (dim=3|id=aaa) (dim=5|id=bbb)
B = randomITensor(k, j)                # ITensor ord=2 (dim=5|id=bbb) (dim=4|id=ccc)
# full contraction
C = A * B                              # ITensor ord=2 (dim=3|id=aaa) (dim=4|id=ccc)
# option 1
A_red = deepcopy(A)
for ii = 1:dim(i) 
    for kk = 1:2
        A_red[i=>ii, k=>kk] = 0
    end
end
C_red = A_red * B
# option 2
k_red = Index(3)
A_red = ITensor(i, k_red)
B_red = ITensor(k_red, j)
for ii = 1:dim(i) 
    for kk = 1:dim(k_red)
        A_red[i=>ii, k_red=>kk] = A[i=>ii, k=>kk]
    end
end
for jj = 1:dim(j) 
    for kk = 1:dim(k_red)
        B_red[j=>jj, k_red=>kk] = B[j=>jj, k=>kk]
    end
end
C_red = A_red * B_red
2

I think the best way right now would be to convert the ITensors to Julia Arrays, slice the Julia Arrays, and then make new ITensor indices of the size of those Julia Arrays and convert back to ITensors. That would only work with dense arrays, however (i.e. the procedure would be more complicated for block sparse/QN-conserving ITensors, which is part of the reason why we don’t support this kind of slicing operation yet, though we plan to).

3

Thanks! I’ll do that then