When you use lanczos algorithm, do you employ krilovkit package or do you have one of your own?
P.S: I still need to try this library out.
I already used krilovkit and I stumbled upon convergence issues
Yes, we use KrylovKit.eigsolve in ITensorMPS.dmrg and also ITensorTDVP.dmrg.
I have a big chunk of code which leads to convergence issue and it’s probably an issue of mine.
(when there is no convergence issue (for smaller sizes or smaller bond sizes) the results are right. Putting an initial guess for the groundstate of the matrix helped pushing the convergence issues a bit further, but then they reappeared again)
But I noticed several strange issues with that program, I didn’t stumble upon when using mkl and arpack with fortran and C or eigen and spectra with C++.
Can I sum them up, with that annoyingly long code (explained), and see if they are errors I introduced by myself or not?
I can do this only from the first days of july
It is better to not share a very long code and instead try to reduce your issue down to a minimal example. As you can imagine, we are all quite busy (often debugging our own long complicated codes) so helping us out as much as possible is ideal. To be blunt, the shorter and more self-contained that the code you share is, the more likely it will be that we will read it and find the time to help you.
If you suspect it is an issue with KrylvKit.eigsolve, it is easy enough to insert checks into ITensorMPS.dmrg to see that it is outputting the correct eigenvector. For example, you can see the call to eigsolve here: ITensors.jl/src/lib/ITensorMPS/src/dmrg.jl at v0.6.14 · ITensor/ITensors.jl · GitHub and you could check how close the result is to being an eigenvector with a line like norm(PH(vecs[1]) - vals[1] * vecs[1]). I think that would be a good place for you to start, and that’s what I would do if I was seeing an issue in DMRG that I suspected was caused by the eigensolver. If KrylvKit.eigsolve isn’t outputting the correct eigenvector, then you can save the effective Hamiltonian and report it as a bug to KrylovKit.jl.
What you are describing sounds more like an issue with DMRG convergence, see DMRG FAQs · ITensors.jl for some hints about that. We’ve been using KrylvKit.eigsolve for quite a while and have found it to be very reliable, though of course it is possible there is some issue we haven’t come across yet.
The fact is that it actually didn’t converge at all.
One issue I definitely incurred in was the fact that multiplication of sparse matrix lead to quite full matrices with near zeros value, due to numerical errors. This didn’t happen this way in C or fortran (numerical errors where distributed among non zero values for some reasons). So I had to manually remove all the non zeros under a defined tolerance I put (found by reaching a saturation in the results), in order to both increase speed and reduce the amount of useless memory the matrix was storing.
This is not an issue of krilovkit by itself, rather the standard multiplication algorithm.
I don’t know if you encountered it. I found it when fully diagonalizing and it gave me wrong results. Then I tried to pick just the upper triangular of what was supposed to be a diagonal matrix, and it converged to quite the good results. Then removed zeros and bam, exact results.
After deploying this tolerance cut, adding a small magnetic field to choose which of the 2 ground states of the Hamiltonian I wanted, lead me to really wrong results. I had to couple the magnetic field only to the first site of the chain.
Then converge issues arose and where they appeared was based on which random seed I used. The results when they converged were really good.
As explained above I wrote the algorithm to use a ground eigenvectors guess, and it pushed the issue to bigger sizes. But it was still there
This is an example:
B^T A B = A
if B is unitary
I expect A to have some little noise, but have we see C++ distributes it in an entire different way than julia:
C++ program:
#include <iostream>
#include <iomanip>
#include "../../Eigen/eigen-3.4.0/Eigen/Dense"
using namespace Eigen;
using namespace std;
int main(){
double a = 1/sqrt(2);
Matrix4d A;
A << 0.5, 0, 0, 0,
0, 0.5, 0, 0,
0, 0, 0.5, 0,
0, 0, 0, 0.5;
Matrix4d B;
B << 0, 1, 0, 0,
a, 0, 0, a,
-a, 0, 0, a,
0, 0, 1, 0;
Matrix4d C = A * B;
Matrix4d D = B.transpose() * C;
// Print the result
cout << setprecision(19) << D << endl;
return 0;
}
compiled with:
g++ filename.cpp -o filename
output:
0.499999999999999889 0 0 0
0 0.5 0 0
0 0 0.5 0
0 0 0 0.499999999999999889
while on julia:
using LinearAlgebra
function myfun1()
a = 1/sqrt(2)
A = Float64[0 1 0 0; a 0 0 a; -a 0 0 a; 0 0 1 0]
B = Float64[0.5 0; 0 0.5]
C = Float64[1 0; 0 1]
A1 = kron(C, B)
display(A1)
A2 = A1 * A
print("\n")
display(A2)
display(A')
B = A' * A2
print("\n")
for i in 1:4
for j in 1:4
#pretty print
if i == 1 && j!=4 i==4 && j!=1
print(round.(B[i,j]; sigdigits=19))
else
print(round.(B[i,j]; sigdigits=2))
end
print(" ")
# pretty print
if i==2 && j==1 i ==3 && j==1
print(" ")
elseif i==4 && j==1
print(" ")
end
end
print("\n")
end
end
myfun1()
the output is:
4×4 Matrix{Float64}:
0.5 0.0 0.0 0.0
0.0 0.5 0.0 0.0
0.0 0.0 0.5 0.0
0.0 0.0 0.0 0.5
4×4 Matrix{Float64}:
0.0 0.5 0.0 0.0
0.353553 0.0 0.0 0.353553
-0.353553 0.0 0.0 0.353553
0.0 0.0 0.5 0.0
4×4 adjoint(::Matrix{Float64}) with eltype Float64:
0.0 0.707107 -0.707107 0.0
1.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0
0.0 0.707107 0.707107 0.0
0.4999999999999999 0.0 0.0 -1.1e-17
0.0 0.5 0.0 0.0
0.0 0.0 0.5 0.0
-1.1e-17 0.0 0.0 0.4999999999999999
I used full matrices, but in reality they should be sparse. So off diagonal values, in this case, are saved in memory.
This means that on one hand we are capable of removing a lot of numerical error, which is redirected on the non zeros values, which we know should be zeros. On the other hand this leads to incorrect results if left untreated, as in my case, and leads to increased computing time if treated, as we need to prune this values out of a sparse matrix, as they occupy most of the memory in the computer.
Thanks for the minimal example, that’s helpful for trying to understand the issue you are discussing. As you probably know, this appears to just be an artifact of how floating point arithmetic works.
If you are hoping for the results to be elementwise sparse, what about making the matrices sparse to begin with?
It was just the tip of the iceberg, not the issue of convergence.
I cannot provide an understandable minimal code of that, as I need 2 days to rewrite it in such a way.
This piece of code I wrote now created a lot of issue and was the first of the problems I encountered on Julia.
Matrices will be really big, as you know in DMRG. These little non zeros just go everywhere and pollute most of the sparse matrix with artifacts, leading to increased memory usage and increased time spent (multiplication of sparse matrices are faster cause they don’t multiply zeros).
At the end I just put the matrix in an exact diagonalizator and the results were wrong since the first digit (not just decimal).
Picking the upper triangular probably suggested the program to use a symmetric diagonalization (which are more stable) and results were much better. Pruning this little stuff out made them exact.
This is the first stuff I did and I both wanted to know if you encounter the same issue and solved somehow or if you leave it untreated, and prepare some groundwork for what I want to write about the lanczos convergence issue with krilovkit (I really hope to use Julia in future, so trying to solve bugs is really important to me)
I still don’t quite understand the issue you are purportedly seeing, I think you have to be a little bit more concrete, perhaps we can try to take a few steps back. First of all, are you trying to run DMRG with ITensorMPS.jl? Is the issue that you are seeing that the function ITensorMPS.dmrg does not converge to the proper ground state given a certain Hamiltonian, and you believe that the reason it isn’t converging is because the Krylov solver isn’t finding the correct eigenvector in one or more steps of DMRG?
Normally, floating point round-off errors (especially ones that are on the order of double precision) don’t lead to convergence issues in DMRG so I can’t say I’ve seen issues like the one you are describing, though I may not be fully understanding the issue you are reporting.
I never used Itensor, but if I’m going to do it, is when I’m really sure to choose Julia over C with MKL and arpack libs.
I just wrote down a DMRG for Ising and XXZ models with kron function, mul!, the standard Julia full eigen solver, krilovkit lanczos and sparse matrices and vectors.
The first issue I found was due to not pruning these artifacts out. Others languages manage them (different algorithm implementations) differently, so there was no need to prune them out.
Regarding convergence of krilovkit (which you use): I don’t know what really is the problem, as the converged results are right. However this is the first thing I did before encountering convergence issues: pruning out artifacts. I’m not rounding value of 0.7 + 1e^-17. I’m rounding values under e^-13 or e^-15, so what by eye you recognizes as clear artifacts. Results saturate when removing them, up to a certain threshold value, so even physics agrees with this.
I don’t know how you handle the sparse/full matrices/tensors/arrays/objects, so this could be one of the things which leads to different results. I know from a colleague of mine Itensor never gave him issues. However he stayed on fortran after trying it out and didn’t try the models I mentioned above.
Could you share the Hamiltonian (say as a formatted \LaTeX formula) you are trying to find the ground state of where these kinds of round-off issues affect the convergence of DMRG? I find that a bit surprising but I suppose it is possible, say in the case of a nearly degenerate ground states, in which case strategies like those outlined in DMRG FAQs · ITensors.jl are generally sufficient (or those situations are just fundamentally difficult for DMRG, since if your answer is sensitive to numerical round-off that indicates something about the fragility of your ground state to begin with).
If you are concerned about KrylovKit.eigsolve (which again, we have found to be very reliable), Julia has a wrapper for Arpack here: GitHub - JuliaLinearAlgebra/Arpack.jl: Julia Wrappers for the arpack-ng Fortran library, which you could compare against as a reference.
In ITensor, if you don’t conserve quantum numbers, we just use dense tensors (or diagonal matrices, in the case of storing singular values or eigenvalues). If you conserve quantum numbers, we store block sparse tensors, where the blocks are generally dense, but may be diagonal in the case of storing singular values or eigenvalues. Other than that, we just perform normal floating point arithmetic with those dense or diagonal tensors, and the only time we truncate small elements of tensors is when we perform truncated SVDs or eigenvalues, in which case we truncate the singular vectors or eigenvectors corresponding to small singular values or eigenvalues.
2 Likes
I will provide latex and a good code in the first days of July. I want it to be understandable, so I will even provide different runnable programs with plots in an article fashion.
I expect my previous results to be proved wrong. This way everything works fine 
I have some stuff I need to address first, so it may slide of a couple of months.
Sorry for the inconvenience, but a relaxed thread is the best