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

from an InitState defined on a BosonSite - Define
`G = psi(1)*psi(2)`

, and perform its SVD (I have checked that`U*S*V == G`

) - Create some new MPS defined on the same BosonSite, called
`phi`

- Assign
`U`

to`phi(1)`

and`S*V`

to`phi(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