KrylovKit spectrum solver

[aryavorskiy]

364

The Arnoldi algorithm implementation from the KrylovKit.jl package seems to be much more performant than the ARPACK one. Also a DiagStrategy trait was implemented to make selecting the correct method on the fly (or defining a new one) more flexible

[codecov[bot]]

Codecov Report

Merging #367 (b9dca8c) into master (eff22f1) will decrease coverage by 0.30%. The diff coverage is 88.37%.

@@            Coverage Diff             @@
##           master     #367      +/-   ##
==========================================
- Coverage   98.11%   97.82%   -0.30%     
==========================================
  Files          18       18              
  Lines        1488     1516      +28     
==========================================
+ Hits         1460     1483      +23     
- Misses         28       33       +5     

Impacted Files	Coverage Δ
src/spectralanalysis.jl	93.42% <88.37%> (-6.58%)	:arrow_down:

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

[aryavorskiy]

@Krastanov thanks for the review! I will implement these requests by Tuesday, 11th, I think. I also have some new ideas:

• Maybe remove the Arpack diagonalization routine? KrylovKit is much faster (also it is universal and can be applied even to CUDA arrays). Benchmarks attached:

  
  julia> using Arpack, KrylovKit, SparseArrays, BenchmarkTools

julia> A = sprand(1000, 1000, 0.001); v0 = rand(1000);

julia> @benchmark eigs(A, v0=v0, nev=5, which=:SR) # Arpack BenchmarkTools.Trial: 589 samples with 1 evaluation. Range (min … max): 7.974 ms … 12.672 ms ┊ GC (min … max): 0.00% … 33.56% Time (median): 8.431 ms ┊ GC (median): 0.00% Time (mean ± σ): 8.492 ms ± 514.293 μs ┊ GC (mean ± σ): 0.64% ± 3.72%

 ▄▆▆█▅▃

▃▄▇██████▆▄▄▃▃▂▁▂▂▂▁▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▂▁▁▂ ▃ 7.97 ms Histogram: frequency by time 12.3 ms <

Memory estimate: 702.83 KiB, allocs estimate: 949.

julia> @benchmark eigsolve(A, v0, 5, :SR) # KrylovKit BenchmarkTools.Trial: 2167 samples with 1 evaluation. Range (min … max): 1.931 ms … 6.259 ms ┊ GC (min … max): 0.00% … 50.37% Time (median): 2.162 ms ┊ GC (median): 0.00% Time (mean ± σ): 2.302 ms ± 457.034 μs ┊ GC (mean ± σ): 2.63% ± 8.01%

▄█

▃▇███▇▇▆▆▅▄▄▃▃▄▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▂▁▁▂▁▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂ ▃ 1.93 ms Histogram: frequency by time 4.52 ms <

- Maybe create an `EigenSystem` struct that will store computed eigenvalues and eigevectors? In this format they can be easily dumped to a file or used to construct a density matrix or Green's function.
- Maybe move this functionality to QuantumOpticsBase? Other quantum physics packages depending on it might want to use this as well. I am myself now in progress of creating such one :-) Or is this a generally bad idea?

[Krastanov]

• It seems relatively conclusive that KrylovKit is better. I am comfortable with the idea of removing arpack. Here are other comments from the community https://discourse.julialang.org/t/suggestions-needed-diagonalizing-large-hermitian-sparse-matrix/96580
• I personally do not think an EigenSystem type would be better than simply an array or tuple of results. My reasoning is that an extra type makes sense only if we need some API on top of it, i.e. a bunch of useful methods that are more ergonomic to use with EigenSystem than with a vector or tuple of results. For the serialization question: I think a Tuple is as easy or easier to serialize than a custom type. Do not hesitate to elaborate on an API for EigenSystem that would be useful for you -- that would be a good reason to reverse my opinion.
• Could you share a bit more about the package you are constructing? Separation like that might make sense (and has happened before):
  • QuantumOpticsBase was split out from QuantumOptics, so that there is a base library for data structures and a separate library for dynamics over these datastructures (and the eigensolver is kinda in between, but in my mental model a bit closer to QuantumOptics than to QuantumOpticsBase)
  • QuantumInterface was split away from QuantumOpticsBase, so that the abstract types and bases are in an interface library that can be reused by simulators that do not care about state vectors and density matrices and master equations. E.g. thanks to that split, the QuantumClifford library for simulating Clifford circuits and the QuantumSymbolics library for symbolic expressions can share an API with QuantumOptics without depending on it
  • As you are contributing more to this ecosystem, we should consider how to make it easier for you, so some rearrangement might make sense

[aryavorskiy]

Done. One final idea: I don't like how two separate keywords are used for nonhermitian warning and info message about the diagonalization method. My proposal is using a verbose keyword for both. This change is breaking, but only a little, so I don't know if it is a good idea. @Krastanov about my package: it is called LatticeModels.jl, it is designed to perform various computations on tight-binding lattices. I decided to use the bases & operators infrastructure from QuantumOpticsBase when I discovered that it is a separate package, and am currently working on its integration. My package is still far from finished, I will definitely rename it sometime, the docs are a little outdated, but most important functionality is already implemented.

[Krastanov]

@aryavorskiy , there was a minor merge conflict that I resolved, but that ended up pushing to your master branch. Be sure to do a pull before you perform other edits. Also, it is a good idea to do your development not to the master branch of your repo, but rather on a separate development branch -- it makes other modifications down the line easier.

I will go through the current changes and review. Thanks for improving this part of the library!

[Krastanov]

I made two small changes to your PR:

• The strategies are not exported anymore. I try to be very minimal with what is exported, because the convention is that exported symbols will be supported in the future, which makes refactoring more difficult.
• I turned your info keyword argument in a @debug statement as this is a more universal julia-wide way to track this type of user messages. That way any function that provides such stdout information can be controlled in the same way (instead of the user needing to learn the conventions of every single library)

I will wait on the tests, give it a final read through, and I should be able to merge it and release it today.

[Krastanov]

My bad, it seems info was something already done to this library so you were just following the established convention. Give me a minute to fix this.

[Krastanov]

It seems the info string was meant to specifically warn about defaulting to the sparse eigsolver when the dense one might be better. I moved it to that location and cleaned it up a bit.