hi, thanks for the great work on itensors! I recently started playing around with it and love its features, though I have to say I’m a bit puzzled by the way tensor indices are handled,
in particular, it seems to me that there’s no straightforward way of accessing and playing around with the indices of a given tensor - I understand that for a given tensor I can access those by calling inds(A) , but sometimes I’d like to be able to call an index with a given label, I thought the concept of ‘tags’ would be helpful there but I’m not sure what is the “recommended” way to use them in this context.
A silly example, if I do the SVD of a matrix A_{ij} and get a decomposition as U_{iu} S_{uv} Vd_{vj}, is there a way to access, eg. to relabel or prime it, the “u” index in the matrix U ?
The only way I could think of is to get both indices of U, and then “know” that “u” is the second one - but since labels can in principle float around (as long as they’re attached to the proper leg), I’m afraid that calling inds(U)[2] could give unpredictable results. Is there a way to call the index containing a “u” in its tag instead? Or say I want to extract the index of the physical leg of an MPS matrix, it would be the one containing a “Site” tag.
I was digging through the documentation and I saw the getfirst() function, but its syntax is not clear to me (in particular what the matching function should look like, but maybe that’s my bad since I’m new to julia as well) - is that the intended way to deal with this use case?
Good question. Our recommended way of accessing the u index is through the commonind function, for example:
i, j = Index.((2, 2))
A = randomITensor(i, j)
U, S, V = svd(A, i)
u = commonind(U, S)
Then you can use it like:
dag(U) * prime(U, u) # Check if `U` is an isometry
More generally, you can think about sets of ITensors as networks of tensors, and usually you can uniquely determine a certain index or set of indices through the network connectivity, such as which indices are unique to a certain tensor or shared between pairs of tensors (which is how the commonind function is being used above). I find that I can almost always just use functions like commonind[s] or uniqueind[s] to access the indices that I want, and that is generally our goal when designing ITensor library code to make it as general as possible (as opposed to relying too much on tags to access a specific index, since then that relies on a certain convention for the names of tags which is more fragile). Often to do this properly it is helpful to have a drawing of the tensor network so you can think visually about which indices are connected to which tensors, our visualization packages like https://github.com/ITensor/ITensors.jl/tree/main/ITensorVisualizationBase can help with that.
Thanks for the feedback!
I see, indeed I had missed these [unique|common]inds functions, those are definitely helpful!
But let’s say (sorry for coming up with awkward use cases) I want to compute some kind of weird “left” environment from two identical MPS, by contracting eg. an A^{[i]}_{LRp} tensor (here I labelled the virtual Left, Right and physical legs) with itself over the legs L and p. How would I manage to single out the left (or the right) legs here ? Even if I started for simplicity from one edge where I’d have only one physical and one virtual indices, if I prime a copy of A, then all indices would get primed, even the physical one (which is something I don’t want). Would there be a clever way here to isolate a specific index without using eg. inds(A)[1] here ?
Again, I realize this is a bit of a special scenario so probably not really your everyday use case, still I wonder if there’s some intended way of getting it in this case.
It’s a good question, it’s something that is really at the core of the philosophy of the ITensor library and Index system, and something we have discussed a lot internally when designing the library to come up with nice code patterns to handle those kinds of operations.
How I think about it is that you shouldn’t really think of a tensor in isolation, but as part of a network. The indices are therefore defined through how they connect to neighboring tensors. So if you have an MPS A, you can perform various contractions of a tensor of A like this:
n = 10
s = siteinds("S=1/2", n)
A = randomMPS(s)
i = 4 # Site we are interested in
# Get indices of `A[i]`
L = commonind(A[i - 1], A[i])
R = commonind(A[i], A[i + 1])
p = uniqueind(A[i], A[i - 1], A[i + 1]) # Could just use `p = s[i]` since we have `s` already
# Contract over `L` and `p`
dag(A[i]) * prime(A[i], R)
# Contract over `L` only
dag(A[i]) * prime(A[i], (R, p))
Note that we provide shorthand functions siteind and linkind, so you could do:
L = linkind(A, i - 1)
R = linkind(A, i)
p = siteind(A, i)
but those just use commonind and uniqueind internally.