Converting ITensors to matrices in Julia

Question

Hi,

Sorry if this has been answered already but I am trying to do some post-processing with the MPS Tensor components and it would be really helpful to just convert them into regular matrices. Specifically, I have some MPS tensor "T" on N sites on a spin-1/2 chain. I can grab the two local components of the MPS using commands like -

using ITensors
N =10
sites = siteinds("S=1/2", N)
T = randomMPS(sites, 20) 
T5_1=T[5]*setelt(sites[5]=>1)
T5_2=T[5]*setelt(sites[5]=>2)

However, I would like to be able to perform matrix multiplications like (T51).\rho.(T51)^{\dagger} (where \rho is a suitable density matrix) and for this it would be best to have them converted into matrices.

PS: I know that the other way of doing this is to encode the density matrix into an ITensor object and change the link indices suitably to get the multiplication to work out.

Answer 1

The short answer is that you can use the Array or Matrix function:

i = Index(2, "i")
T = randomITensor(i', dag(i))
A = Array(A, i', i)
# Same thing, for an order-2 ITensor
M = Matrix(A, i', i)

where the indices are used to specify the memory ordering of the Array/Matrix. You can also use the lower case function A = array(A) which uses the current index ordering of the ITensor for the memory ordering of the Array, and makes the Array a view of the ITensor when possible. However, generally I would recommend working with ITensors when possible since it can help minimize mistakes with connecting the wrong indices together.

I'll also point out that we are planning on adding more functions to help with these conversions: https://github.com/ITensor/ITensors.jl/issues/374

Comments

Hi Matt,

Thank you for your answer! I will try this out today. I was actually first attempting to do everything in ITensor format but I am not very familiar with the underlying data-structure to do all the things I need. For example, I am a little puzzled that even when I make two links match the contraction doesn't happen because of different ids?(where rho has indices i,j)

M = T[finalsite] * setelt(sites[finalsite] => 1) * delta(i, linkind(T, finalsite)) * rho = ITensor ord=4
Dim 1: (dim=6|id=738|"Link,fact")
Dim 2: (dim=7|id=618|"i")
Dim 3: (dim=7|id=63|"i")
Dim 4: (dim=7|id=846|"j")

---

It's difficult to answer since I don't know how you are constructing rho, but it looks like you have two different `i` indices that have different ids so you would have to fix that in some way. Could you show a minimal working example of the code you are running?

---

Taking your suggestion to heart, I got things to work by just working with ITensors. I am running a loop and that was changing the "i" index earlier . A somewhat unrelated question regarding ITensors: is there automatic ways to set the elements of a randomITensor to (i) all zero to start with? (ii) complex values in Julia? I presume it follows the C++ version closely?
https://www.itensor.org/docs.cgi?page=classes/itensor&vers=cppv3

---

If you use the simple constructor `ITensor(i', i)` it will make an ITensor filled with zeros. To make an ITensor filled with zeros with complex elements, you can use `ITensor(ComplexF64, i', i)`. Similarly, you can fill it with random complex elements with `randomITensor(ComplexF64, i', i)`.

---

The constructors are documented here: https://itensor.github.io/ITensors.jl/stable/ITensorType.html#Dense-Constructors-1 though they should probably have some simple examples to explain how they are used in common cases.

---

Hi Matt,
Thank you for your answers. The documentation was helpful to resolve my issue. Sorry for asking so many basic dumb questions. I am not sure what exactly is going wrong when I try to access the 2nd order components of the 3 dimensional tensors using "setelt". Here is a code that generates an ensemble of 30 random density matrices of dimension 5*5.
using ITensors
using QuantumInformation
using LinearAlgebra
function random_ensemble(size,en_size)
i = Index(size,"i")
j = Index(size,"j")
en = Index(en_size,"en")
rho = ITensor(i,j,en);
for en_index=1:en_size
h = HilbertSchmidtStates{1}(size)
ρ = rand(h);
for x=1:size
for y=1:size
rho[x,y,en_index]= ρ[x,y];
end
end
end
return rho
end
rho_new=random_ensemble(5,30)
However, when I call "rho_new*setelt(en=>5)" to  obtain the 5th element of the ensemble, I instead get an order 4 tensor with two "en" index.

---

Again, I don't think I have quite enough information to answer since I don't know where the Index `en` is coming from when you are trying `rho_new*setelt(en=>5)`. However, likely it is not the same as the one inside the function `random_ensemble`. Probably instead you want to write the function `random_ensemble` to accept an Index `en`, so then you would do:

    function random_ensemble(i::Index, en::Index)
        size = dim(i)
        j = settags(i, "j")
        en_size = dim(en)
        [...]
        return rho # rho was made using the input Index en
    end

    i = Index(5, "i")
    en = Index(30, "en")
    rho_new = random_ensemble(i, en)
    rho_new * setelt(en=>5) # This should now work since the `en` is the same one used to construct `rho_new`

It can take some getting used to working with ITensor indices, but a common pattern is that you will need to pass indices to functions that construct ITensors so then you have access to those indices outside of the functions.

---

I see. In my script when I now return en,rho in the function then rho*setelt(en=>5) indeed gives me the correct result. I will keep this in mind in the future.