Hi Miles,
I wonder whether there is a way to take an element-wise square root of a MPS? Namely, given a non-negative tensor \psi_{i1,\dots,iN}, I want to obtain a new tensor with the element \psi'_{i1,\dots,iN}= \sqrt{\psi_{i1,\dots,iN}}.
This is actually a follow-up question to my previous question: Compute the sum of norm of an MPO
Currently, I have some thoughts, but I am not sure whether this is the optimal/right way. Therefore, your comments would be greatly appreciated.
My ultimate goal is to compute the sum of the norm of a (pure state) density matrix \rho, i.e., \sum_{ij}|{\rho_{ij}}|.
The summation can actually be done by contracting with an âall-one tensor.â So the problem boils down to how to find the tensor |{\rho_{ij}}|.
Since in my problem, \rho comes from a pure state, i.e., \rho =|\psi\rangle\langle\psi|.
Assume |\psi\rangle = [a_1,a_2,\dots, a_i, \dots]^\intercal , clearly, |\rho_{ij}| = |{a_i}| |{a_j}|.
This is to say, if I can find a way to convert the MPS of |{\psi}\rangle to a |{\psi}'\rangle, where the âphase of each component is droppedâ, then I can just construct an outer product of |{\psi'}\rangle, and take the sum of all the elements.
Then, this problem boils down to how to âapply a norm to an MPS.â
(assume the new MPS can be represented by MPS, even if not, I still want to find an algorithm that will be exact if allowing exponential growth of bond dim, such that I can later do some truncation to approach the exact results)
In order to do so, the norm of a is actually just \sqrt{a_i a_i^*}, which can be further broken into 2 steps:
- Element-wise product of
psiandconj(psi), which we callpsi_abs2 - Element-wise square root of
psi_abs2.
The first step can be obtained by the element-wise product as follows
function elementwise_product(mps1::MPS, mps2::MPS;cutoff::Float64=1e-10,maxdim::Int=25,method::String="densitymatrix")
site_idx=siteinds(mps1)
mpo1=MPO(site_idx)
for i in 1:length(mps1)
mpo1[i]=mps1[i]* delta(site_idx[i],prime(site_idx[i],1),prime(site_idx[i],2))
end
mps_prod=apply(mpo1,prime(mps2;tags="Site"),method=method,cutoff=cutoff,maxdim=maxdim)
return noprime!(mps_prod)
end
following Sec. IIIB in https://journals.aps.org/prx/abstract/10.1103/PhysRevX.13.021015
My attempt for the second step is as follows:
a_0=S
c_0=S-1
a_{n+1}=a_n-a_n c_n/2
c_{n+1}=c_n^2(c_n-3)/4
Here, we take S to be the MPS we want to take the square of. Namely, S=psi_abs2, and then iterate a_i, until it reaches convergence (e.g., the update is smaller than a threshold)
Then the final tensor a_n is the square root of the tensor S.
The iterative method is from Methods of computing square roots - Wikipedia
Mathematically, this method converges given 0<S<3, which is naturally satisfied here as the amplitude of each component cannot be larger than 1.
Because each iteration rule here only involves âelement-wise productâ and âadd/minus a constant MPS,â it, in principle, can be efficiently performed.
My minimal implementation is
function sqrt_mps(mps::MPS;eps::Float64=1e-5,maxiter::Int=100,cutoff::Float64=1e-10,maxdim::Int=20)
L=siteinds(mps)
allone=CT.all_one_mps(siteinds(mps))
a=mps
c=mps-allone
diff_a = 100
while diff_a > eps && maxiter > 0
a_new = -(a,elementwise_product(a,c,cutoff=cutoff,maxdim=maxdim)/2,cutoff=cutoff,maxdim=maxdim)
truncate!(a_new)
c = elementwise_product(elementwise_product(c,c, cutoff=cutoff,maxdim=maxdim),-(c,3*allone,cutoff=cutoff,maxdim=maxdim), cutoff=cutoff,maxdim=maxdim)/4
truncate!(c)
diff_a_mps=a-a_new
diff_a = abs(inner(conj(diff_a_mps),diff_a_mps))
a=a_new
maxiter -= 1
end
return a,c
end
function all_one_mps(sites::Vector{Index{Int64}})
mps=MPS(sites)
for i in 1:length(mps)
mps[i]=ITensor([1,1],sites[i])
end
return mps
end
I experimented on several examples and found that for random MPS with bond dim ~6 and system size ~ 12. this can give you an approximation of the âelement-wise norm of an MPS.â
However, I do also find numerical stability to be an issueâ sometimes, when the amplitude of the wavefunction is too small, c_i could numerically diverge, and therefore ruin the next a_{i+1}.
Another thing is that I found the memory allocation is very large in my code, which I guess is due to the âelement-wise productâ, where it returns a new MPS each time. So, it also takes quite a long time to perform the iteration.
Do you think my method is feasible, or is there another better way?
In general, is it an open question to ask how to apply an arbitrary function f(x) to an MPS or QTT, and what is the function space such that the new MPS aftering the element-wise mapping is also representable by MPS?