This Issue has arise from this previous question MPO for Anyon Hubbard model, the user in question wrote a non hermitian hamiltonian by mistake, but he was able to notice that thanks to a warning from KrylovKit. When I tried to reproduce the error I didn’t get any warning, that’s because the library has changed the settings about the verbosity and now to enable such warning one should pass eigsolve_verbosity=KryilovKit.WARN_LEVEL, but it would enable a bunch of other warnings that are not useful.
Since it is not difficult to make errors, either on the paper or in the code, and write the wrong hamiltonian, I think it would be better to add a check previous to actually starting the algorithm. something on the line of
Change ishermitian::Bool=true to ishermitian::Union{Bool,Nothing}=nothing
add another keyword check_hermitian::Bool=true
In this way the user could bypass the check if he is absolutely sure about the hermitianity of the MPO. The combination of the 2 parameters would be of the form:
!check_hermitian && isnothing(ishermitian) && @error "ishermitian can't be nothing if check_hermitian is false"
if isnothing(ishermitian)
ishermitian=ITensorMPS.ishermitian(H)
elseif check_hermitian
@assert ishermitian==ITensorMPS.ishermitan(H)
end
This also shows other vulnerabilities, like the absence of a method ITensorMPS.ishermitian(::MPO) and of an utility function to compute the adjoint (since adjoint(M::AbstractMPS)=prime(M))
We’ve considered this in the past but I wasn’t convinced there was a cheap way to check if general MPOs are (approximately) Hermitian, say if they have large bond dimensions. As you mentioned, I think ideally the check would be in the Krylov solver.
I’m not convinced there is a cheap way too, for this reason I want to give the chance to bypass the check, but I still think it is necessary.
Krylov library gives a warning but it comes with many others and it gets lost. Other solutions would be to either work with LoggingExtras and filter out the warning, or ask to KrylovKit dev to raise the level at which that warning gets printed, like pass it as a warning but printing it on the error level
One other hypothetical thing we could do in the future would be to take the terms of an OpSum and define Hermitian conjugation rules on them at the symbolic level. Like Adag mapping to A and so on. Then applying these rules and checking whether the transformed OpSum has the same terms afterward or not.
I like this idea, it wouldn’t be completely general but probably a person competent enough to don’t construct the Hamiltonian from an OpSum would also be able to check things by himself. For the moment I’d say to add a sort of warning or note to the documentation that would warn the user that no check is performed (in theory it wouldn’t be strictly necessary since already the documentation says that ishermitian=trueassumes the MPO to represent an hermitian matrix, but I think about a middle/low competent user that might miss this)
Can you maybe elaborate a bit on what you mean with the other warnings? For KrylovKit, we have definitely been struggling to strike a balance between not spamming your terminal, and also still providing warnings. It turns out that it is very tricky to find the correct threshold values for fully generic eigenvalue problems, given that from the viewpoint of KrylovKit, we don’t actually know anything about the problem, except that it is defined by a function handle f.
For the current (v0.9) iteration, we were thinking that the default value of 1 should be good for most use-cases, while still allowing it to be explicitly turned off when you really are sure about that. IIRC the current warnings are:
warn if we detect complex eigenvalues
warn if we detect degenerate subspaces which can get the algorithms stuck (mathematically that is not well-defined)
warn if the solver did not converge in the number of iterations.
We are definitely happy with any and all feedback on that though!
Hi @lkdvos,
Good to have you here, I’ll try to explain the point hoping to don’t say anything wrong since I’m the last arrived.
I think the point is that we are doing many subsequent calls to eigsolve, if the user supplied an MPO that does not represent an hermitian operator, but mistakenly supplied ishermitian=true, we would get both the warning for not converging and the warning for complex eigenvalues, the problem is that given the high number of calls, and that complex eigenvalues are not always detected, while the not converging warn is always present, we obtain a lot of warnings due to the not convergence of the algorithm and the warning about the complex eigenvalues gets completely lost (consider O(10^2/10^3) calls).
@mtfishman another think that I’ve noticed is that even when the matrix is hermitian, we are requesting by default a tolerance of 10^{-14} with a eigsolve_maxiter=1, if think it would be more coherent to pass as default eigsolve_tol=1. and add a warning in case the user increases eigsolve_maxiter but does not decrease eigsolve_tol, something like
if eigsolve_maxiter>10 && eigsolve_tol>KrylovDefaults.tol[]
@warn "when requesting eigsolve_maxiter>10 it is suggested to set eigsolve_tol\geq $(KrylovDefaults.tol[])" eigsolve_maxiter eigsolve_tol
end
and at this point set the default verbosity to KrilovKit.WARN_LEVEL
@VinceNeede is your suggestion meant to try to decrease the number of warnings output by KrylovKit? I think the part you are missing is the krylovdim parameter, even with maxiter=1, if eager=true the tolerance also controls if the entire Krylov space is used during that iteration, and we don’t want to build the entire Krylov space if we don’t need to.
@lkdvos would you consider giving more fine-grained control over what warnings get output? I could understand if you didn’t want to since there are endless things like that one could potentially add, which adds code complexity, maintenance burden, etc.
@mtfishman Yes that was the aim. Anyway, I’m not sure I understood what you are saying, I’m going to copy paste what I think is the relevant part from eigsolve
if K == krylovdim || β <= tol || (alg.eager && K >= howmany)
U = copyto!(view(UU, 1:K, 1:K), I)
f = view(HH, K + 1, 1:K)
T = rayleighquotient(fact) # symtridiagonal
# compute eigenvalues
if K == 1
D = [T[1, 1]]
f[1] = β
converged = Int(β <= tol)
else
if K < krylovdim
T = deepcopy(T)
end
D, U = tridiageigh!(T, U)
by, rev = eigsort(which)
p = sortperm(D; by=by, rev=rev)
D, U = permuteeig!(D, U, p)
mul!(f, view(U, K, :), β)
converged = 0
while converged < K && abs(f[converged + 1]) <= tol
converged += 1
end
end
if converged >= howmany || β <= tol
break
elseif alg.verbosity >= EACHITERATION_LEVEL
@info "Lanczos eigsolve in iteration $numiter, step = $K: $converged values converged, normres = $(normres2string(abs.(f[1:howmany])))"
end
end
if K < krylovdim # expand Krylov factorization
fact = expand!(iter, fact; verbosity=alg.verbosity)
numops += 1
else ## shrink and restart
if numiter == maxiter
break
end
Now, alg.eager is defaulted to false and the dmrg doesn’t modify it, furthermore if tol is high, either converged>=howmanyor β <= tol would yield true, exiting from the loop without expanding the subsapce. So I’d would say it is the opposite, by selecting a low value of tol we are forcing the algorithm to construct the whole krylov subspace of size krylovdim.
You can test it by yourself by running a small dmrg with eigsolve_verbosity=KrilovKit.STARTSTOP_LEVELyou can see that with eigsolve_tol=1. most of the times number of operations = 1, signaling that the algorithm stopped before constructing the whole subspace. If instead you ask for a lower tolerance, you would get almost always a warning with number of operations = 3 which is the default krylovdim passed, which means that the whole subspace has been constructed cause it was not able to converge to the desired tolerance.
@mtfishman mhh the point is that we are asking for a very low tolerance but then discarding every information about if such tolerance has been reached. Probably a solution would be to ask @lkdvos if it would be possible to add another parameter to the alg_rulewhich would me a minimum number of operations
I think it is purposeful that we are requesting a small Krylov space and also setting a small tolerance (maybe it doesn’t need to by quite that small, but that’s something we would need to investigate). Normally we aren’t worried about reaching the tolerance and will saturate the Krylov space, but then sometimes we might reach the tolerance and ideally we wouldn’t build the entire space to avoid unnecessary operations (so ideally we would set eager=true). So I think the settings and behavior are basically what we want as long as we set eager=true, it seems to me the main issue is that we want to tell KrylovKit that we don’t generally care that we aren’t reaching the desired tolerance so it doesn’t print those warnings. (@miles please correct me if I’m wrong about this, is there a reason we haven’t set eager=true or is that just a mistake in ITensorMPS.dmrg?)
For the moment I have created a branch with a logger filter GitHub VinceNeede/ITensorMPS/FilterKrilovKitWarnings
Essentially eigsolve is always called with at least WARN_LEVEL of verbosity, but if eigsolve_verbosity=0 the warnings are filtered and passed only if the file that sent the warning ends with “factorizations/lanczos.jl”, additionally, the warning is shown only one time per sweep, to don’t spam the terminal. I have also added a test that would allow us to recognize if KrylovKit changed something and we are not able to filter the warning anymore.
This way of filtering is not bullet proof, but if you like the idea, I can try to create a pull request on KrilovKit where I’d add a custom id to the warning, allowing us to filter the warning via that specific id.
I’d still think at this as a temporary solution, KrylovKit doesn’t always recognize complex eigenvalues, or it does that after some sweeps (when I was trying I had to wait 20 sweeps), so I think it is still better to add a pre check on the MPO.
Anyway let me know what do you think about the filter idea.
@mtfishman I think the lack of eager=true in ITensorMPS is an oversight and we should set it to true by default. Probably neither of us appreciated that flag when we were putting that code together (I know I didn’t), though clearly you or someone else did catch its importance in ITensorNetworks (modulo the typo where it wasn’t being passed through).
That’s an interesting idea to filter the logger to only see the warning you want. Does this require modifying ITensorMPS.dmrg itself or can you do the same thing by wrapping that filter around the call to ITensorMPS.dmrg?
I’m not sure I fully understood your idea, could you elaborate better?
What I think you mean is to create a new step in the call stack of dmrg that would look like this:
function dmrg(...[usual parameters]...; ... logger=filter_krylovkit_warnigns()
res = nothing
with_logger(logger) do
res = dmrg(...;... eigsolve_verbosity = max(eigsolve_verbosity, 1)) #without passing logger
end
return res
If this is what you meant, it would be possible in principle, I’m just not really sure how to don’t make it spam the terminal and print it once per sweep, probably a solution would be to use a counter and make the message print every 2(N-1) times (where N=length(H)), which should be the number of iterations in a single sweep. It wouldn’t be exactly the same as printing once per sweep but it would be quite close. If we want to pass the warning at most one time per sweep then the logger has to be internal to the dmrg function, otherwise I can’t see how it would be possible to track if the warning has already been printed in that sweep
I think the custom logger is also what I would suggest. In principle, the logging system in Julia is designed for the log id to remain stable between different versions, but I think we could definitely agree to “hardcode” them in KrylovKit such that they become reliable, thus giving any user the option to have very fine-grained control over what messages are printed.
The main advantage of the logging system is that this can be handled at the level of dmrg, but a user can also simply set the global logger or overwrite the logger themselves, so it does not require any modification of package code.
All of this aside, one thing that might be interesting (and possibly worth rethinking in KrylovKit) is that we typically think of tol as the thing you want to achieve, while maxiter is a failsafe to avoid looping forever. In some sense, we consider reaching maxiter a “failed run”. This reason I’m saying this is that we scale some of our checks and warning messages etc by tol, and setting this to be 1 effectively destroys that behavior: you would no longer get warnings about degenerate subspaces either.
@miles and @mtfishman , do you think there is an alternative strategy where you simply set the tolerances for example to 1e-5 and don’t rely on the maxiter? That might align the messages more with the assumptions in KrylovKit too.
As an aside: up to some memory usage, you can very often get away by simply using Arnoldi, even for hermitian systems. It might not be unreasonable to just default to ishermitian=false, and only set that to true in cases where we can explicitly check hermiticity, for example when the Hamiltonian is a sum of all hermitian terms.
d = 2
D = 64
A = rand(d*D^2, d*D^2);
A += A'
x = rand(size(A, 1));
julia> @time eigsolve(A; ishermitian=true);
0.381823 seconds (181 allocations: 2.154 MiB)
julia> @time eigsolve(A; ishermitian=false);
0.391688 seconds (181 allocations: 2.154 MiB)
For time evolution using KrylovKit.exponentiate (inside of TDVP) I am experimenting right now with ways to speed it up and yet keep global accuracy. I’ve been a bit cautious about raising the tol parameter much so far but it’s the next thing I want to try. Since accuracy is crucial for time evolution I do ideally want to do as you say and keep things like maxiter high and control everything by tol.
On the other hand for DMRG, I think it’s more subtle because there is a tradeoff in the sense of: is it worth it to spend more time inside KrylovKit.eigsolve or just go on to the next site or sites and optimize it some more to continue to improve the basis? There is a lot of evidence that for easier problems (e.g. say if the gap is large) that such a strategy works well and one should just exit eigsolve after say only building two or three Krylov vectors.
But we have also seen cases where it’s better to run eigsolve a bit harder at each step. I think you see my point though that since DMRG has a global aspect and is just doing optimization one may not want to do high-accuracy eigensolving until the basis is good.