Return truncation error and singular values from replaceBond and factorize

[orialb]

Currently replaceBond, factorize, svd and eigen don't return the truncation error which we would like to track in many cases. It is also necessary in some cases to have the actual singular values from replaceBond and factorize.
It would be simple to just propagate this through all these functions, but might make the code a little more cumbersome in cases when one doesn't need this information (although not much more).
Another option is to have a factorization object returned from factorize, svd and eigen functions which will contain all the information.
A third option is to return this information from factorize, svd and eigen only when it is asked for by the caller with a keyword argument.

Any thoughts?

[mtfishman]

It seems like a good idea. We already pass the newly introduced indices from those functions, which you can grab or not just by including them in your output list:

U,S,V = svd(A,(i,j))
U,S,V,u,v = svd(A,(i,j))

We can just add the truncation error as another output after the indices, it wouldn't affect users who don't want it. So people could get it like this:

U,S,V,u,v,truncerr = svd(A,(i,j))

The list could get a bit long, but I think that is ok (and we can always pack things after U,S,V into an object if we end up outputting a lot of stuff).

[emstoudenmire]

I second the idea of putting things into an object beyond a certain point. A good question would be whether the u and v inds should be included in that object too.

In C++ ITensor there is a Spectrum object returned by the SVD which includes not only the truncerr, but all of the singular values (these are already in memory anyway, so no overhead incurred by returning a reference to them). It could also return the number of singular values kept and discarded (these can be recomputed but it's just convenient to have them at hand).

Using an object instead of a longer list of return values would make existing code more "future proof" since we could add things to the object without breaking anyone's code.

[orialb]

Yes, it would be good to return also the singular values so one can track entanglement entropy during sweeps.

A Spectrum object would be useful for that, since, depending on the internal factorization used, factorize would return either the singular values, the square of the singular values or the square root of the singular values. Having a dedicated object will allow to store the information about the type of the returned spectrum, so user could just call entropy(spec) (for example) without having to think whether they need to square/sqrt the spectrum before. But this might be also overcomplicating things.

Otherwise we could just transform the spectrum such as to always return singular values (or square of them if that's better), but I think this might incur some extra overhead depending on the internal factorization used.

[emstoudenmire]

Thanks, these are some good considerations.

In the C++ version, we just return a vector of the squared singular values (density matrix eigenvalues). In many cases, these are computed anyway because they are the input passed to the truncate subroutine called if a truncation of the SVD is performed. But yes this would be extra overhead if no truncation is performed. Maybe if the truncate=false keyword is passed, these won't be included in the returned object, unless another keyword argument is passed asking for them? For large matrices, the overhead is negligible either way of course.

One smallish issue with returning the squares is a loss of precision due to floating point behavior. But for most purposes, like computing entanglement entropy, this should not typically be a serious issue I don't think.

Either way, having this type defined would allow us to try this different designs behind the scenes without breaking any, or not much user-facing code in the future (if absolutely needed).

[orialb]

ok, I could implement something along these lines. It would be probably smart to wait that #109 is merged first because it seems like a big change that affects all the factorization functions (if I understand correctly).

[mtfishman]

Yes, the factorization functions are reworked quite a lot in #109. If this is not urgent, it would be easier to wait for that to land.

[orialb]

Closing this issue, as it was implemented by #139.