# Confusing behavior of SVD

**URL:** <https://itensor.discourse.group/t/confusing-behavior-of-svd/120>\
**Category:** ITensor C++ Questions\
**Created:** [May 5, 2022, 6:19pm UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120 "2022-05-05T18:19:13Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![vik](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vik/32/50_2.png) [@vik](https://itensor.discourse.group/u/vik)\
**Post date:** [May 5, 2022, 6:19pm UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/1 "2022-05-05T18:19:13Z")

</div>

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!

```auto
#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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [May 23, 2022, 1:57pm UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/2 "2022-05-23T13:57:39Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![vik](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vik/32/50_2.png) [@vik](https://itensor.discourse.group/u/vik)\
**Post date:** [June 1, 2022, 11:26pm UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/3 "2022-06-01T23:26:32Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [June 2, 2022, 4:53pm UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/4 "2022-06-02T16:53:49Z")

</div>

Glad you had figured out the main issue already!

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.

---

<div class="post-metadata">

**Author:** ![vik](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vik/32/50_2.png) [@vik](https://itensor.discourse.group/u/vik)\
**Post date:** [September 10, 2022, 12:25am UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/5 "2022-09-10T00:25:43Z")

</div>

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)

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

orthogonalize!(psi,3)

U,S,V = svd(psi[3], (siteind(psi,3),linkind(psi,2)))

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

#The following fail
S*dag(S)
S*prime(dag(S);tags="Link,v")

```

The error I receive is  
`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.

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [September 11, 2022, 3:35am UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/6 "2022-09-11T03:35:33Z")

</div>

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:

```auto
S = dense(S)

```

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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [September 11, 2022, 3:36am UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/7 "2022-09-11T03:36:38Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![vik](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vik/32/50_2.png) [@vik](https://itensor.discourse.group/u/vik)\
**Post date:** [September 12, 2022, 11:07pm UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/8 "2022-09-12T23:07:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [September 13, 2022, 7:29am UTC](https://itensor.discourse.group/t/confusing-behavior-of-svd/120/9 "2022-09-13T07:29:44Z")

</div>

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.