Hessenberg matrix size in Krylov subspace in the expokit routine

[Vipin Varma]

Dear expokit users/developers,

I am currently using a C++ version of the expokit's Matlab code for performing unitary time evolution of large sparse matrices, with linear dimension between 1000 and 1 million; the routine is proving very helpful for my tasks.

However the following two points regarding the size of the Hessenberg matrix are not fully clear to me; any inputs will be appreciated:

1) The online html documention explains that the H matrix being exponentiated is an (m+1)x(m+1) size matrix, where m is the size of the Krylov subspace. But the codes [and the pseudocode in the original paper] implement the exponentiation of an (m+2)x(m+2) matrix with H(m+2,m+1) = 1. Could you please explain what is the point of initialising this element to 1, and exponentiating this larger (m+2)x(m+2) matrix instead of the smaller (m+1)x(m+1) matrix? Perhaps this is for better numerical stability?

2) My understanding is that the time-step can be arbitrarily large for a given Krylov subspace dimension, say m=30; but then the variable mxrej would have to be increased so that the required tolerance can be met. Does that sound right? I ask because I see a couple of papers which gurantee an upper bound on the error only if the Krylov subspace size and the norm or eigenvalues of the matrix satisfy some relation.

Thanks in advance.

Best,
Vipin

[expokitfan]

Dear Vipin,

1) You are quite sharp to have noticed that line 

"H(m+2,m+1) = 1"

 in the Alg. 3.1 pseudocode of the original paper!  :)  If you read further on in the same paper, you'll also find the answer to your question: the vector c* of Corollary 1 of Theorem 1 is the culprit, having been set to e_{m+1} in the pseudocode by means of the line "H(m+2,m+1) = 1".   I suppose the missing h_{m+1,m} factor in c* alluded to in the text following Corollary 1 is somehow  taken care of after the exponentiation of H_{m+1}.   So indeed, all this is required for stability, as you rightly guessed.

 

2) I'm afraid I cannot figure out what you mean by "the variable mxrej".  A reference to exactly where it is mentioned in the original paper, or elsewhere, might be more helpful.

Cheers,

expokitfan 

[Vipin Varma]

Dear expokitfan,

(i) Great, makes sense; that was very helpful.

(ii) That was sloppy of me: I was referring to the maximum number of internal time-step rejections permitted before an error is thrown out; it is named mxrej in the Matlab code and mxreject in the Fortran code. My understanding, which may be incorrect, is the following: there are mathematical relations, as described in certain papers, between the Krylov subspace size and the properties of the matrix A so that a given tolerance can be met in the matrix exponentiation. But in expokit it seems that this problem is circumvented by the time-stepping procedure, which will ensure that the tolerance will always be met independent of how large the time 't' requested is, even if the Krylov subspace is kept fixed at some low and constant value, say m=30. That is, for larger requested times 't' and fixed 'm' we just need to increase the mxrej/mxreject variable to a larger value. Does that analysis sound right?

Thanks!

[expokitfan]

Dear Vipin,

You're welcome.

(ii)  So I googled a matlab version on the web, and your reasoning seems sound to me, since increasing mxrej will enable further attempts to adjust the step-size so that the local truncation error err_loc meets the tolerance. I guess the matlab expm routine would prefer a not too large Hessenberg matrix as an argument, so it's convenient to put an upper bound on m.  Nevertheless, the original paper suggests yet another way (apparently not implemented in the code) to adjust m, if the need arises (see page 11, just before section 4).  But see the discussion after Equation 12 for a perhaps more complete answer to your queries, especially regarding the spectrum of the matrix.

Cheers,

expokitfan

[animesh....@gmail.com]

Hi Vipin,

I am also trying to use expokit for unitary time evolution of large sparse matrix. In all the examples they exponentiated hermitian matrix. How do I exponentiate i*hermitian matrix which is itself anti-hermitian. I shall be thankful for any help.

best
animesh

[Vipin Varma]

I have been using the generalised version of the matrix exponentiation i.e. using Arnoldi procedure instead of Lanczos. That should handle anti-Hermition as well.