[ENHANCEMENT] Expose `truncerr` in calls to `truncate!`

[NuclearPowerNerd]

Hello,

I'd be interested on thoughts here. I think it would be useful to be able to get the truncation error corresponding to a call to truncate!.

Currently truncate! as defined in abstractmps.jl only returns the MPS, but each truncation that is performed is a call to svd which returns the TruncSVD type which contains the truncation error.

Maybe one solution then would be to allocate a small vector with as many elements as there are bonds of the MPS. Then in the right to left sweep just collect the truncation error into each element and have that be optionally returned (could make it a kwarg).

If this functionality is already available then please just disregard and close this feature request.

Here is some psuedocode for what I had in mind (I've not actually tested this to see if it works, just what I was thinking). If this is something that sounds reasonable, I can prepare the pull request. I can make the corresponding updates to truncate as well.

function truncate!(
  ::Algorithm"frobenius", M::AbstractMPS; site_range=1:length(M), return_truncation_error=false, kwargs...
)
  N = length(M)
  truncerr = zeros(N-1)

  # Left-orthogonalize all tensors to make
  # truncations controlled
  orthogonalize!(M, last(site_range))

  # Perform truncations in a right-to-left sweep
  for (idx, j) in enumerate(reverse((first(site_range) + 1):last(site_range)))
    rinds = uniqueinds(M[j], M[j - 1])
    ltags = tags(commonind(M[j], M[j - 1]))
    U, S, V, spec = svd(M[j], rinds; lefttags=ltags, kwargs...)
    truncerr[idx] = spec.truncerr
    M[j] = U
    M[j - 1] *= (S * V)
    setrightlim!(M, j)
  end
  if return_truncation_error
  return M, truncerr 
  else
   return M
   end
end

[emstoudenmire]

I think this is a good feature for us to pursue, in the sense of making this information available through a return value. Probably, though, we should do it by having truncate return the MPS and a NamedTuple with various things in it, including these truncation errors. The nice thing about a NamedTuple is it would let us add other things in the future and also users could pass keyword arguments to opt in to (possibly) more expensive or technical things that could be returned but which we wouldn't want to return by default.

[NuclearPowerNerd]

Hey Miles - sounds good. When I originally posted this I also was thinking of truncate! in isolation. But I can see how it would be useful if you could return this on any call that eventually ends up calling truncate, such as contract or add (when using densitymatrix).

For the named tuple approach, were you thinking of something like this? If not, do you have any examples within the current ITensor code I could take a look at for reference?

function truncate!(
  ::Algorithm"frobenius", M::AbstractMPS; site_range=1:length(M), return_spec=false, kwargs...
)
  N = length(M)
  spec = (;truncerr=zeros(N-1)) # calling it spec to keep it generic? Since other things might be added in the future

  # Left-orthogonalize all tensors to make
  # truncations controlled
  orthogonalize!(M, last(site_range))

  # Perform truncations in a right-to-left sweep
  for (idx, j) in enumerate(reverse((first(site_range) + 1):last(site_range)))
    rinds = uniqueinds(M[j], M[j - 1])
    ltags = tags(commonind(M[j], M[j - 1]))
    U, S, V, spec = svd(M[j], rinds; lefttags=ltags, kwargs...)
    spec.truncerr[idx] = spec.truncerr
    M[j] = U
    M[j - 1] *= (S * V)
    setrightlim!(M, j)
  end
  if return_spec
  return M, spec
  else
   return M
   end
end

[NuclearPowerNerd]

Hello! Has there been any additional thoughts on this, or the interface it should follow?

I just came across a situation again where it would have been nice to somehow recover the truncation error computed by TruncSVD.

I'm happy to develop a PR if we could agree on what this should look like.

[mtfishman]

What I would have in mind is:

function truncate!(
  ::Algorithm"frobenius", M::AbstractMPS; site_range=1:length(M), (truncation_error!)=nothing, kwargs...
)
  # Code
  if !isnothing(truncation_error!)
    truncation_error![] = truncerr
  end
  return M
end

So if a user wants the truncation error they can call it like this:

truncation_error_ref = Ref{Vector{Float64}}()
truncate!(M; (truncation_error!)= truncation_error_ref)
truncation_error = truncation_error_ref[]

Basically, the design is that you pass a reference that is able to store the truncation errors that gets written into so the user can access the truncation errors.