Add support for LAPACK SVD routines - #333
Conversation
|
Thanks for a great looking PR. This is something we needed, I agree. It's something I had been planning to add for a while, but other things just got in the way. I'm not too surprised about the difference in accuracy of the methods to our home-made one: I recently implemented the same algorithm in Julia and found it gets a lot more accurate if we substitute a QR factorization for the "orthog" step, which was always a bit suspicious. What I hadn't quite known how to do until recently was use QR properly here, but I figured out that they key was to make the diagonal of R non-negative to ensure uniqueness of the result. Then it gives a way to recover the singular values (and V) more accurately. I'd be curious how the accuracy then compares to the LAPACK methods. That being said, I'm certain the LAPACK methods are faster. The reason we didn't just use them from the very beginning is that we've found many times that they occasionally crash. So we've had our home made one in place because it's essentially never had any crashes. So that's why I agree with your choice of keeping the default as ITensor; we can recommend to users who find the SVD to be a bottleneck to try the other two methods as a possible optimization. Happy to discuss more - |
|
Great, OK, glad the PR was helpful and interesting to hear that it was the stability issue in the LAPACK versions which were the issue. I guess now QR is in ITensor it shouldn't be too hard to update the current default implementation to your Julia one. [There's also a method using Jacobi preconditioning for extreme accuracy in the s.v.s in LAPACK but my guess was this was beyond the needs of typical ITensor users.] Just to make sure, did you see the part at the end about the thresh argument being removed, and the value for the default? |
|
Hi, thanks for emphasizing that. I did see it but I agree it's worth thinking over the impact to the user. Overall, I'd say that replacing it with an Args argument is a better design anyway, and we can still allow it as an optional arg. So it's a slight breaking change, but I bet hardly anyone is using it themselves in practice. Appreciate it! I'll try to get around to improving the ITensor home-made SVD soon with that QR idea (which we now have in ITensor thanks also to your previous PR), and then it would be interesting to compare the accuracy again. I'm not expecting it to beat LAPACK, but I'm hopeful it will be comparable. |
|
Its nice to see this functionality got integrated into ITensor in the end !
|
Add support for LAPACK SVD routines
3 participants
Adds an optional argument to svd "SVDMethod", allowing the selection of either the ITensor implementation "ITensor", the divide-and-conquer LAPACK method "gesdd" or the QR LAPACK method "gesvd". Effectively #179 updated to ITensor v3, with gesvd, and slightly reorganised to avoid code duplication when calling multiple different SVD routines, so full credit to @jurajHasik for the majority of the code (I only didn't build on top because pulling their branch onto v3 merged in some unrelated v2 stuff as well).
Use Cases: Benchmarking, both of the ITensor routine and against other MPS/DMRG codes which often implicitly use gesdd through their languages' SVDs (e.g. numpy, Julia). [This is why I personally needed this.] Also, both gesdd and gesvd seem much faster and more accurate (see attached graphs based on 1000 runs of the "Accuracy Stress Test" unit test in ITensor on a matrix with rapidly decreasing s.v.s with a single thread (MKL)), so unless the matrices are pathological it would seem to make sense to use gesdd, or gesvd if that fails (see discussion about gesvd vs gesdd e.g. here: https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603)
Notes: Currently would break any user code which called the thresh argument explicitly, but could easily be fixed with an overload to not do so if this is a problem -- I'm not sure if that argument is user-facing anyway? On the subject, the default in the old code is currently SVD_THRESH = 1e-5 for matrices and 1e-3 for ITensors (from auto thresh = args.getReal("SVDThreshold",1E-3)) -- which is more sensible?