Truncation error of the svd with maxdim on MPS

ITensor Julia Questions
1

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

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

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

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

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!

2

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:

    \epsilon = \frac{\sum_{j=\chi'+1}^{\chi} p_j}{\sum_{j=1}^{\chi} p_j} = \frac{\sum_{j=\chi'+1}^{\chi} s^2_j}{\sum_{j=1}^{\chi} s^2_j}

    Above, \chi refers to the original (or maximal) rank of the tensor you’re SVD’ing, while \chi' refers to the truncated rank or bond dimension set by the maxdim argument in your case.

    So a key point is that the reported truncation error is divided by the sum of all of the untruncated eigenvalues. Below I’ll show you some code that demonstrates this.

  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:

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
3

Thank you a lot! That’s a very helpful answer :slight_smile: 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

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