Extending expv using Saad to general phiv?

[MSeeker1340]

@jagot I have finally went through the paper by Saad. While in theory the algorithm can be extended to arbitrary phi orders and non-Hermitian matrices, I still have some doubts and would like to discuss them with you.

Applicability to non-hermitian operators?

Our idea is to controll the Krylov subspace size not by a priori error estimates ("happy-breakdown" being the simplest version), but by using a posteri error estimates as the computation proceeds (btw, I don't think Saad put forward this idea himself; can I assume it's your invention then ;)).

Now a problem with this approach is that the error estimates should not be expensive. For hermitian operators this is fine because tridiagonal Hm is fairly cheap to work with (direct diagonalization); on the other hand if Hm is just upper Hessenberg then having to calculate exp(Hm) every step is very expensive, and we're better off just doing a few more Arnoldi iterations with a fixed m.

If you agree with this, then we can restrict ourselves to implementing the algorithm only for Hermitian operators.

Error estimates

Now the good news is that all of Saad's derivations (in particular, the error formula (42) and (43)) can be easily adapted to arbitrary phis, and we can get the phi versions of Er1 and Er2. In fact, Er1 is the same error estimate used by Niesen & Wright in their phipm algorithm.

phiv! already has Er1 implemented and it is used by phiv_timestep!. The interface is a bit ugly at the moment: since to get Er1 estimate for phi_k you need phi_(k+1), the way I went is actually computing the error estimate of the second-to-last phi using the last term (so, as an example, if you want to compute phi_3(tA)b you need to actually call phiv! with k=4; this is really inconvenient and I may rewrite it in the future). But nevertheless it is present, so creating a termination condition using Er1 estimate is very much doable.

The problem I have is Er2 estimate, which the current expv! with Saad is using. Replacing exp with phi_1 is... kinda ok? But we certainly cannot do the same for phi_1's error estimate, since phi_1 and phi_2 don't even converge to the same value. Something has to be done to generalize Er2, and I'm not sure how.

[jagot]

First off, I must say that I am in no way an expert with Krylov methods.

I do think that Saad intended this, since if you look in section 6, numerical experiments, he does tabulate the actual errors as well as the estimates, as a function of m. In any way, I learned the trick from a friend at Stockholm University who, incidentally, did work with non-Hermitian matrices and did diagonalize the Hessenberg matrix H_m at every m. Compared to calculation the action of H^m on v, the diagonalization of H_m is cheap. However, Saad recommends rational approximations to the subspace exponential (see section 2.4), which are fairly cheap to compute.

I agree that the Er2 estimate is only applicable to phi_0, but if you already have the better estimate implemented, do you want one corresponding to Er2 simply for speed?

[MSeeker1340]

I agree that the Er2 estimate is only applicable to phi_0, but if you already have the better estimate implemented, do you want one corresponding to Er2 simply for speed?

That's a good question. Short answer is I have no idea since I don't have a grasp of how much time different portions of the code requires. That's probably why Chris asked me to do a detailed benchmark in the first place (#12). After that we can probably have a more meaningful discussion.

And yeah using rational approximations is probably the best way to go for quick error estimates (and for small |Hm| getting a truncated series expansion for phi_k(Hm) is very easy).