# Canonical MPS representation given an initial state in Julia

+1 vote

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?

commented by (51.6k points)
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
commented by (530 points)
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
commented by (51.6k points)
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
commented by (530 points)
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).

selected by

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 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

commented by (530 points)
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
commented by (10.7k points)
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
commented by (530 points)
Hi Matt,

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

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
commented by (10.7k points)
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.
commented by (10.7k points)
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.
commented by (530 points)
edited by
I set the cut-off to 1e-5 which (I  think) is not terribly strict.

norm(A - prod(psi_optimised)) gives me 0.006543262674151699.
commented by (530 points)
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!
commented by (10.7k points)
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.
commented by (530 points)
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.
commented by (10.7k points)
Maybe:

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