Exponentiation algorithms, roadmap

[jagot]

This issue collects the ideas that go into the roadmap for Expmv.jl. The basic problem is to calculate v = kappa(dt*A)*u), where kappa(z) is some approximation to exp(z). This should be done in-place, and a matrix exponential should never be calculated (unless explicitly requested).

Different exponentiation algorithms

• [ ] Padé approximants: these could be automatically generated from Polynomials.jl. For constant matrices, the denominator polynomials can be prefactorized using lufact!, or preconditioned, if using an iterative method for the solver step.
• [ ] Taylor methods; could also be automatically generated to arbitrary order. Some caution should be exercised.
• [ ] Krylov methods; maybe some help from IterativeSolvers.jl?. These methods also need a way to calculate the subspace approximation to the matrix exponential (Padé, Chebyshev, or full diagonalization).
  • [ ] Hermitian Lanczos
  • [ ] Implicitly restarted Arnoldi
• [ ] Chebyshev method via the Clenshaw algorithm (this requires a priori knowledge of the extent of the operator spectrum).

Base changes

It could be that for an operator L = A + B + C + ..., A is diagonal or at least easy to exponentiate in one basis, whilst B is simpler in another base. A prime example is H = T + V, where the kinetic operator T is diagonal in momentum space (being just k^2/2), whereas V is usually diagonal in position space. The basis transform is simply the Fourier transform (in Cartesian coordinates, in cylindrical coordinates, one wants to use the Hankel transform). Some nice way of tagging suitable bases for each operator and providing transforms, would be useful.

Traits

Caches

References

General

@ARTICLE{Adsuara2016,
  AUTHOR = {Adsuara, J.E. and Cordero-Carrión, I. and Cerdá-Durán, P. and Aloy, M.A.},
  PUBLISHER = {Elsevier BV},
  URL = {http://dx.doi.org/10.1016/j.jcp.2016.05.053},
  DATE = {2016-09},
  DOI = {10.1016/j.jcp.2016.05.053},
  ISSN = {0021-9991},
  JOURNALTITLE = {Journal of Computational Physics},
  PAGES = {369--413},
  TITLE = {Scheduled Relaxation Jacobi method: Improvements and applications},
  VOLUME = {321},
}

@ARTICLE{Alvermann2012,
  AUTHOR = {Alvermann, A and Fehske, H and Littlewood, P B},
  PUBLISHER = {IOP Publishing},
  URL = {http://dx.doi.org/10.1088/1367-2630/14/10/105008},
  DATE = {2012-10},
  DOI = {10.1088/1367-2630/14/10/105008},
  ISSN = {1367-2630},
  JOURNALTITLE = {New Journal of Physics},
  NUMBER = {10},
  PAGES = {105008},
  TITLE = {Numerical time propagation of quantum systems in radiation fields},
  VOLUME = {14},
}

@BOOK{Lubich2008,
  AUTHOR = {Lubich, Christian},
  PUBLISHER = {European Mathematical Society},
  DATE = {2008},
  ISBN = {3037190671},
  TITLE = {From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis (Zurich Lectures in Advanced Mathematics)},
}

Krylov methods

@ARTICLE{Saad1986SJosasc,
  AUTHOR = {Saad, Yousef and Schultz, Martin H},
  PUBLISHER = {SIAM},
  DATE = {1986},
  DOI = {10.1137/0907058},
  JOURNALTITLE = {SIAM Journal on scientific and statistical computing},
  NUMBER = {3},
  PAGES = {856--869},
  TITLE = {{GMRES}: A generalized minimal residual algorithm for solving nonsymmetric linear systems},
  VOLUME = {7},
}

@ARTICLE{Saad1992SJoNA,
  AUTHOR = {Saad, Yousef},
  DATE = {1992},
  JOURNALTITLE = {SIAM Journal on Numerical Analysis},
  TITLE = {{Analysis of some Krylov subspace approximations}},
}

[ChrisRackauckas]

Can we reserve the symbol phi to refer to the phi functions from exponential integrators? I mean the ones discussed in this review:

https://na.math.kit.edu/download/papers/acta-final.pdf

like phi1 = (exp(z)-1)/z etc. since these are common in exponential Runge-Kutta. I think we should include ways to efficiently compute these as well. ExpoKit's paper is a nice reference:

https://www.maths.uq.edu.au/expokit/paper.pdf

[ChrisRackauckas]

This should be done in-place, and a matrix exponential should never be calculated (unless explicitly requested).

This is key for my applications, where usually that's a large PDE discretization where the matrix exponential would produce a large dense matrix which doesn't even fit into memory. So yes, this should be a requirement for an expmv! function.

[ChrisRackauckas]

Note that we can target getting just one method at first so that way the integrators can be written, and leave alternatives as a future GSoC project. Some of these are great self-contained problems for a student.

[mforets]

Hello, just to mention that an expmv! function is provided in Expokit.jl. there is some current ongoing work to add more functions.

[ChrisRackauckas]

@mforets does ExpoKit.jl now have all of this?