Balancing & time scaling

[andreadelprete]

Time scaling is currently giving great results, showing a gain of about 40% in computation time, but the way in which we are tuning the scaling parameter T is not general. It would be nice to design an algorithm to find T automatically and that can work in any case.

BALANCING

Moreover, time scaling corresponds to a special case of balancing, that consists in using a similarity transformation given by a diagonal matrix D, so that the norm of B = D^-1 * A * D is much smaller than the norm of A. Then rather than computing directly e^A, we compute it as:

e^A = D * e^B * D^-1

Time scaling corresponds to setting diagonal elements of D to 1 for the first half and to T for the second half. Clearly this is suboptimal because we are using only one degree of freedom, whereas D has many more parameters.

The problem of balancing for the matrix exponential has been discussed in this mathworks post, coming to the conclusion that he didn't know whether to recommend balancing in general, and that many problems could be solved just by avoiding overscaling using a new version of the scaling-and-squaring algorithm from 2009.

BALANCING FOR THE MATRIX EXPONENTIAL

The balancing algorithm of LAPACK that was used for the matrix exponential preconditioning, was actually originally thought for eigenvalue problems. This algorithm minimizes the Frobenious norm of B, which is not really what we care about when computing matrix exponentials. It'd be better to design a new balancing algorithm that minimizes the 1-norm of B, which is what affects the computation of e^B.

[andreadelprete]

Reading better this paper I see that the current balancing algorithm of LAPACK (called GEBAL) actually minimizes the 1-norm of the vectorized matrix, that is the nXn matrix represented as an n^2 vector. This is not exactly what we need for the matrix exponential (which would be to minimize the max absolute column sum), but it seems rather close. We should investigate if we could modify the balancing algorithm to minimize this norm.

[olimexsmart]

I just finished reading this last paper you shared, and I have to say I found it very informative.

What are your toughs on chapter six? It seems to me that it reaches the same conclusions of the mathworks post

Currently I am tinkering with the balancing algorithms in MATLAB

[andreadelprete]

I think they're not considering that balancing can be much faster than squaring because you multiply times diagonal matrices rather than dense matrices.

[andreadelprete]

I've written a small python script to compare balancing with time scaling. Unfortunately they seem to give very similar results...

[olimexsmart]

I just finished to finally stabilize in C++ Algorithm3, in reference to the paper. But I needed to make a modification: the last row of matrix A1 are just zeros, hence any norm of that row is zero. This means that the execution gets stuck in the while loop at row 10. I added a condition to the inner while loops checking if the norm is zero.

I also tried to implement Algorithm1, but it is still unstable. Possibily because the norm should be computed disregarding the diagonal element, but it's not an immediate operation with Eigen.

For now the results are these:

*** PROFILING RESULTS [ms]                                  min        avg        max        lastTime   nSamples   totalTime  ***
BalancedComputeIntegralXt                                   0.068      0.075      0.211      0.072      5000       375.855   
BalancedComputeIntegralXtTV                                 0.047      0.052      0.157      0.050      5000       260.103   
ComputeIntegralXt                                           0.090      0.098      0.235      0.095      5000       490.439   

*** STATISTICS (min - avg - max - last - nSamples - total) ***
BalancedComputeIntegralXt                                   5.00        5.00    5.00    5.00    5000    25000.00
BalancedComputeIntegralXtTV                                 5.00        5.00    5.00    5.00    5000    25000.00
ComputeIntegralXt                                           13.00       13.00   13.00   13.00   5000    65000.00
ERROR2___BalancedComputeIntegralXt                          0.00        0.00    0.00    0.00    5000    0.03
ERROR2___BalancedComputeIntegralXtTV                        0.00        0.00    0.00    0.00    5000    0.03
ERRORINF_BalancedComputeIntegralXt                          0.00        0.00    0.00    0.00    5000    0.02
ERRORINF_BalancedComputeIntegralXtTV                        0.00        0.00    0.00    0.00    5000    0.02

Yes, they are similar to the timescaling, but the overhead introduced by the balancing is not negligible and the number of squarings is not lower.

[andreadelprete]

I just finished to finally stabilize in C++ Algorithm3, in reference to the paper. But I needed to make a modification: the last row of matrix A1 are just zeros, hence any norm of that row is zero. This means that the execution gets stuck in the while loop at row 10. I added a condition to the inner while loops checking if the norm is zero.

This is weird. I've used the balancing algorithm of scipy, which should be just a wrapper on LAPACK’s xGEBAL routine, and it was able to handle to row of zeros without any problem.

Yes, they are similar to the timescaling, but the overhead introduced by the balancing is not negligible and the number of squarings is not lower.

How much time does balancing take? If the results of balancing and time scaling are comparable, the only advantage of TS could be its computation time.

[olimexsmart]

Balancing takes about 3uS, which is coherent with the slight increase in time we are seeing in respect to TS

*** PROFILING RESULTS [ms]                                  min        avg        max        lastTime   nSamples   totalTime  ***
BalancedComputeIntegralXt                                   0.067      0.071      0.161      0.071      5000       356.878   
BalancedComputeIntegralXtTV                                 0.046      0.049      0.138      0.048      5000       246.356   
ComputeIntegralXt                                           0.089      0.094      0.219      0.094      5000       470.893   
DeltaX___TSON__TVON                                         0.040      0.046      0.118      0.045      5000       229.262   

*** STATISTICS (min - avg - max - last - nSamples - total) ***
BalancedComputeIntegralXt                                   5.00        5.00    5.00    5.00    5000    25000.00
BalancedComputeIntegralXtTV                                 5.00        5.00    5.00    5.00    5000    25000.00
ComputeIntegralXt                                           13.00       13.00   13.00   13.00   5000    65000.00
ERROR2___BalancedComputeIntegralXt                          0.00        0.00    0.00    0.00    5000    0.03
ERROR2___BalancedComputeIntegralXtTV                        0.00        0.00    0.00    0.00    5000    0.03
ERROR2___DeltaX___TSON__TVON                                0.00        0.00    0.00    0.00    5000    0.03
ERRORINF_BalancedComputeIntegralXt                          0.00        0.00    0.00    0.00    5000    0.02
ERRORINF_BalancedComputeIntegralXtTV                        0.00        0.00    0.00    0.00    5000    0.02
ERRORINF_DeltaX___TSON__TVON                                0.00        0.00    0.00    0.00    5000    0.02
SQUARINGS_DeltaX___TSON__TVON                               4.00        4.99    5.00    5.00    5000    24925.00

Sometimes the squarings drop to 4 too with TS

I may try to change the Beta parameter, or try to fix Algorithm1. Maybe a more accurate choice of T in our TS could gain us a definitive advantage.

[andreadelprete]

3 us is really little. It's about 6% of the total computation time (49 us). Maybe we could figure out a new version of balancing that is i) customized for the matrix exponential (tries to minimize the number of squarings) and ii) exploits the structure of our matrix to be faster (the top half of our matrix contains only zeros and ones).

What's Beta?

[olimexsmart]

What's Beta?

That would be the radix base used in Algorithm3. At the moment is set at 2 as suggested by the paper, because multiplications and divisions by powers of 2 intrinsically do not produce numerical errors (they are just shifts). Increasing this value should make the algorithm converge faster.

[andreadelprete]

I think it's better to leave Beta at 2.

[andreadelprete]

Hi @olimexsmart , in my latest commit bb74e36 you can find a new balancing algorithm, which (tries to) minimize the 1-norm of the given matrix. It's a more advanced version of the one you implemented in Matlab, and it is described in detail in this pdf.

Balancing for the matrix exponential.pdf

The algorithm seems to work well, in the sense that it converges in a reasonable number of iterations (always less than n^2, where n is the size of the matrix).

The script I've pushed also contains a comparison between the new balancing, the LAPACK balancing and the time scaling approach.

• In most cases the three approaches result in the same norm of A
• Sometimes time scaling performs a bit worse, and that could result in an additional squaring when computing the matrix exponential
• Sometimes the new balancing results in a slightly smaller norm but that doesn't affect the number of squarings

The tests I've carried out are rather limited, so I believe that a more exhaustive benchmark is necessary to really assess the new balancing algorithm. The fact that sometimes the new balancing manages to reduce the norm a bit more than the old one seems promising. Maybe there exist some matrices for which the new balancing would be beneficial (and we might have to make an effort to find them for the paper).

Another aspect to analyze is the computational cost of the new algorithm compared to the old one. I also believe that this new algorithm could be optimized. In particular, two possible strategies are evaluated at each iteration, but then in practice I've seen that the first one is always picked (which actually corresponds to what you had originally implemented). We should try to understand if that's always the case, so we could save some useless computations.

[andreadelprete]

I've also implemented the "classic" balancing algorithm (Parlett and Reinsch) following algorithm 3 of this paper to understand better how it differs from the one I've implemented earlier today. I'm not so sure it's exactly what's implemented in LAPACK though because it had problems due to a row of zeros in the matrix to balance, so I had to add some checks to handle this case. I think we should go and see exactly what's implemented in LAPACK.

[olimexsmart]

I had to add some checks to handle this case

Same problem I was referring to here

But I needed to make a modification: the last row of matrix A1 are just zeros, hence any norm of that row is zero. This means that the execution gets stuck in the while loop at row 10. I added a condition to the inner while loops checking if the norm is zero.

[andreadelprete]

Matrix balancing - Case study.pdf
I've found some cases in which minimizing the norm of the matrix (i.e. our algorithm) leads to better results than minimizing the norm of the vectorized matrix (i.e. the standard balancing). The most interesting case seems to be for this 2x2 matrix:

A = [[M, 1], 
     [M, 1]]

where M must be a large value (e.g., at least 10^6). I haven't verified this with the code, but the math seems solid. The next steps are:

• verify that the cases I've found are really interesting (i.e. we get better results with our algorithm)
• try to optimize the computation time of our algorithm because I think it's quite inefficient at the moment

[andreadelprete]

Balancing for matrix exp - Case study.pdf
Hi @olimexsmart, here you can find the updated analysis for the 2x2 matrix case, comparing our algorithm with the classic one. I've found a case study in which we gain a bit more (w.r.t. the one I told you yesterday) in terms of relative norm (wrt the classic balancing): the classic algorithm gets a norm that's 25% bigger than ours. But in the end we're still saving at most 1 squaring. This is the matrix:

A = [[1.5*M, M],
     [M, 0]]

It's the first case study, you can find it at the end of page 4 in the attached pdf.

The remaining open question is: "Could we gain more than 1 squarings for matrices bigger than 2x2?". I think we should try to figure that out by running tests with random matrices, rather than via mathematical analysis, because it seems too complex at a first glance.