Max number of matrix-matrix multiplications in expm

[andreadelprete]

It could be useful to allow the user to specify a maximum number of matrix-matrix multiplications in the computation of the matrix exponential.

[olimexsmart]

In compute times vector it makes sense to have this limiting factor on top of the number of squarings determined as: mat_squarings = squarings - vec_squarings?

[andreadelprete]

https://github.com/andreadelprete/consim/blob/master/script/consim_py/utils/exponential_matrix_helper.py

[olimexsmart]

I've commited the changes, my tests pass but I still need to confirm the new functionalities. I thinking about confronting them with the python version.

All the TimeScaling stuff, like the second order LDS class, is still relevant?

[andreadelprete]

For the test, an option is to check that by decreasing the max number of mat. mult. you get a graceful degradation of the accuracy of the exponential matrix.

Regarding time scaling, we kind of replaced it with balancing because we've figured out that TS is just a special case of balancing. I don't know if it could still be useful. Maybe not.

[andreadelprete]

@olimexsmart as soon as this feature is tested, please notify @hammoudbilal so he can start using it in consim. Especially, it would be nice to be able to change the max number of mat. mult. from Python for running tests.

[olimexsmart]

I've run some tests, gradually incrementing the max number of multiplications. The first set of results are relative to a matrix with a "big" norm, 4 squarings required.

Correct result norm: 4.74453e+07
Max Multiplications: -1 error: 0
Max Multiplications: 0 error: 4.74453e+07
Max Multiplications: 1 error: 4.74453e+07
Max Multiplications: 2 error: 4.7445e+07
Max Multiplications: 3 error: 4.74454e+07
Max Multiplications: 4 error: 4.74462e+07
Max Multiplications: 5 error: 4.74802e+07
Max Multiplications: 6 error: 1.33422e+07
Max Multiplications: 7 error: 0.0227784
Max Multiplications: 8 error: 1.14632e-07
Max Multiplications: 9 error: 6.88169e-08
Max Multiplications: 10 error: 0
Max Multiplications: 11 error: 0
Max Multiplications: 12 error: 0

The second set for a smaller norm matrix, 2 squarings normally required:

Correct result norm: 306.622
Max Multiplications: -1 error: 0
Max Multiplications: 0 error: 315.838
Max Multiplications: 1 error: 290.501
Max Multiplications: 2 error: 539.902
Max Multiplications: 3 error: 3.55689
Max Multiplications: 4 error: 0.0046275
Max Multiplications: 5 error: 2.20303e-06
Max Multiplications: 6 error: 3.51989e-12
Max Multiplications: 7 error: 2.67098e-13
Max Multiplications: 8 error: 0
Max Multiplications: 9 error: 0

Does it makes sense? My doubts are about the saturation of the error to the norm of the correct result. Basically meaning that the results obtained with a small number of multiplications are just small in norm

[andreadelprete]

These results do not make much sense to me, but maybe it's because I've not understood how you computed the errors. In my mind, errors should be computed as the norm of the different between the "exact" expm (computed without any bound on the number of mat. mult.) and the "approximated" expm (computed with a specific bound on the number of mat. mult.). With this definition of the error, I expect it to grow monotonically as max-mult decreases.

[olimexsmart]

The errors were computed exactly how you imagine: error = (reference - result).norm()

The point I was trying to make is that if you have a matrix which exponential's norm is about e+07, forcing the pade aprox to 1, 2, 3, ecc produces results that have much smaller norms, like e+00, e+02.

This means that in the error the norm of the reference result dominates over the results with reduced multiplications. The error decreases exponentially, not linearly.

Is this right?

[andreadelprete]

Ok, now that I'm taking a second look, the results do make sense. The phenomena you're observing is simply due to the fact that you're testing unstable matrices, that is with positive eigenvalues. If you test stable real matrices (as we should already have in consim) the results will look different.

[olimexsmart]

Ok, I will look into a method to generate random matrices with given eigenvalues. In the mean time we could assume the implementation is correct?

@hammoudbilal the usage is very simple, I've added a get and set properties to change the max number of multiplications both in LDSUtilities.h and MatrixExponential.h:

    void setMaxMultiplications(int mm);
    int getMaxMultiplications();

[andreadelprete]

For generating a stable matrix you can simply do this:

U = random(n,n)
A = -U*U.T

[olimexsmart]

The formula works, it produces much bigger exponential results.

But I don't notice a change in the behavior of the error. The eigenvalues now are:

 (116.526,0)
 (75.5768,0)
 (40.7053,0)
 (23.1863,0)
 (17.5936,0)
 (13.4599,0)
 (4.18339,0)
 (6.45416,0)
(0.107136,0)
(0.241162,0)

And the error:

Correct result norm: 4.04277e+50
Max Multiplications: -1 error: 0
Max Multiplications: 0 error: 4.04277e+50
Max Multiplications: 1 error: 4.04277e+50
Max Multiplications: 2 error: 4.04277e+50
Max Multiplications: 3 error: 4.04277e+50
Max Multiplications: 4 error: 4.04277e+50
Max Multiplications: 5 error: 4.04277e+50
Max Multiplications: 6 error: 4.04277e+50
Max Multiplications: 7 error: 4.04277e+50
Max Multiplications: 8 error: 4.04277e+50
Max Multiplications: 9 error: 4.97511e+48
Max Multiplications: 10 error: 1.80335e+40
Max Multiplications: 11 error: 9.39615e+36
Max Multiplications: 12 error: 0
Max Multiplications: 13 error: 0
Max Multiplications: 14 error: 0

[andreadelprete]

The eigenvalues should be negative. Did you forget the minus?

[olimexsmart]

Did you forget the minus?

I forget to turn on the brain

(-1294.74,0)
(-839.742,0)
(-452.281,0)
(-257.626,0)
(-195.485,0)
(-149.554,0)
(-46.4821,0)
 (-1.1904,0)
(-71.7129,0)
(-2.67957,0)

Correct result norm: 2.224e-05
Max Multiplications: -1 error: 0
Max Multiplications: 0 error: 3.02292
Max Multiplications: 1 error: 2.88648
Max Multiplications: 2 error: 2.80452
Max Multiplications: 3 error: 2.70771
Max Multiplications: 4 error: 10.4774
Max Multiplications: 5 error: 6.15498
Max Multiplications: 6 error: 40.956
Max Multiplications: 7 error: 21.6585
Max Multiplications: 8 error: 74.4735
Max Multiplications: 9 error: 26.1126
Max Multiplications: 10 error: 0.0161914
Max Multiplications: 11 error: 4.26435e+16
Max Multiplications: 12 error: 4.37617e-07
Max Multiplications: 13 error: 7.49102e-11
Max Multiplications: 14 error: 2.94241e-13
Max Multiplications: 15 error: 1.61602e-15
Max Multiplications: 16 error: 8.08984e-17
Max Multiplications: 17 error: 5.41386e-18
Max Multiplications: 18 error: 0
Max Multiplications: 19 error: 0

Well this is strange, the decrease is there but is not monotonic. I've tried different values and it can have multiple peaks before a very big one and then a dramatic descend. That peak corresponds with the value before the normal number of squarings.

i.e. in this case the algorithm would normally perform 18 multiplications, 6 in pade13 and 12 in squarings

I have not looked into it yet

[andreadelprete]

This looks like a big fat bug! :)

[olimexsmart]

I've corrected the bug we have talked about, but that influenced only test matrices that have a very small norm, that normally would be computed with pade smaller that 13. None my test and results above were performed with matrices with such small norm.

My results does not change from my previous post.

What I notice now is that the norms increases stepping from pade5 (MaxM=3) to pade7 (MaxM = 4), for example. Again between pade9 and pade13.

That is something I have no control on, I just call the given pade approximations, which I haven't modified.

I am starting to think either the input matrices from your tests had different features, or the pade approximations and/or other given code is different between Eigen and Python

[andreadelprete]

Maybe you can try to reproduce this test in Python to see what results you get there. The test I did was only looking at the error of the simulator, so that was different.

[olimexsmart]

I've reproduced the same test, with the same input data, in Python here using the class ExponentialMatrixHelper.

And the method _solve_P_Q explodes at maxMul = 5 saying that the matrix is singular.

This means that the difference is both in the code and the input data. Eigen just solves, giving back results of ill-conditioned operations, which lead to bigger errors. While Numpy has some run-time checks and quits. At the same time the matrices generated by the simulator does not trigger the Numpy checks.

If you provide me a matrix produced by the simulator I could try re-running all the tests

[andreadelprete]

That's really weird. I've read that in theory the matrix Q should always be well-conditioned. Regarding the matrices computed by the simulator, didn't you already log some of them in the past?

[olimexsmart]

With the data from the Python simulator it works fine, I guess:

Correct result norm: 232.042
Max Multiplications: -1 error: 0
Max Multiplications: 0 error: 872.096
Max Multiplications: 1 error: 726.695
Max Multiplications: 2 error: 535.24
Max Multiplications: 3 error: 297.684
Max Multiplications: 4 error: 140.568
Max Multiplications: 5 error: 53.0464
Max Multiplications: 6 error: 3.88075
Max Multiplications: 7 error: 1.23547e-06
Max Multiplications: 8 error: 3.36535e-10
Max Multiplications: 9 error: 3.10329e-10
Max Multiplications: 10 error: 3.10199e-10
Max Multiplications: 11 error: 3.05411e-10
Max Multiplications: 12 error: 2.90114e-10
Max Multiplications: 13 error: 3.90878e-10
Max Multiplications: 14 error: 7.40558e-10
Max Multiplications: 15 error: 1.82847e-09
Max Multiplications: 16 error: 4.70624e-09
Max Multiplications: 17 error: 1.65604e-08
Max Multiplications: 18 error: 6.23341e-08
Max Multiplications: 19 error: 0
Max Multiplications: 20 error: 0
Max Multiplications: 21 error: 0

A little bit of ramp up at the end but no serious spikes

[andreadelprete]

That's ok, it's quite small. I think this test validates the code.