Add back restarts to applyExp

[mtfishman]

There may be some physically relevant situations where the Krylov dimension required to reach a certain timestep may be too large for the current implementation of applyExp. The current implementation may fail to converge because of orthogonality problems among the Krylov vectors, since we do not have reorthogonalization or restarts (for example in this support question). It would be good to add back restarts to make the implementation more robust. @mingruyang, what was the restart strategy used in the previous version?

I was thinking that a strategy could be to add a flag like MaxKrylovDim. If MaxKrylovDim iterations is reached and the error is not satisfied, one can decrease the timestep and recalculate the error, and keep doing so until the error goal (or some fraction of it) is reached. Then, one can restart the calculation and attempt to reach the remaining time, and repeat the same procedure. Is that roughly what was being done? I suppose the main question is the best procedure for figuring out what partial timestep should be attempted when the full timestep fails to be reached.

[mingruyang]

Sorry for the late reply.

Yes, that is roughly what was being done.

My original implementation has included a procedure to figure out the partial time step according to this paper. First we build a Krylov subspace of dimension max_iter (or less if beta_tol has been reached), and then do the exponential and estimate the errors, if the error goal is not reached, then adjust the time step to a smaller value and repeat the above procedure.

I tested this old implementation and found its performance issue is due to the inaccurate estimate of norm of the operator to be exponentiate. Since the old implementation will adjust the time step according to this norm before start, this inaccuracy will give a too small adjusted initial time step. An accurate estimate will be too expensive, so I now simply discard this initial time step adjustment and set the initial t_new = 1.0, which is also what Jutho did in the expintegrator of the KrylovKit. After fixing that, by setting the max_iter equal to the dimension of the Krylov subspace which is actually used in your new implementation and matching the tol, I got the same performance for both implementations.

There are two major difference of my implementation and yours.

The first is that my implementation only do the exponential of the small matrix after building the Krylov subspace (Jutho did the same way), but yours do that after each iteration. Unlike the Davidson, I think this is not necessary for the applyExp. Instead, if a smaller max_iter is already enough, one can simply adjust the max_iter to get less computation time.

Another difference is the definition of tol. In my implementation, tol = err/time_step. In yours, tol = err. So when comparing the performances, one needs to match this two definition.

I think after my fixing, the old implementation works fine (but of course need to be made more compact like what you did in your new implementation). The time-step adjustment can be used automatically if the results after building a Krylov subspace of max_iter is not satisfying. But we can output a warning message telling that the time-step adjustment is used and if this slow-down the calculation the user can try to set the max_iter larger. Also we can set a maximum allowed number of times the Krylov subpace can be restarted, just like Jutho did.

A further note is that the Lanczos often has a more serious issue of loss of orthogonality than Arnoldi. We can bring back the Arnoldi, which has more vector-vector multiplication but the same number of matrix-vector multiplication, so their complexity differs little.

How do you think?