# How to give local measurement for MPS and sample the measurement bitstrings from the calculated probabilities.

I do a tebd time evolution with MPS state. I want to make a local measurement for MPS state during each time slot. I do not know how to take local measurement on MPS state. I want to obtain the final samples, i.e. if I give X\otimesY\otimesZ\otimesX\otimesY… measurements over N qubits, then i will obtain a bitstring like 01010101… And i will give many shots over the specified measurement setting. So, how to implement this based on itensor. Could some one give me a snipt code? Thanks.

Hi, Xiao!

Do you want the expectation value of certain operator string for your MPS? Then it can be directly calculated in ITensor.
Let’s take a domain wall state for illustration. So it looks like |1111100000> where 1 & 0 correspond to “Up” & “Down” state for spin=1/2.
The state is initialised below for 10 sites.

using ITensors

N=10
sites=siteinds("S=1/2",N)
psi = productMPS(sites, n -> n>5 ? "Up" : "Dn")


Now for example lets take an operator string that has Sz_j \otimes Sz_{(j+1)} .
This operator will help us to detect the domain wall position. You can choose any other string of local operator in place of Sz .This operator will give you < Sz_j Sz_{(j+1)}> when applied on sites j,j+1.

A function that measures the expectation value is given below.

function dw_pos(psi,n)
psi=psi*(1.0/norm(psi))
sites=siteinds(psi)    #Get the site indices of MPS to construct suitable operator
psi = ITensors.orthogonalize(psi,n)   # Improves accuracy
ket=psi[n]*psi[n+1]
# Here you need to put your desired operator. sites[n] implies on which site to apply the operator"Sz"
zzop = op("Sz",sites[n])*op("Sz",sites[n+1])
bra=(dag(prime(ket,"Site")))

dw=scalar(bra*zzop*ket)
return dw
end


To test this,

for j=1:length(psi)-1

println("Sz_",j,"_Sz_",j+1,"   ",dw_pos(psi,j))
end


You will get -0.25 @ bond 5-6 indicating anti parallel spin . All other bonds will give you 0.25.

Notice this line,

    zzop = op("Sz",sites[n])*op("Sz",sites[n+1])



If you want something like X\otimes Y \otimes Z \otimes X on qubits 1,2,3,4 respectively,just put

observ = op("Sx",sites[1])*op("Sy",sites[2])*op("Sy",sites[3])*op("Sx",sites[4])



The rest is same.You don’t need to put the identities,they are automatically assumed.

I am not sure this is what you are asking for. Hope this helps.

Have a good day.

Hello, sandipanmanna, first of all, thank you very much for your carefule reply. You indeed give me a good turorial on how to calculate the expecation of local observables. That is very useful.

My further questions are more concrete. In quantum computation, we need to simulate the measurement bitstrings from measuring the final quantum state. For instance, suppose the final evolved MPS state is denoted as \rho, and we use the unitary operator U = X\otimesY\otimes Z\otimes Z … to tranforme the quantum state, denoted as \rho’ = U\rhoU^dagger. Then, we calculate the probability of projecting to 0 for each site i, i.e. P_0(i) = \langle 0 | \rho’ |0\rangle. The probability of measuring 0 once is calculated, then the probability of obtaining 1 can also be obtaineed. Then, we sample the 0 or 1 according to the probability we have obtained. And we will repeate the measurement many times to obtain a series of measurement records. The bitstring length is equal to the number of sites (#qubits). The core task is to calculate the probability of projecting the final transformed state to 0 or 1 state for each site and then according to the probability to sample the bitstrings.

Since current density matrix simulator in classical computer cannot simulate large size of quantum systems, I want to try MPS and its evolution.

Thanks again.

Currently, there is another choice. I can save the MPS state into the disk, and then read the file from disk. Then, I can process the read MPS information to accomplish the sampling process. But I do not know how to save the MPS state in the fixed format of the following,

Each Matrix product state is described by two files:

index.txt lists for each of the N tensors in the MPS, its rank, the number of non-zero elements and the dimension of each index:

rank Number non-zero elements dim(1) … dim(rank)

tensor.txt contains all the non-zero values of each tensor T (from 1 to N), with their index configuration:

index1 index2 T(index1,index2) … …

index1 index2 index3 T(index1,index2,index3) …

The order in which the indices appear for each tensor is:

for the left-most tensor: (right,physical)
for “bulk" tensors : (right, physical, left)
for the right-most tensor: (physical, left).

If you can teach me how to save the MPS into two files, index.txt and tensor.txt and I will be very appreaciated. The index.txt and tensor.txt files can be seen in Github source code example. POVM_GENMODEL/data/training_data_1D_TFIM_model/Matrix_product_state at master · carrasqu/POVM_GENMODEL · GitHub.

Thanks.

Hi Xiao,
To measure an MPS in the Z basis (computational basis), you can use the function sample! provided by ITensor:

bits = sample!(psi)
@show bits


You can call sample! many times on the same wavefunction to get multiple samples from it.

If you want to do the measurement in a different basis, you can first make a copy of the wavefunction:

psic = copy(psi)


then apply the X or H or other operators to it as desired following this code example:
https://itensor.github.io/ITensors.jl/stable/examples/MPSandMPO.html#Applying-a-Single-site-Operator-to-an-MPS
then call sample! to measure in the (new) computational basis.

Miles

P.S. @sandipanmanna thanks for posting the helpful code about computing expectation values

1 Like

Thanks very very much, miles.

1 Like

hello, miles. I tried your tutorial. That is what i need. I have another question that if i want to make a computation basis rotation for N sites (qubits) with U= X \otimes Y \otimes Z … \otimes Y, then how to achieve this. One way is to separately add the local operator into the MPS and update it. The process iterates N times and there will be N updates.

I do not know whether orthogonalizing the mps does affect the final results during applying the local gates to each MPS site.

I suppose there will be a convenient way to simultaneously apply N local gates to all sites and I can sample it. Can you give some suggestions?

Thanks.

Hi Xiao, if I understand correctly, are you asking how to apply an operator like

X_1 \otimes Y_2 \otimes Z_3 \cdots \otimes Y_N

to an MPS?

And also whether the process of orthogonalizing the MPS changes the final results in some way?

Please let me know if this is your question (it’s a good question) and we can discuss further.

Yes, I’m trying to write the code of N local tensor product. It seems that when N is 20 or larger, the warnning as the following

Contraction resulted in ITensor with 20 indices, which is greater than or equal to the ITensor order warning threshold 14. You can modify the threshold with macros like @set_warn_order N, @reset_warn_order, and @disable_warn_order or functions like ITensors.set_warn_order(N::Int), ITensors.reset_warn_order(), and ITensors.disable_warn_order().

will be given. I do not know wether this warning will affect the performance. That is the first question.

The code looks like:

ops = [op(string(paulistring[i]),sites[i]) for i=1:N]
unitary = reduce(*, ops) ## give the tensor product of N ops
new_psi = unitary*psi
noprime!(new_psi)
bits = sample!(new_psi)


To reply your question, firstly, I can iteratively give one local operator to each one site until the last site is given the local operator X or Y or Z. As the tutorial suggest (apply one-site local operator to MPS), I need to update the MPS each time i apply the local operator. Since the target is to do sampling from MPS state with given pauli transformations XYZYX…Y, I will do it as the following code:

for j =1:10000
define a paulistring  with length N
psi_ = copy(psi)  ## here deepcopy is necessary?
for i=1:N
local_op = op(string(paulistring[i]), sites[i])
psi_ = local_op*psi_  ## here may have errors
noprime!(psi_)
end

bits = sample!(psi_)
println(bits)
end


Here psi i want to sample is a deepcopy of another psi, and the latter may be used to evolution. As long as i make a deepcopy of time evolution psi, sampling process will not affect the psi in time evolution. My question is that when i give local gates to each sites one by one, if i do not give orthogonalization of the psi during each local gate, whether the final sampling results are affected?

Thanks.

I see thanks for explaining some more. Here are some responses and thoughts that may help you:

• you should definitely not try to make the entire operator as a single tensor (i.e. by calling reduce(*,ops)). This will result in a tensor with 2N indices and will have 2^N elements so will scale exponentially and cost a lot of resources.
• instead the approach in your second block of code is much closer to what you should do, which is to apply the Pauli operators one at a time to your MPS. This will be equivalent to applying the whole tensor product operator.
• to apply the Pauli operators, simply repeat the code in the Code Example (https://itensor.github.io/ITensors.jl/stable/examples/MPSandMPO.html#Applying-a-Single-site-Operator-to-an-MPS) multiple times.
• Please pay close attention to the details of the code in the example. If you look at it, you’ll see that when applying the operator you have to apply it to a specific MPS tensor, so like local_op*psi_[j] (similar to the line newA = G * psi[3] in the example)
• For these single-site operators, since no truncation of the MPS happens, you can ignore any details of the orthogonality of the MPS while applying the operators. After all the operators are applied, you can call orthogonalize!(psi,1) to restore the orthogonal form of the MPS, though also if you just directly call sample!(psi) it internally orthogonalizes the MPS (as mentioned in the docstring for sample!) so technically you don’t need to do this yourself.
• You shouldn’t have to call deepcopy here, as I believe copy would do but it’s always safe to call deepcopy so can’t hurt if you aren’t sure.

Hope that helps & let me know if you have more questions.

1 Like

Glad it helped. So basically the lines in your for loop should look like this:

psi_[i] = local_op*psi_[i]
noprime!(psi_[i])


The rest of that code looks good to me

1 Like

Hahaha, thanks for you check.

1 Like