How to create and add MPS

Question

Hello there,

Sorry for the newbie question (and if that has been answered elsewhere but I haven't found approaching solutions) but here it is:

I want to create 2 MPS in Julia and add them. The first question is that I would like to customise the core MPS. As such, my first and last core MPS are different from the others in the indices they hold. Is there a way to do that ? I.e. do I have to create a MPS with say N sites and then customise one by one like MPS1 = CustomTensor1; MPS[2] = CustomTensor2, etc... ? Then, I am not sure how would ITensor resolve the indexing then.

The other question is how to change the link dimension ? I heard about pluss* but not sure it was for Julia.

Finally, if I solve question 1 I assume I'd be able to sum them together. However, is that sum equivalent to the operations described in that article (bottom of page 7, top page 8) ?

Thanks for your help, that is really appreciated.

Kind regards,

Roland

Comments

Hi Roland, thanks for the question. We will need you to clarify a few things before we can answer.

1. when you say your first and last core MPS are different from the others, what others do you mean? And what do you mean by "core"? Finally, how are they different specifically? Maybe you could answer by saying it like this: I want to make an MPS where the first tensor of the MPS has n site indices and the last has m site indices, and you could tell us what n and m are i.e. how many indices.

2. when you ask how to change the link dimension, we have a function in Julia called `truncate!` which can truncate all the links of the MPS based on various truncation parameters. So is your question about the  existence of this function, or whether it would still work for your custom MPS or both?

3. regarding whether you can sum the MPS, it is a good question and we'd have to review our code in ITensor to check whether it is robust to having multiple sites. So we may have to help you with the first part of your question then chat more. As to how a sum is defined, it is defined in the usual way like in the paper (otherwise I'm not really sure what you are asking there).

Thanks,
Miles

---

Hi Miles,

Thanks for your answer. Let me clarify few points then:

1. In few places I've seen the terminology "core" to refer to MPS tensors at a given site. In my case, the tensors at the endpoints of my train have different sets of indices like in this https://arxiv.org/pdf/1802.07259.pdf Fig. 5. Essentially I am developing a MPS to polynomially approximate a function. Tensors at the endpoints have only bond dimensions equal to the degree of the approximator plus 1. All the other tensors have a physical bond of dimension 2 to accept a binary representation of a number.
2. The question may simply be rephrased: in the randomMPS case, I can set the linkdim (btw, if I do that, the tag is not updated to the dimension). It doesn't seem to be possible to do so while creating from the MPS method as (as far as I've seen) the linkdim is always 1. But maybe I missed something there.
 3. OK that sounds good. Happy to provide feedback/help as much as I can.

BW,
Roland

---

Not sure if you are aware, but we have a built-in MPS addition method:

https://itensor.github.io/ITensors.jl/stable/MPSandMPO.html#Base.:+-Tuple{MPS,MPS}

which works with MPS with any number of site indices. In order to add 2 or more MPS, they will need to have the same site indices. If they don't, you can use functions like `replaceinds` on the tensors of the MPS (https://itensor.github.io/ITensors.jl/stable/ITensorType.html#ITensors.replaceinds-Tuple{ITensor,Vararg{Any,N}%20where%20N}) to make the site indices match.

The current implementation in the Julia version is not quite what is in that paper, right now we use a density matrix approach where the local density matrices of each of the MPS are summed up and diagonalized to create the new MPS, similar to the algorithm here: https://tensornetwork.org/mps/algorithms/denmat_mpo_mps/ since it is simple to code and generalizes easily to summing multiple MPS as well as operations like MPS + MPO * MPS (not implemented yet but planned).

If you need very high precision for your MPS, this approach is limited to single precision accuracy for truncating the MPS, so let us know if it isn't sufficient.

When you say that you want an MPS with link dimensions that isn't one, but not randomMPS, do you want it to be filled with zeros? If so you could do something like:

s = siteinds("S=1/2", 4)
psi = randomMPS(s, 4)
for A in psi
  A .*= 0
end

Is it the case that you have a functional form for the elements of the MPS, so you know what values you will fill it with?

---

Dear @MattFishman,

Thanks the clarification.

OK then, it looks like I have implement my own MPS summation to match the article definition. I just wanted to make sure the ITensor MPS sum did the job before actually using it.

Yes, I know functionally the elements of the MPS. Since I follow what is in the article, I am at the point where I can construct several MPS but I need to perform some arithmetic, hence the questions about the sum.

OK let me proceed a bit further and come back to you if necessary.

Thanks again for the support.

BW,

Roland

---

Hello again,

Another stupid question maybe, can the block diagonalisation (i.e. the sum in the article above) be achieved using QNTensors ?

---

Hi Roland,
Regarding your earlier question about summing and whether ITensor implements the sum you want, you can always use the `inner` function to compute inner products between various MPS to check that a summation has succeeded, and to the accuracy you require. I think you get my point, but what I mean is if you write out the expression for the squared norm of the difference of two MPS, then if you expand that expression, you can compute that norm efficiently like inner(A,A) + inner(B,B) - 2*inner(A,B)  (here assuming all elements are real-valued. So if A is the result of the sum and B = C+D is the formal summation of two MPS C+D whose sum should give A, then you can expand the above expression to get a precise check of if the sum succeeded.

Regarding the block diagonalization, are you asking if you can construct an MPS whose entries are such that each tensor is block diagonal, and formed from two other MPS such that the resulting MPS equals the sum of those two MPS? If so, do you want to do this in order to comptue the sum? Or do you need this precise block-diagonal form for some other purposes? (I ask because if it's just to compute the sum, there are multiple ways to do that and not all necessarily involve constructing the block-diagonal form explicitly).

Best,
Miles

---

Hi Miles,

Thank you for your answer.

Actually yes, what I am after now is to implement the sum of two MPS as in the article I am following (ref above). In there, the sum of two MPS is given by constructing a resulting third whose entries are block diagonal tensors formed by putting the corresponding tensors from the other two. in short, I want to do M3[i] = M1[i] + M2[i] at a site i where M3[i] = diag(M1[i], M2[i]) where diag means block diagonal form. This is a preliminary step towards compression next.

If such a method is not present in ITensor yet, I am planning to implement my own by using multiple dispatch (if my basic understanding of Julia is right). But any other and better suggestion on how to do it is appreciated.

Thanks !

BW,

Roland

---

Hi Roland, it sounds like your goal is to sum two MPS. If so, then are you ok with the summation being carried out using a different technique from the one you read about as long as the result is the same? The technique we currently use does not explicitly form the block-diagonal MPS but gives the same final MPS as if you had done that and then compressed it.

Miles

---

Hi Miles,

Oh OK, yes sure then. Happy to go for it really. Part of that discussion was also for me to understand if there were equivalent classes of arithmetic operations built-in ITensor that I could use in lieu of reimplementing the ones in the paper.

Thanks !

BW,

Roland

---

Hello again,

Sorry for the hassle really.

Another question. If I want to contract an MPS over an input bit string, should I actually convert that bit string into an MPO to use the `contract` method ?

Thanks for your help

BW,

Roland

---

Hello,

Another (hopefully last) question, I am generating `randomMPS` like so:

```
 sites
4-element Array{Index{Int64},1}:
 (dim=2|id=138|"S=1/2,Site,n=1")
 (dim=2|id=322|"S=1/2,Site,n=2")
 (dim=2|id=961|"S=1/2,Site,n=3")
 (dim=2|id=555|"S=1/2,Site,n=4")

psi = randomMPS(sites, 2)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=2|id=139|"Link,l=1")
[2] IndexSet{3} (dim=2|id=139|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=2|id=979|"Link,l=2")
[3] IndexSet{3} (dim=2|id=979|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=869|"Link,l=3")
[4] IndexSet{2} (dim=2|id=869|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")
```
But if I increase `linkdim` it looks like the last bond dimension stays to 2

```
julia> psi = randomMPS(sites, 3)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=3|id=896|"Link,l=1")
[2] IndexSet{3} (dim=3|id=896|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=3|id=732|"Link,l=2")
[3] IndexSet{3} (dim=3|id=732|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=271|"Link,l=3")
[4] IndexSet{2} (dim=2|id=271|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")


julia> psi = randomMPS(sites, 4)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=4|id=261|"Link,l=1")
[2] IndexSet{3} (dim=4|id=261|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=4|id=677|"Link,l=2")
[3] IndexSet{3} (dim=4|id=677|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=418|"Link,l=3")
[4] IndexSet{2} (dim=2|id=418|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")


julia> psi = randomMPS(sites, 5)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=5|id=893|"Link,l=1")
[2] IndexSet{3} (dim=5|id=893|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=4|id=986|"Link,l=2")
[3] IndexSet{3} (dim=4|id=986|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=62|"Link,l=3")
[4] IndexSet{2} (dim=2|id=62|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")


julia> psi = randomMPS(sites, 6)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=6|id=226|"Link,l=1")
[2] IndexSet{3} (dim=6|id=226|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=4|id=841|"Link,l=2")
[3] IndexSet{3} (dim=4|id=841|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=562|"Link,l=3")
[4] IndexSet{2} (dim=2|id=562|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")


julia> psi = randomMPS(sites, 7)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=7|id=5|"Link,l=1")
[2] IndexSet{3} (dim=7|id=5|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=4|id=433|"Link,l=2")
[3] IndexSet{3} (dim=4|id=433|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=586|"Link,l=3")
[4] IndexSet{2} (dim=2|id=586|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")


julia> psi = randomMPS(sites, 8)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=8|id=120|"Link,l=1")
[2] IndexSet{3} (dim=8|id=120|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=4|id=831|"Link,l=2")
[3] IndexSet{3} (dim=4|id=831|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=292|"Link,l=3")
[4] IndexSet{2} (dim=2|id=292|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")


julia> psi = randomMPS(sites, 9)
MPS
[1] IndexSet{2} (dim=2|id=138|"S=1/2,Site,n=1") (dim=8|id=845|"Link,l=1")
[2] IndexSet{3} (dim=8|id=845|"Link,l=1") (dim=2|id=322|"S=1/2,Site,n=2") (dim=4|id=910|"Link,l=2")
[3] IndexSet{3} (dim=4|id=910|"Link,l=2") (dim=2|id=961|"S=1/2,Site,n=3") (dim=2|id=214|"Link,l=3")
[4] IndexSet{2} (dim=2|id=214|"Link,l=3") (dim=2|id=555|"S=1/2,Site,n=4")
```

Is there a limit or something ?

Thanks for your help

BW,

Roland

---

To contract with a bit string you can just create a `productMPS`, i.e.:

s = siteinds("S=1/2", 4)
psi = randomMPS(s, 4)
psi2 = productMPS(s, [1, 2, 1, 2]) # bit string but with 1-based indexing
inner(psi, psi2)

Indeed, there is an exact upper bound that can be imposed on the link dimensions at the edges of the MPS:

julia> s = siteinds("S=1/2", 10);

julia> linkdims(truncate(randomMPS(s, 5); cutoff = 1e-10))
9-element Vector{Int64}:
 2
 4
 5
 5
 5
 5
 5
 4
 2

This is true because the bond dimension is the rank of the matrix formed by reshaping the state into a matrix, for example the first dimension is the rank of the matrix:

psi_{s1, (s2,s3,s4,s5,...)}

which is always the dimension of the first site index, in this case 2. This may be inconvenient for analytically constructing an MPS, you could use the following function instead:

function emptyMPS(sites::Vector{<:Index}, linkdim = 1)
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(linkdim, "Link,l=$ii") for ii=1:N-1]
  for ii in eachindex(sites)
    s = sites[ii]
    if ii == 1
      v[ii] = emptyITensor(l[ii], s)
    elseif ii == N
      v[ii] = emptyITensor(l[ii-1], s)
    else
      v[ii] = emptyITensor(l[ii-1], s, l[ii])
    end
  end
  return MPS(v)
end

s = siteinds("S=1/2", 10)
psi = emptyMPS(s, 3) # make a bond dimension 3 MPS

In fact this would be useful functionality to have in ITensors.jl.

---

Thanks @MattFishman

That is really helpful.

A last point maybe and I think I'm good to go :)

If I set a MPS site to be any constructed tensor, say in your example above if I set psi[1] = randomTensor(Index(3)), I mess up with the indexing recognition system. Is there a way to assign a site to a tensor and respect the bond structure of the MPS ?

Thanks really, this is greatly appreciated.

BW,

Roland

---

It depends on what you want to do, but in general you would have to first extract the indices that the MPS already has, such as:

s = siteinds("S=1/2", 10)
psi = emptyMPS(s, 3)

# Want to modify MPS tensor 4 in some way
n = 4
sn = siteind(psi, n)
ln = linkind(psi, n-1)
rn = linkind(psi, n)
psi[n] = randomITensor(sn, ln, rn)

Note that you would have to be more careful with the Index arrow directions if you want your code to work with quantum number conservation.

---

Dear @Miles and @MattFishman

Thank you for your help, I finally went to the bottom of it now.

A maybe last quesion, you mentioned earlier than the truncation was limited to single precision. If I understand correctly, the sum operation performs such a truncation under the hood and you can control the precision using the cutoff keyword. Do you plan to increase this to double precision in a foreseeable future ?

BW,
Roland

---

The precision limitation is based on the algorithm we are using to perform the sum. It wouldn't be too hard to add a different algorithm that could get higher precision, but I haven't come across many cases where it is necessary yet. If you find the current algorithm is insufficient for one of your use cases let us know.

---

OK, sounds good to me.

You can consider this thread closed.

Thanks again for your precious help.

Roland

---

Hello again,

I am afraid I have to ask few more MPS questions :)

I am now trying to implement a gate extractor from the MPS I generate following the procedure in this article https://arxiv.org/pdf/1908.07958.pdf. Specifically the Fig. 1 a) at the top of page 2. This is realised by creating an MPD (D for disentangler) out of the MPS. It turns out that this MPD is constructed by adding new indices to the MPS. Hence the questions:

- Is the MPD equivalent to the MPO in ITensor ?
- How can I achieve the above construction ? (adding a vector of indices to a MPS) My naive solution would be to extract the indice set for each MPS site and add an index accordingly but maybe there is a more efficient built-in method that does it.

Thanks for your help !

BW,
Roland

---

Yes, the MPD would be represented as an MPO in ITensor.

Here would be some steps to add on a site index to an MPS:

s = siteinds("S=1/2", 4)
psi = randomMPS(s, 2)
Umps = copy(psi)
Umps .= Umps .* setelt.(prime.(dag.(s)) .=> 1)
U = convert(MPO, Umps)

Perhaps you are aware, but it is not quite as straightforward as just adding an additional index to the MPS to get the MPD/MPO (since the above will be filled with zeros for the new states that get introduced). You will have to follow along the arguments in Eq. 5-9, and make an orthonormal basis of the gauged MPS tensors and extend that basis by computing the kernel/null space of the MPS basis. This can be done in a variety of ways such as with a full QR decomposition or a full SVD, but unfortunately those kinds of tools aren't readily available in ITensors.jl right now so you'll have to put them together yourself (which shouldn't be too hard).

---

OK thanks a lot @MattFishman.

Is it possible to get back to you for some further advice in case ?

---

Sure, no problem. If you put it into a nice form, it would be a good algorithm to make widely available, such as in conjunction with the PastaQ.jl package (https://github.com/GTorlai/PastaQ.jl).

---

Hi @MattFishman,

Just a follow-up question on the MPS/MPD discussion above. I understand how to construct the MPD/MPO from the MPS by filling the MPD tensor elements with the corresponding MPS tensor elements. But I still fail to practically understand how to get the vectors from the kernel of the MPS tensors. If I get correctly what you suggested above, getting the MPD tensors from MPS tensors can be done by
1. Left-orthonormalising the MPS
2. Compute a full decomposition (SVD or QR) from which the matrix U is what I am looking for.
However, SVD is implemented in ITensors.jl ? Or is it not the full SVD ?

Thanks for your time and help.

BW,

Roland

---

The full SVD/QR is not currently implemented in ITensors.jl, so you will have to do something custom for your case. The Julia LinearAlgebra function `nullspace` would give you what you need (https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.nullspace). Here is a small script that gets the null space/kernel of an orthogonalized MPS tensor:

using ITensors
using LinearAlgebra

function LinearAlgebra.nullspace(A::ITensor, left_inds)
  Cl = combiner(left_inds)
  cl = combinedind(Cl)
  AC = A * Cl
  n = nullspace(array(AC))
  ln = Index(size(n, 2), "null")
  NC = itensor(nullspace(array(AC)), cl, ln)
  return NC * dag(Cl)
end

N = 4
s = siteinds("S=1/2", N)
psi = randomMPS(s, 2)
orthogonalize!(psi, N)

n = 2
A = psi[n]
left_inds = (linkind(psi, n-1), siteind(psi, n))
N = nullspace(A, left_inds)

@show N * A

The nullspace(::ITensor, ...) function simply wraps the Julia one that works on the underlying Array data of the ITensor. This is similar to how we would implement it internally (I have been meaning to add a nullspace(::ITensor, ...) function but it is a bit more complicated when there are QNs involved and I haven't had the time). Note you should check it does the correct thing for complex numbers (I didn't test it for that), and also note it won't work for QN conserving ITensors/MPS.

---

Hi @MattFishman,

Ok thanks I get it. It works fine as my tensors are only real-valued.

Last point that puzzles me a little is when computing the nullspace of end tensors. In a simple case describe in the article, these are 2x2 matrices and if the nullspace is non-empty, I should get 2x1 vectors. However, in most cases, I get 2x0 from LinearAlgebra.nullspace. I tried using rtol=2 but the returned matrix doesn't nullify the tensor. Any idea maybe ?

The other point, I am not sure what is meant by "Then choose (d^2 −d) orthonormal vectors in the kernel of A[n]". From the above procedure, I can compute the kernel of As which are 2x2x2 tensors but when constructing G, filling the elements for the added index should be done by choosing a vector in the kernel of A (as stated above). However, is this choice arbitrary ? What if the kernel is empty ? Sorry for the slightly off question. I am still trying to understand the whole process as being quite new to it.

Thanks a lot.

Kind regards,
Roland

---

The only case when there are 2x2 matrices is the first tensor at the edge of the MPS (tensor A[1] of Fig. 1 (a)), which indeed will always have a null space of size 2x0 (keep in mind because of the gauge it is already an orthogonal matrix, which always has a null space of size 0). So there will be no subspace to expand by in that case (and no index added on to the first tensor in the MPD), which is in Fig. 1(a).

For tensor A[2] to A[N-1], as you say the tensors are size 2x2x2 (dxdxd) which get reshaped into matrices of size d^2xd. Because you pick the orthogonal gauge with the gauge center at the end of the MPS (site N), the columns of this matrix are orthogonal to each other and therefore linearly independent. So the (left) nullspace is always the maximum it can be, i.e. (d^2-d) = 2 orthonormal vectors of length d^2 = 4. The left null space of that matrix cannot be empty since it is a rectangular matrix (the columns don't span the entire space).

---

Hi @MattFishman,

Thank you for you answer, that's clear really.

I am now able to compute the nullspace and place the resulting vectors in G. However, I can't seem to fullfil orthogonality conditions provided by Eq. (8). Any ideas maybe ? I tried several options without much success tbh.

Thanks again !

Roland

---

The two options are that you are putting the null space vectors into G incorrectly, or your MPS tensors weren't properly orthonormalized to begin with. Are you sure your starting MPS fulfills Eq. 2-4? Something you can do is check Eq. 8 component by component, i.e. if you fix `i = 0` it should reduce to Eq. 2-4, and then if you fix `i = 1` it should reduce to the orthonormality condition of the null space. Both of those should be "trivially" satisfied.

---

Dear @MattFishman and @Miles,

I guess I manage to successfully achieve what I wanted to do now.

I wanted to thank you once again for your amazing support and help. All this wouldn't have been possible without it.

I guess I will have no more comments there. If there is anything I can help with please let me know.

Kind regards,

Roland

---

Hi Roland,
Glad to hear you achieved what you were trying to do! (Thanks Matt for helping so much on this thread.)

Since you asked about helping with things, there are a few ways you can contribute:
(1) answering questions by others on this forum
(2) submitting improvements to the documentation, such as code examples (full code examples or small code example snippets after documentation of specific functions).
(3) If you find there are some specific features you've made which would be helpful for others, please consider making a pull request on our Github for us to include them in the library. If you do want to do this, though, please contact us first to ask if the contributions are in line with the overall design of the library.

Best,
Miles