Skip to content
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

WIP: Some improvements for Hermitian operators. #1

Closed
wants to merge 0 commits into from

Conversation

jagot
Copy link
Contributor

@jagot jagot commented Sep 27, 2018

I have made some improvements to the exponentiation for the (skew-)Hermitian case:

To consider:

  1. PR install make target JuliaLang/julia#6 Stopping criterion à la Saad; when using the Krylov iteration to approximate the exponential, rather than relying on the happy breakdown to occur (which I found to be rather unstable, i.e. it is very easy to get a too large, ill-conditioned Krylov space), it is better to successively calculate the exponential and use one of the "Practical error estimates" found in section 5.2 of Saad (1992). This is what I did in https://github.com/jagot/Magnus.jl/blob/master/src/lanczos.jl#L39 (I don't have the nice infrastructure as is available here).
  2. Operator/vector agnosticism à la KrylovKit.jl; specifically, I am interested in having a state vector that is not most optimally represented as an AbstractVector. If one could leave it up to the user to provide function for operator–vector multiplication, dot products, etc, I think it would be very useful.
  3. PR install make target JuliaLang/julia#6 Caching of subspace exponential scratch space: if one uses these exponentiators for time-stepping, one would prefer not to reallocate scratch space necessary for diagonalizing the subspace operator. In the Hermitian case, a real SymTridiagonal matrix is diagonalized, which under-the-hood calls dstegr from LAPACK. Looking into the Julia standard library, every call to the dstegr wrapper calls dstegr twice, the first time to query for necessary workspace size and allocating this. I hacked around this is Magnus.jl sub_exp.jl, but this should be tackled in a nicer way. I did open a ticket, but I never did much about it 😳
  4. We should get rid of opnorm, since it can be tricky to implement for complicated linear maps and costly to calculate, especially if you want to repeatedly exponentiate as you do in time-stepping. Admittedly, one could send in a bogus function as a kwarg that always returns a precomputed norm, but I don't really like it. This is conjunction with the stopping criterion in 1.

Comments and ideas are welcome.

@codecov
Copy link

codecov bot commented Sep 27, 2018

Codecov Report

Merging #1 into master will increase coverage by 0.01%.
The diff coverage is 96.82%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master       JuliaLang/julia#1      +/-   ##
==========================================
+ Coverage   96.85%   96.87%   +0.01%     
==========================================
  Files           4        4              
  Lines         191      192       +1     
==========================================
+ Hits          185      186       +1     
  Misses          6        6
Impacted Files Coverage Δ
src/phi.jl 100% <100%> (ø) ⬆️
src/krylov_phiv.jl 93.33% <93.33%> (ø) ⬆️
src/arnoldi.jl 93.75% <93.47%> (+0.13%) ⬆️
src/krylov_phiv_adaptive.jl 98.73% <98.71%> (ø) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update cedf70f...eb68cf6. Read the comment docs.

@ChrisRackauckas
Copy link
Member

We should get rid of opnorm, since it can be tricky to implement for complicated linear maps and costly to calculate, especially if you want to repeatedly exponentiate as you do in time-stepping. Admittedly, one could send in a bogus function as a kwarg that always returns a precomputed norm, but I don't really like it. This is conjunction with the stopping criterion in 1.

Before we do this, we should start benchmarking and profiling to know if this is the problem.

I think it would be best if each of these points are different PRs just to keep things clean. Definitely worth considering.

@ChrisRackauckas
Copy link
Member

Note: please do not use Project and Manifest TOMLs for now. These are not compatible with the current METADATA repository.

@ChrisRackauckas
Copy link
Member

Changing spacing should be a different PR. From a quick check of the diff I have no idea what code was actually changed.

@jagot
Copy link
Contributor Author

jagot commented Sep 29, 2018

Ok! I will split it into a couple of PRs starting with JuliaLang/julia#2

@MSeeker1340
Copy link
Contributor

Thanks for the work! Some comments:

Allowing a complex time-step, which enables exp(-imdtH)*psi, where H is Hermitian.

This would be very nice to have indeed. My original plan is to add support for skew-Hermitian H and keep time steps real, but if this is workable then great.

Stopping criterion à la Saad

Having a switch between this and happy-breakdown would be great. If I remember correctly Neisen & Wright also uses happy-breakdown in their matlab code so keeping the old way can make more meaningful comparison/benchmarking.

In the Hermitian case, a real SymTridiagonal matrix is diagonalized, which under-the-hood calls dstegr from LAPACK. Looking into the Julia standard library, every call to the dstegr wrapper calls dstegr twice, the first time to query for necessary workspace size and allocating this.

The other side of this issue is that for non-Hermitian cases, exp! is used instead of diagonalization which is also a Julia function that allocates every time. You can also take a look at this if you have the time ;)

We should get rid of opnorm, since it can be tricky to implement for complicated linear maps and costly to calculate, especially if you want to repeatedly exponentiate as you do in time-stepping. Admittedly, one could send in a bogus function as a kwarg that always returns a precomputed norm, but I don't really like it. This is conjunction with the stopping criterion in 1.

When we started out this summer vector and operator norms were not yet separated in Julia. Back then I used a custom norm function for both these roles (copying what Expokit.jl did), and during the v0.7 update we decided that custom vector norms aren't really necessary so I just changed norm to opnorm and restricted it to operators.

I agree with your comment that this can be inconvenient and since the operator doesn't really change we're better off passing in its norm as a value. I think in the kwargs we can default it to nothing and compute using LinearAlgebra.opnorm if no value is passed in.

@jagot
Copy link
Contributor Author

jagot commented Oct 1, 2018

Allowing a complex time-step, which enables exp(-im_dt_H)*psi, where H is Hermitian.

This would be very nice to have indeed. My original plan is to add support for skew-Hermitian H and keep time steps real, but if this is workable then great.

Yeah, so I will make a separate PR for this, once we are finished with JuliaLang/julia#4.

Stopping criterion à la Saad

Having a switch between this and happy-breakdown would be great. If I remember correctly Neisen & Wright also uses happy-breakdown in their matlab code so keeping the old way can make more meaningful comparison/benchmarking.

My latest idea about this is that we should have two separate expv! interfaces, one that does the full Krylov iteration (until happy breakdown) and then does the exponentiation, and another that tries to exponentiate early and stops when Saad is happy. To this end, I would factor out the Krylov step into its own function (arnoldi_step!/lanczos_step!), which is called in a loop from arnoldi!/lanczos! for the first interface, and directly from expv! for the second interface.

In the Hermitian case, a real SymTridiagonal matrix is diagonalized, which under-the-hood calls dstegr from LAPACK. Looking into the Julia standard library, every call to the dstegr wrapper calls dstegr twice, the first time to query for necessary workspace size and allocating this.

The other side of this issue is that for non-Hermitian cases, exp! is used instead of diagonalization which is also a Julia function that allocates every time. You can also take a look at this if you have the time ;)

Well, about this one, I have less of an idea what we can do to improve caching, but that the Hermitian case is non-allocating is vital to me.

We should get rid of opnorm, since it can be tricky to implement for complicated linear maps and costly to calculate, especially if you want to repeatedly exponentiate as you do in time-stepping. Admittedly, one could send in a bogus function as a kwarg that always returns a precomputed norm, but I don't really like it. This is conjunction with the stopping criterion in 1.

I agree with your comment that this can be inconvenient and since the operator doesn't really change we're better off passing in its norm as a value. I think in the kwargs we can default it to nothing and compute using LinearAlgebra.opnorm if no value is passed in.

If one uses the stopping criterion proposed by Saad, the norm of the operator is not so important, but rather the norm of the subspace vectors (the betas). Hence why I prefer that stopping criterion.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants