Balancing for the 2009 Scaling & Squaring matrix exponential algorithm

[andreadelprete]

I start a new discussion about a balancing algorithm tailored for the 2009 scaling and squaring (S&S) matrix exponential algorithm, to avoid polluting the already existing issue where we discuss balancing for the 2005 algorithm.

[andreadelprete]

I've taken a quick look at the new 2009 version of the matrix exponential S&S algorithm. Basically, instead of choosing the number of squarings s based on the 1-norm of A, it is based on the 1-norm of A^k, with k = 4, 6, 8, 10. This suggests that an ideal balancing algorithm for this algorithm should try to minimize these norms. I'm not saying we should do it, for now I'm just trying to get ready for potential critics from the reviewers in this direction.

[andreadelprete]

Looking at the implementation of the algorithm in scipy, if we wanna focus on the hard cases where you end up using a Pade-13 approximant, the crucial part of the code is this:

    eta_3 = max(h.d6_tight, h.d8_loose)
    ...
    eta_4 = max(h.d8_loose, h.d10_loose)
    eta_5 = min(eta_3, eta_4)
    theta_13 = 4.25

    # Choose smallest s>=0 such that 2**(-s) eta_5 <= theta_13
    if eta_5 == 0:    # Nilpotent special case
        s = 0
    else:
        s = max(int(np.ceil(np.log2(eta_5 / theta_13))), 0)

where dk is ||A^k||^(1/k). For k=6 the exact (i.e. tight) value is computed (because anyway A^6 must be computed to evaluate the Pade approximant), whereas for k=8, 10 they use an approximated (i.e. loose) value (because A^8/10 need not be explicitly computed). In the end, the quantity that a balancing algorithm should minimize is eta_5 = min( max(d6,d8), max(d8,d10) ). This seems extremely challenging at first sight because we first need to figure out which of these three values (d6, d8, d10) we really care about, with the additional challenge that for d8 and d10 we need to work with approximations because we don't have the corresponding matrix.

A possible strategy could be to simply minimize d6, which is the only one for which we have access to the corresponding matrix power. We should look for a diagonal similarity transformation D to apply to the matrix A:

B = D^-1 * A * D

such that the resulting d6 of B is minimized:

|| B^6 || = || D^-1 * A * D * ... * D^-1 * A * D || = || D^-1 * (A^6) * D ||

So we should first compute A^6, and then apply the matrix balancing to it. Therefore the balancing algorithm should work inside the S&S algorithm. This is rather easy to do, because it doesn't affect the balancing algorithm per se, we just need to apply it to A^6 rather than to A. We know that minimizing d6 does not necessarily minimize eta_5, but it seems better than minimizing ||A|| for this S&S algorithm.

[andreadelprete]

instead of choosing the number of squarings s based on the 1-norm of A, it is based on the 1-norm of A^k, with k = 4, 6, 8, 10
That's the conclusion I've reach too but if that norm calculation takes a lot of time I don't see the point

The key point is that inside the S&S algorithm you need anyway to compute A^6, so balancing A or balancing A^6 takes exactly the same computation time. The same is not true for A^8 and A^10 though, which would require an additional matrix-matrix multiplication.

[olimexsmart]

If I got it right, are we saying we should basically take the code of the current scipy.linalg.expm and modify it including inside the balancing of A6?

A6 can be balanced with the same logic we worked so far? i.e. minimizing the columns?

This also means we need to implement the 2009 algorithm in C++ sooner or later

[andreadelprete]

No, sorry, there was a misunderstanding. This discussing about the 2009 S&S algorithm is hypothetical. I'm not suggesting you to implement this. We are just considering this option for the future. Most likely, it's just something that we can discuss in the "future work" section of the paper.

[olimexsmart]

Ok that's reassuring, now I am sure on what to do