Svd truncate with minblockdim

[manuschneider]

In Svd_truncate and Gevd_truncate, I added an optional argument minblockdim. With this, one can set a minimal bond dimension for each block.
Additionally, some cleanup and problems were fixed in the original implementation with mindim. The standard was set to mindim=1.

[codecov]

Codecov Report

Attention: Patch coverage is 4.89297% with 311 lines in your changes missing coverage. Please review.

Project coverage is 16.69%. Comparing base (72453ca) to head (7c62492).
Report is 31 commits behind head on dev-master.

Files with missing lines	Patch %	Lines
src/linalg/Gesvd_truncate.cpp	3.65%	154 Missing and 4 partials ⚠️
src/linalg/Svd_truncate.cpp	6.13%	142 Missing and 11 partials ⚠️

Additional details and impacted files

@@              Coverage Diff               @@
##           dev-master     #491      +/-   ##
==============================================
+ Coverage       16.60%   16.69%   +0.09%     
==============================================
  Files             221      221              
  Lines           48490    48302     -188     
  Branches        20272    19900     -372     
==============================================
+ Hits             8053     8066      +13     
+ Misses          36143    35973     -170     
+ Partials         4294     4263      -31     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[IvanaGyro]

What is the difference between keepdim and mindim? Can we leave only one of them? There should be only two use cases:

1. truncated with the low error bound
2. truncated with the minimal dimensions

Please correct me if I am wrong~

[manuschneider]

I added a detailed explanation of the truncation order in the API documentation:
The truncation order is as following (later constraints might be violated by previous ones):

1. Keep the largest min_blockdim singular values in each block; reduce keepdim and mindim by the number of already kept singular values
2. Keep at most keepdim singular values; there might be an exception in case of exact degeneracies where more singular values are kept
3. Keep at least mindim singular values
4. Drop all singular values smaller than err (no normalization applied to the singular values)

Concerning the input parameters: currently, input arguments are taken as const and, if they need to be changed in a program, will be copied. This ensures two things:

1. The inputs will not be changed, and thus no unexpected behavior outside the function
2. The user can use a simple API and does not need to deal with constant pointers etc

This creates some overhead, but I think it is negligible in the cases here.

[kaihsin]

I think the meaning of mindim (specifically for with symmetrty) is different than the the keepdim. one is upper bound and one is lower bound (with the truncation rate). @yingjerkao could explain it more on this topic. For pybind issue, I am happy to chat on this, just shoot me a msg.

[IvanaGyro]

Concerning the input parameters: currently, input arguments are taken as const and, if they need to be changed in a program, will be copied. This ensures two things:

1. The inputs will not be changed, and thus no unexpected behavior outside the function
2. The user can use a simple API and does not need to deal with constant pointers etc

This creates some overhead, but I think it is negligible in the cases here.

The API document is clearer after the change!!

When a container is passed as a const-value or a non-const-value, the container is copied. Modifying the value argument inside the function will not cause effects outside the function. This means just passing the argument as a non-const-value if we decide to pass the argument as a value instead of a reference. See this StackOverflow discussion and this demo.

And the reason I suggest expanding minblockdim in the upper function is that we can reduce one or two copies. Here is the demo.

[ianmccul]

Two brief comments:
I generally agree with the comments above on pass by value vs pass by reference. To be clear, in modern C++, if you are going to modify the object locally and you don't want those changes to appear as a side-effect, then pass by value is fine. In many cases the compiler will be able to avoid a copy anyway, and 'move' the parameter. Because many of the LAPACK functions destroy the input object, the gesvd function itself is a good candidate for taking the input matrix by value; if you want to keep the original matrix around then you need to copy it anyway, and if you don't need the input matrix anymore then you can 'move' it and avoid the copy. So it could be

std::vector<Tensor> Svd_truncate(Tensor Tin,  .... // pass by value here

Currently, Svd_truncate() ends up calling Sdd_internal_xxx(), which always makes a copy of the input. This isn't necessary; you could pass the Tensor by value and leave it up to the caller as to whether they want to make the copy.

   Tensor X = ....;
   ...
   Svd_truncate(X);   // will copy X; it needs to do this anyway since Lapack overwrites the storage
   // X is unmodified here

   // ...

   Svd_truncate(std::move(X), ...);   // if we don't need X anymore
   // X is destroyed here

   // Calling Svd_truncate via a temporary (eg the return value of another function call) will automatically move rather than copy
   Svd_truncate(ConstructTensor(...), ...);    // assuming some function Tensor ConstructTensor(...).  No extra copy

Currently, in those last two cases Cytnx will make a copy of the tensor, and then immediately destroy it.

Secondly, in practice you might want to do something moderately complicated with the selection of singular values to keep, especially in the UniTensor case where there are many quantum number sectors. In my Toolkit, depending on the algorithm, there is an overall minimum number of states, minimum number of states per quantum number sector, a flag to indicate whether this minimum number of states per quantum number sector should include states with zero singular value or not (depending on the type of bond dimension expansion, both choices are possible), maximum number of states, a cutoff, a truncation error, a flag to indicate whether the truncation error is absolute or relative. All of these options are actually used in various different situations. I implement this in two stages; there is a class object that does the actual SVD and constructs a list of the singular values, and it appears as a (sorted) container of singular values. The caller can then select which singular values to keep via a helper function that iterates through the list. In simple cases this is just a subset of the original list (eg if you want to keep the N largest singular values), but in more complicated cases (eg you want to keep a minimum number per quantum sector) the states you want to keep are not simply the largest N so you need to assemble a list or vector of the states that you want to keep.

[manuschneider]

@IvanaGyro You are right, the vector is copied already anyways when passed to the function, so we can make it non-const and change it internally as in your code example SvdExpandInside, only that the dimensions are determined in the function and don't need to be given as an argument. I wanted to avoid the other way, SvdExpandOutside. In this case, the user can only pass a vector of the correct length. In many (most) cases, one wants to have the same min_blockdim for all sectors though. In this case the user does not have to think about the number of sectors in S after the SVD (the number is not equal to the number of sectors in the input tensor!). So I want to do this expansion inside the SVD_truncate function.
I changed the code accordingly, please check.

@ianmccul I agree, there are more ways to truncate the SVD and since this is the crucial step in many algorithms many options might be important (and I missed a relative truncation error as well). In addition to your suggestions, another topic is how to deal with degeneracies. Itensor provides an option for this as well.
Currently, if the singular values around the cutoff are exactly degenerate, then the bond dimension of the resulting tensor can be larger than keepdim.
The current implementation works as follows:

1. keep the min_blockdim singular values in each block and add the others into a vector Sall
2. Find the smallest singular value that is still kept according to the other criteria
3. Truncate all singular values in each block that are smaller than this value

If one wants to have keepdim or do it similar to your function, @ianmccul , I think the function has to be organized differently. For example, the singular values can be ordered and the permutation saved, such that one can relate each value to a sector. This would mean major changes, since our sort function does not return the permutation.

In general, I think we can go in two ways:

1. Implement all meaningful truncation criteria
2. Provide an easy way for the user to implement custom truncation functions

[ianmccul]

@manuschneider, depending on the context, I have an overall minimum number of states (eg in DMRG), or a minimum number of states per sector (which the bond expansion code uses, to make sure that all viable quantum number sectors are included. Usually this number is either 1 or 0, but you could set it to > 1 if you want). Other than that, I've never been concerned about the actual number of states in each sector -- I can't imagine the circumstances where I'd want to construct an array for the number per sector (except as some implementation detail). But maybe this has some uses in special cases?

I'd forgotten about degeneracies - I've never worried about this in my own code since it usually means there is some non-abelian symmetry that you should be using (so just use it!), and for large bond dimension it doesn't have much effect. But yes it is better to keep degenerate blocks, and I might add that to my code sometime. Just need a good algorithm for finding degeneracies. It looks like itensor just checks to see if the next singular value is within 1E-3 (or probably eigenvalue [squared singular value] rather than singular value itself). That seems like a bit of a random choice. This is https://github.com/ITensor/ITensor/blob/ea88f512d258ef23b5554950639c5acedf4df1f0/itensor/decomp.cc#L437 I cannot see anything corresponding in the Julia version of itensor, it looks like they got rid of it?

But something that does arise with non-abelian symmetries is how to count the weight of a state. Do you order an N-degenerate multiplet, each of weight w, by w itself or N*w? This makes a difference for which states you keep at the bottom of the spectrum, eg do you keep some degeneracy 1 state with weight epsilon, or do you instead keep a degeneracy 3 state with weight epsilon/2 ? Minimizing the truncation error corresponds to including the degeneracy in the weight, but that gives a different selection of states than if you use w itself (and therefore a different selection versus if you didn't use the symmetry). in practice, for finite DMRG, I include the degeneracy but in iDMRG I do not (since it makes it slightly unstable, to favor higher spin representations you can end up with the calculation going haywire). So this needs to be another option when doing the truncation!

Adding some functions to separate the SVD itself from the selection of states to keep doesn't mean throwing away the existing Svd_truncate() functions. Maybe you want to keep them for simple cases anyway, and they could be implemented in terms of a more general 2-part SVD followed by truncation. iTensor already has a large number of possible arguments for their truncate, and even then it doesn't cover all of the possibilities that I actually use.

[manuschneider]

The use case I had in mind was the following:
If we have a quantum number like the magnetization, and we know that the state of interest has a low total magnetization. Then I could force the blocks with low quantum numbers (say, -1, 0, 1) to have at least one singular value kept in order not to loose an important sector. But a quantum number of 10 could be dropped if the singular values become small because I do not expect it to be relevant. Allowing for a vector gives the user more flexibility and it did not make much difference to implement it. But not sure if one would really use this, otherwise just one number as an argument would be fine as well.

I agree that it would be best to have an implementation of Svd_truncate for common cases, and allow for customized truncations with not so much effort. For this, it would be nice to move the truncation step itself to a method of UniTensor that the user can access. And to have a version of Svd that sorts the singular values (possibly with further information like which block does each singular value correspond to).

[IvanaGyro]

I didn't go deep into the algorithm. Except for the algorithm, LGTM.

[IvanaGyro]

Does this line remove the largest min_blockdim[b] numbers?

[IvanaGyro]

@IvanaGyro You are right, the vector is copied already anyways when passed to the function, so we can make it non-const and change it internally as in your code example SvdExpandInside, only that the dimensions are determined in the function and don't need to be given as an argument. I wanted to avoid the other way, SvdExpandOutside. In this case, the user can only pass a vector of the correct length. In many (most) cases, one wants to have the same min_blockdim for all sectors though. In this case the user does not have to think about the number of sectors in S after the SVD (the number is not equal to the number of sectors in the input tensor!). So I want to do this expansion inside the SVD_truncate function.
I changed the code accordingly, please check.

@manuschneider The if-clause is in _svd_truncate_Block_UT() not Svd_truncate() now. If we move the if-clause branch into Svd_truncate(), both functions can consume const-references and avoid copy when the size of min_blockdim equals the number of blocks. Under this change, users can still pass only one min_blockdim value. This is just a suggestion; it's up to you whether to adopt it. If the number of blocks is usually not large, it shouldn't have a significant impact.

What will happen if I dismiss this review?

[ianmccul]

@manuschneider selecting states with largest singular values will automatically select the important sectors, it shouldn't normally be necessary to force it. And if something has gone wrong (eg some sector that you think should be important has zero singular value) then forcing it into the basis probably won't work anyway.

My code for truncations is at https://github.com/mptoolkit/mptoolkit/blob/main/mps/truncation.h and - https://github.com/mptoolkit/mptoolkit/blob/main/mps/density.h I wouldn't necessarily copy this style in new code, it grew somewhat organically and isn't exactly what I'd do if I was writing it from scratch. But the basic approach, of an SVD that gives a container of eigenvalues (the EigenInfo values in truncation.h) is quite flexible. In simple cases we just choose the first N elements according to the criteria in the StatesInfo struct (MinStates, MaxStates, TruncationCutoff, EigenvalueCutoff, and a flag to indicate whether the truncation is relative or absolute). In more complicated cases where you want to keep a selection of states that is not in the same order (eg because you want to keep some minimum number of states in each quantum number sector, even if there are larger singular values that you throw away) then you need to assemble your own list of states to keep, eg the TruncateExtraStates() function is used by the bond expansion code, and returns an std::list<EigenInfo> rather than an end iterator.

I also have (at least) 3 variants of the SVD, depending on what you want to do with structurally zero singular values. i.e. given an $n\times m$ matrix, do you want the dimension of the diagonal D matrix of singular values to be size $m$, or size $n$, or $\min(m,n)$? If you are forcing states into the basis with a minimum number per quantum number sector (or the truncation settings allow for strictly zero singular values) then you want it to be $m$ or $n$ (depending on which direction you are moving), and include singular values that are strictly zero. I do this in the bond expansion code, where having a zero singular value of the expansion tensor doesn't mean that the state will necessarily have zero weight in the wavefunction! (The CBE paper also does this, in fact I think it is even more important in their scheme.)

@manuschneider , @IvanaGyro , I can't imagine any circumstances where the cost of copying the vector of block dimensions is significant; the amount of work is insignificant compared with doing the SVD and performing the truncation. Much more significant is avoiding making a extra copy in the case where you don't need the tensor anymore (and LAPACK dgesvd can safely overwrite it, or even better, re-use that memory directly into the U/V matrix). My QR decomposition code in the Toolkit does exactly this - the tensor is passed by value and the returned Q matrix uses the same memory, so if you do a QR decomposition of a temporary or an std::move'd object then it doesn't need to allocate new memory for Q. I haven't updated my SVD functions to do the same yet, but in principle it is fairly straight forward. But I guess in Python there is no way of controlling this anyway? I guess you could have some in-place variant? Or a variant that destroys the original tensor? that might be a bit obscure.

[ianmccul]

@manuschneider if you want to avoid the copy of the min_blockdim array, then this line¹ needs to be

-        _svd_truncate_Block_UT(outCyT, Tin, keepdim, min_blockdim, err, is_UvT, return_err, mindim);
+        _svd_truncate_Block_UT(outCyT, Tin, keepdim, std::move(min_blockdim), err, is_UvT, return_err, mindim);

Footnotes

1. https://github.com/Cytnx-dev/Cytnx/blob/7c62492a3d80aeb5c871c0d4e716749ebf040b6d/src/linalg/Svd_truncate.cpp#L585C9-L585C31 ↩