# Confusing behavior of SVD

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:

1. Create an MPS psi from an InitState defined on a BosonSite
2. Define G = psi(1)*psi(2), and perform its SVD (I have checked that U*S*V == G)
3. Create some new MPS defined on the same BosonSite, called phi
4. 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

Hi, sorry for the slow reply here. This post got caught up in the forum spam filter and I’ve just adjusted that filter to hopefully not flag as many posts in the future.

Thanks for the question and great to see you are getting your hands on some of the lower level parts of ITensor and testing out various functions.

Do you still have this question? One comment is that some of the things you are wanting to do, such as compute \Tr[S^2] would be easier to do in the Julia version of ITensor since the interface there is a bit more flexible than the C++ version. Have you considered using the Julia version?

But also I found the issue with your code above. It’s that your implementation of my_norm has a bug in it: when you are contracting the tensors of the MPS and its conjugate together, they will have the same link or bond (or virtual) indices, so in the line auto C = dag(wf(1))*wf(1); those link indices are being contracted together, which is not the behavior you want. What you need to do there is to call prime(dag(wf(1)),"Link") to prime any indices which have the tag “Link” to prevent that contraction from happening, and only contract over the site indices there.

Please let me know if you have any questions about that - it’s helpful when writing code like that to print out each intermediate tensor to verify it has the indices you expect and also to diagram on paper what is happening at each step.

Best regards,
Miles

Hi Miles,

Thanks so much. I had replied to the Automod message asking that the question not be reinstated (when I first received the flag notification) since I had figured out that each contraction was giving me an order 1 tensor because of the lack of priming, and assumed that this message had been seen because I couldn’t find the post after that.

I have not considered using the Julia version, and didn’t think I would need it thus far, beyond computing \text{tr}(S^2) or \text{tr}(S^n) a little more “diagrammatically”. Printing out the actual matrices S and dag(S) would have me believe that their indices should contract, but am I right in inferring that this does not apply generally to QDiagReal objects?

Regarding QDiagReal, those are still ITensors but just have a different storage type internally. So they still follow exactly the same contraction rules as for other ITensors.

1 Like

Hi Miles,

Sorry to come back to this again, and after so long. I took your advice and shifted to using Julia, and while this isn’t a prohibitively concerning issue, I’m still running into trouble while calculating tr(S^n) “diagrammatically”, when using QNs. This is the MWE (in Julia)

sites = siteinds("Boson",4;dim=3,conserve_qns=true)
states = [1 2 1 1]

orthogonalize!(psi,3)

U*dag(U)    #Works okay
V*dag(V)     #Works okay

#The following fail
S*dag(S)


MethodError: zeros(::Type{NDTensors.DiagBlockSparseTensor{Float64, 0, Tuple{}, NDTensors.DiagBlockSparse{Float64, Vector{Float64}, 2}}}, ::Tuple{}) is ambiguous.

I am not sure where zeros is used in the code; I’m assuming in initializing matrices in the method in LinearAlgebra or NDTensors, since it isn’t there in decomp.jl (which just calls one of these svd functions for the actual decomposition).

I completely understand if this is not worth the team’s time to sort out or fix, and if there isn’t a solution (or an error in what I’m doing), I would be more than happy to dig a little deeper. Also, please let me know if I should repost this in the Julia section of the discourse group. Thanks a lot!

P.S. I am aware that the actual values that one needs for calculating entropy are the entries of S squared; the n in tr(S^n) isn’t meant to be the Renyi index.

I think the errors you are getting are bugs or incomplete features in ITensor, unfortunately. This because S has a special type of storage that is “diagonal sparse” storage which we use for efficiency. But often that storage has a more limited set of features it supports - basically we would like to fix this so I’ll file a bug report about it.

The good news is there are some things that you can do as a workaround here. The simplest is just to do:

S = dense(S)


then write one of the lines of code you were trying, like S*dag(S) which should now work.

Actually I am finding that S*dag(S) works ok even when S is the diagonal sparse ITensor obtained from the SVD. So I can’t reproduce that bug.

What version of ITensors.jl are you using?

I was using v0.3.19, but upgraded to v0.3.20 today. The problem still persists, however.

I see, thanks. Could you then share a minimal working code that reproduces the error you are getting? I was not able to reproduce it myself.