-
Notifications
You must be signed in to change notification settings - Fork 109
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
KrylovKit spectrum solver #367
Conversation
Codecov Report
@@ 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
📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for taking care of all of this, I greatly appreciate the contribution! Attached here are a few comments and suggestions. Let me know what you think and how you would like to proceed.
@@ -9,6 +9,7 @@ DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" | |||
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" | |||
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" | |||
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" | |||
KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@david-pl , I am in favor of adding KrylovKit
as a direct dependency -- it is a very mature and well supported Julia library, and an important part of the iterative linear algebra ecosystem. And we already depend on it because of SciML https://juliahub.com/ui/Packages/KrylovKit/JPeTr/0.6.0?page=2
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@Krastanov fine by me! @aryavorskiy Thanks for the work here.
src/spectralanalysis.jl
Outdated
DiagStrategy{T}(n::Int) where T = DiagStrategy{T}(n, nothing) | ||
const LapackDiag = DiagStrategy{:lapack} | ||
const ArpackDiag = DiagStrategy{:arpack} | ||
const KrylovDiag = DiagStrategy{:krylov} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
could you consider structuring this a bit differently:
abstract type DiagStrategy end
type LapackDiag{T} <: DiagStrategy
n::Int
v0::T
end
My main reasons for this suggestion:
- It is a bit closer in style to how SciML treats their solver types
- It is a bit closer in style to Holy Traits https://invenia.github.io/blog/2019/11/06/julialang-features-part-2/
- It parameterizes
v0
instead of leaving it of abstract type, which avoids boxing (but to be fair, you can do that to your version either way) - It is not mutable anymore
Please feel free to voice dissent to this suggestion (as always!)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I was wrong on the last one: it still needs to be mutable because you might change n
as I see below
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree, this is a better way. Initially I used this structure to make custom subtyping way easier, so that LapackDiag, ArpackDiag and KrylovDiag inherited the same internal structure.
Probably a good way to avoid the DiagStrategy
struct being mutable is making the DiagStrategy
constructor process the keyword arguments from the eigenstates
method. This explanation is probably kinda messy, will clarify in the nearest commit.
src/spectralanalysis.jl
Outdated
end | ||
DiagStrategy(m::Matrix) = LapackDiag(size(m)[1]) | ||
DiagStrategy(::SparseMatrixCSC) = KrylovDiag(6) | ||
DiagStrategy(m::AbstractMatrix) = ArgumentError("Cannot detect DiagStrategy for array type $(typeof(m))") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could you elaborate on this error message. For the last year or more SciML has undertaken a really great process of writing as helpful as possible error messages. For instance, could you add something along the lines of This error probably is raised because you have tried to ... To fix it, consider doing ...
src/spectralanalysis.jl
Outdated
|
||
function eigenenergies(op::AbstractOperator, n::Union{Int,Nothing}=nothing; kw...) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This seems like the main entry to most of the code paths. Could you consider reorganizing the docstrings a bit, because now they seem to be at least partially outdated. A few options:
- move the information from most docstrings to only here and explain how the dispatch works. Then add very short docstrings to the other methods
- declare a
function eigenenergies end
without a method, and add the docstring to it.
Could you also add a brief sentence to the appropriate docstring explaining what is happening and that the library tries to detect the best choice on its own. And that the choice can be overwritten by using the appropriate type.
src/QuantumOptics.jl
Outdated
@@ -7,6 +7,7 @@ using SparseArrays, LinearAlgebra | |||
export | |||
ylm, | |||
eigenstates, eigenenergies, simdiag, | |||
LapackDiag, ArpackDiag, KrylovDiag, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Given that the library mostly automatically decides which method to use, I think it would be best to not export these. If nonetheless we want them exported, they should at least have a docstring each.
@Krastanov thanks for the review! I will implement these requests by Tuesday, 11th, I think. I also have some new ideas:
|
|
Done. |
@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! |
I made two small changes to your PR:
I will wait on the tests, give it a final read through, and I should be able to merge it and release it today. |
My bad, it seems |
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. |
#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