# Truncation error of the svd with maxdim on MPS

**URL:** <https://itensor.discourse.group/t/truncation-error-of-the-svd-with-maxdim-on-mps/1202>\
**Category:** ITensor Julia Questions\
**Created:** [October 9, 2023, 3:34pm UTC](https://itensor.discourse.group/t/truncation-error-of-the-svd-with-maxdim-on-mps/1202 "2023-10-09T15:34:45Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![scichy](https://avatars.discourse-cdn.com/v4/letter/s/f4b2a3/32.png) [@scichy](https://itensor.discourse.group/u/scichy)\
**Post date:** [October 9, 2023, 3:34pm UTC](https://itensor.discourse.group/t/truncation-error-of-the-svd-with-maxdim-on-mps/1202/1 "2023-10-09T15:34:45Z")

</div>

Hi!

I am sorry if this is a duplicate, as I would expect it to be answered somewhere but I could genuinely not find the answer. Shortly stated: I want to retrieve the sum of discarded squared singular values when using the svd on a MPS with the maxdim argument. Is there a way of doing so?

I can compute the svd without truncation, retrieve the singular values from S and then perform some truncation, but that is probably far from efficient. Digging in the source code I have also figured out that the svd function returns more than just U, S and V, and one can obtain a Spectrum object, which contains the (conserved) singular values and some error term: what is the latter? I could not relate it to the discarded singular values and could not find any documentation on that.

Here an example: I start from some MPS, and apply a unitary to two neighbouring sites

```auto
inds_sites = siteinds("Qubit", 8)
psi = randomMPS(inds_sites, linkdims=8)
q=4
new_T = psi[q] * psi[q+1] * op("RandomUnitary", [inds_sites[q], inds_sites[q+1]])
noprime!(new_T)

```

Then I want to perform a svd with maximum bond dimension. For reference, if I give no truncation indications I get

```auto
U, S, V, spec = svd(new_T, uniqueinds(psi[q], psi[q+1])); spec
>>> Spectrum{Vector{Float64}, Float64}([1.3888894384006771, 1.16280226828166, 1.0076611812149532, 0.9505778312199186, 0.7606731410439445, 0.6710854849239525, 0.595632286840479, 0.46336037208905767, 0.32191293790266035, 0.20206667411714405, 0.1809923124921696, 0.1336780702863157, 0.08791572340819503, 0.055008793515782126, 0.01707487116659689, 0.0006686130965095594], 0.0)

```

I know that the first element of spec is the list of the singular values (which are also on the diagonal of S) squared, so `spec.eigs == diag(S).^2`. Now when imposing a truncation I get

```auto
U, S, V, spec = svd(new_T, uniqueinds(psi[q], psi[q+1]), maxdim=8); spec
Spectrum{Vector{Float64}, Float64}([1.3888894384006771, 1.16280226828166, 1.0076611812149532, 0.9505778312199186, 0.7606731410439445, 0.6710854849239525, 0.595632286840479, 0.46336037208905767], 0.12491474949817141)

```

Now, is there a way to retrieve the same U, S and V, and also to get the sum of discarded singular values squared? In an attempt, I tried to figure out what this `spec.truncerr = 0.12491474949817141` is. I thought it would be truncerr = \sum\_{i=9}^{16} s\_i^2 or truncerr = \sqrt{ 2 \sum\_{i=9}^{16} s\_i^2 } but it seems not to be the case

```auto
discarded_singvals = [0.32191293790266035, 0.20206667411714405, 0.1809923124921696, 0.1336780702863157, 0.08791572340819503, 0.055008793515782126, 0.01707487116659689, 0.0006686130965095594]
sum(discarded_singvals) = 0.9993179959853734
sqrt(2*sum(discarded_singvals)) = 1.4137312304574539

```

Sorry for the long post, I tried to make my doubt as clear as possible.  
Some details about my settings:

- Julia Version 1.8.5
- ITensors v0.3.30

Thanks a lot for any help!

---

<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:** [October 9, 2023, 8:59pm UTC](https://itensor.discourse.group/t/truncation-error-of-the-svd-with-maxdim-on-mps/1202/2 "2023-10-09T20:59:15Z")

</div>

Thanks for the well-considered question and for trying out different test calculations to get a better idea of what’s going on. I think our documentation could be more detailed about both what the `svd` function returns and how each quantity, such as the truncation error, is computed.

But to answer your questions,

1. it’s important to know that the values stored in the `Spectrum` object are the _density matrix eigenvalues_, so the squares of the singular values. I mention that partly because I see you named them `discarded_singvals` but actually those numbers are the discarded squares of singular values. Let’s call the returned numbers the “eigenvalues” for short. You can call `eigs(spec)` to retrieve them as a Vector from the `Spectrum` object.

2. calling `truncerror(spec)` on the `Spectrum` object returned by the `svd` function returns a number which is given by the following formula:

3. One important detail about your code is that the MPS you are using has its “orthogonality center” at site 1, which is the behavior of our `randomMPS` function. So if you then act a gate on sites q and q+1 and SVD, the density matrix spectrum (squares of singular values) won’t add up to 1.0 as it would when acting on the orthogonality center of a normalized wavefunction. This choice can also make the truncation step “uncontrolled”. So I would recommend adding a line `orthogonalize!(psi,q)` to move the orthogonality center to site q before running the rest of the code you wrote.

Here is an expanded sample code for you, based on yours:

```julia
using ITensors
using Random

let
  Random.seed!(1)
  inds_sites = siteinds("Qubit", 8)
  psi = randomMPS(inds_sites, linkdims=8)

  # Apply gate
  q=4
  orthogonalize!(psi,q)
  new_T = psi[q] * psi[q+1] * op("RandomUnitary", [inds_sites[q], inds_sites[q+1]])
  noprime!(new_T)

  U, S, V, xspec = svd(new_T, uniqueinds(psi[q], psi[q+1]))

  eig_scale = sum(eigs(xspec))
  @show eig_scale

  maxdim = 8
  U, S, V, tspec = svd(new_T, uniqueinds(psi[q], psi[q+1]); maxdim)

  terr = 0.0
  for j=(maxdim+1):length(eigs(xspec))
    terr += eigs(xspec)[j]
  end

  @show terr/eig_scale
  @show truncerror(tspec)

  return
end

```

---

<div class="post-metadata">

**Author:** ![scichy](https://avatars.discourse-cdn.com/v4/letter/s/f4b2a3/32.png) [@scichy](https://itensor.discourse.group/u/scichy)\
**Post date:** [October 17, 2023, 10:54am UTC](https://itensor.discourse.group/t/truncation-error-of-the-svd-with-maxdim-on-mps/1202/3 "2023-10-17T10:54:36Z")

</div>

Thank you a lot! That’s a very helpful answer 🙂 I guess I just did not think about normalization  
First of all, I will definitely make sure of adding the `orthogonalize!` method. I have it in most of my implementations but forgot in this example.  
Then, I am interested in `terr`, not in `truncerror(tspec)`. To avoid computing two SVDs, would you say that computing

```auto
terr = truncerror(tspec) * LinearAlgebra.norm(psi)^2

```

is an efficient way of doing so?

I used the fact that the sum of all eigenvalues (when have the orthogonality center placed properly) should be the squared norm of the state

\sum\_{j=1}^\chi p\_j = \sum\_{j=1}^\chi s\_j^2 = \| |\psi\rangle \|\_2^2 \leq 1