Canonical MPS representation given an initial state in Julia

Question

Hi,

I was wondering if there is any inbuilt script to automatically get the canonical MPS representation from a given state using Schmidt decomposition.

If not, how would one go about doing this?

Comments

Hi Arnab, Matt or I will post a longer answer soon, but one question to clarify: by canonical MPS representation do you mean the “gamma-lambda” form of Vidal? Or do you mean a different one? Sometimes the forms where most of the tensors are left- or right-orthogonal are called canonical too.

If you are asking about the gamma-lambda form, the short answer is no we don’t have a function you can call for this at present. But there’s definitely some compact code we could help you write to do it. Also we do have a lot of facilities for making the left and right orthogonal forms of MPS already so if you just need those then we can help you there.

Miles

---

Hi Miles,

Maybe I can change my question a little. My objective here is the following:
Given the amplitudes of a state in some basis, I want to get the most efficient (minimum  bond dimension) MPS representation of it.

I don't really know if some procedure achieves that( my initial guess was the canonical form).

Best,
Arnab

---

Hi Arnab, I see. Yes that is a rather different thing from obtaining a canonical form of a given MPS, but it's a good question.

It may not be too important in terms of changing my answer, but let me ask: do you mean you will input *all* of the amplitudes of the state or just a sampling of some of them?

Miles

---

Hi Miles,

Yes, the state is *completely* known with all its amplitudes. Maybe for reference, let me say what needs to be done in this case is outlined in  Schollwoeck's review (https://arxiv.org/pdf/1008.3477.pdf) in Section 4.1.3: Decomposition of arbitrary quantum states into MPS (it uses different Schmidt decompositions).

---

Hi Miles,

Just saw this post here. May I ask if there's any function in the Julia version which one, say, could initialize an arbitrary MPS with left/right canonical form?

Say I want all to be right-canonical. My approach is that maybe I can orthogolize MPS to site 1, and then perform an additional SVD or QR decomposition on the first site?

Best,
Tianqi

---

Hi Tianqi,
So once you have obtained an MPS from any constructor (so like randomMPS, productMPS, or any other way of obtaining an MPS), then you can shift the orthogonality center using the function:
orthogonalize!(psi,j)
where psi is the MPS and j is the site number you want to be the ortho center. If you choose j=1, then the MPS will afterward be in right-orthogonal form (all tensors right-orthogonal). If you choose j=N then all tensors will be left-orthogonal.

Does that help?

Miles

---

Thanks a lot Miles! Indeed it is helpful for me to understand it better now. In fact, I want to have this kind of canonical form so that I can implement myself a \Gamma-\Lambda notation of MPS for DMRG and TEBD.

Answer 1

Is this what you are looking for?

using ITensors
N = 4
s = siteinds("S=1/2", N)
A = randomITensor(s...)
psi = MPS(A, s; cutoff = 1e-10)
@show maxlinkdim(psi)
@show isapprox(A, prod(psi))

Notice that prod is a general Julia function for multiplying a collection together, so in the case of an MPS it contracts all of the tensors together into a single ITensor to reproduce the (approximation) of the original state. You can also set a maximum dimension with the keyword argument maxdim.

-Matt

Comments

Hi Matt,

Sorry for the late response. I am not quite sure if this does the trick. Let me be a little clearer about my question: I want to feed in the amplitudes of a state (so a bunch of scalars wrt. some basis) and get back an (approximate) MPS description with the smallest bond dimension. So, for example, if I feed in the amplitudes of the cluster state in 1D (in the computational basis) , I want to get back an MPS description with bond dimension 2.

Best,
Arnab

---

Hi Arnab,

I'm probably misunderstanding, but I believe that is covered by the code above. The amplitudes of the initial state are stored in the ITensor `A`, for example a single amplitude would be `A[s[1]=>1, s[2]=>2, s[3]=>1, s[4]=>2]` for the basis state `|1212>`.

Maybe if the above doesn't work, could you give a simple example and what you have tried that hasn't worked?

-Matt

---

Hi Matt,

for the states I am concerned with,  the output is the following:

maxlinkdim(psi_optimised) = 4
isapprox(A, prod(psi_optimised)) = false

I don't know how to interpret that.  The link dimension is pretty small but the approximation seemingly fails. I don't mind the link dimension going up as long as it is the smallest to make the approximation work.

Best,
Arnab

---

What is the output of `norm(A - prod(psi_optimised))`? It may just be that the default cutoff for the `isapprox` function is too strict, but they may be approximately equal up to some floating point amount.

---

You can also adjust the cutoff to something else, like `cutoff = 0` (or just don't specify a cutoff in which case it won't truncate). If you are still finding that the MPS state is not approximately equal to the original ITensor, a minimal example of what you are trying would be helpful so we can help debug what is going on.

---

I set the cut-off to 1e-5 which (I  think) is not terribly strict.

norm(A - prod(psi_optimised)) gives me 0.006543262674151699.

---

I set cut-off = 0 and it seems to work. I think that will do for now. Thanks a ton Matt for bearing with me!

---

No problem, glad that functionality is being used and seems to be working (I wrote it mostly for internal usage in the package but I think it is useful more generally).

Getting a global error of around 6e-3 using a cutoff of 1e-5 doesn't sounds so unreasonable, since the cutoff is used at every bond of the MPS (so the total error would be something like the sum of the truncation error at each bond). Depending on the size of your MPS, the global error measured by `norm(A - prod(psi_optimised))` can be larger than the cutoff you use, as you were seeing.

Also note that the `isapprox` function has a few keyword arguments for setting the tolerance for the approximation, for example `isapprox(A, prod(psi); rtol = 1e-5)` (where `rtol` sets a relative tolerance scaled by the norms of the input tensors). A reasonable choice for `rtol` may be something like `isapprox(A, prod(psi); rtol = cutoff * N)` where `cutoff` is the cutoff you set and `N` is the system size. You would have to play around with this, of course. Doing floating point comparisons is always a bit tricky.

---

Thank you Matt! The keywords in the 'isapprox' function were of great help when running the code for large chain. A somewhat unrelated (and perhaps unimportant) question: Is there an intelligent way of assigning values to an ITensor? For example: I want to assign values like
A[sites[i]=>(arr[i]+1)]= constant- for i=1 to N; where A is an ITensor; arr is a 3 element array.
I can do this brute-force by A[sites[1]=>(arr[1]+1),sites[2]=>(arr[2]+1),sites[3]=>(arr[3]+1)]= constant for N=3 but for straight-forward extension to {for i=1:N do A[sites[i]=>(arr[i]+1)...]= constant end } doesn't work.

---

Maybe:

    A[[sites[i] => (arr[i]+1) for i in 1:N]...] = ...