Hi,
I’m a newcomer to ITensor, and I encountered some confusing behavior when trying to decompose and remake two tensors in an MPS. The steps I follow are:
- Create an MPS
psifrom an InitState defined on a BosonSite - Define
G = psi(1)*psi(2), and perform its SVD (I have checked thatU*S*V == G) - Create some new MPS defined on the same BosonSite, called
phi - Assign
Utophi(1)andS*Vtophi(2)
Since psi is normalized, I expect \text{tr}(S^2) to be 1 (which it is), which would make phi normalized as well. phi however isn’t normalized. If I explicitly set a MaxDim parameter in the svd, then phi is. Is there anything in the behavior of svd or psi.set() that I am missing? I have attached a version of the code below. Thanks!
#include<iostream>
#include "itensor/all.h"
using namespace itensor;
template <typename T>
double sumofschm2(T S)
{
double s,d = 0;
for (auto i : range1(minDim(S)))
{
s = elt(S,i,i);
d+=(s*s);
}
return d;
}
template <typename T>
Complex my_norm(T wf)
{
auto C = dag(wf(1))*wf(1);
C*=wf(2);
C*=dag(wf(2));
return eltC(C);
}
int main()
{
uint N = 7;
auto sites = Boson(2,{"MaxOcc=",N});
auto state = InitState(sites,"2");
auto psi = MPS(state);
auto G = psi(1)*psi(2);
//SVD without MaxDim argument
{
auto [U,S,V] = svd(G,inds(psi(1)));
print(U);
print(S); //Is of order 5x5, with one nonzero element (=1)
auto phi = MPS(sites);
phi.set(1,U);
phi.set(2,S*V);
//Print sum of squares of schmidt values
print('\n',"Sum of squares of schmidt Values: ",sumofschm2(S),'\n'); //Returns 1
//Manually computing the norm of psi
print("The norm of the wavefunction is: ",my_norm(phi),'\n'); //Returns 5
//U*S*V is equal to G
print("Norm of U*S*V - G: ",norm(U*S*V - G),'\n');
}
print("\nSVD with MaxDim\n");
//SVD with MaxDim argument
{
auto [U,S,V] = svd(G,inds(psi(1)),{"MaxDim",7}); //Even though MaxDim is greater than 5
print(U);
print(S); //Is of order 1x1, with one nonzero element (=1)
auto phi = MPS(sites);
phi.set(1,U);
phi.set(2,S*V);
//Print sum of squares of schmidt values
print('\n',"Sum of squares of schmidt Values: ",sumofschm2(S),'\n'); //Returns 1
//Manually computing the norm of psi
print("The norm of the wavefunction is: ",my_norm(phi),'\n'); //Returns 1
//U*S*V is equal to G
print("Norm of U*S*V - G: ",norm(U*S*V - G),'\n');
}
return 0;
}
Additionally, if I manually set the number of Schmidt values to n (by setting MaxDim and MinDim in svd), |\ket{\phi}|^2 ends up being equal to n.
P.S. Is there an easier way to calculate \text{tr}(S^2) ? Attempting to contract S*dag(S) throws an error