BUG: linalg.expm: converting NaN to int

[andreas-schwab]

pick_pade_structure tries to convert NaN to int which is a big no-no and
results in undefined behaviour

[ilayn]

where does it do that?

[andreas-schwab]

lm = (int)ceil(log2(temp/normA/coeff[0])/6);

[andreas-schwab]

Test case from python-control:

from numpy import array, zeros, linspace
from control.matlab import ss, step
A = array([[-81.82, -45.45], [ 10., -1. ]])
B = array([[9.09], [0. ]])
C = array([[0, 0.159]])
D = zeros((1, 1))
sys_ss = ss(A, B, C, D)
t = linspace(0, 1000, 1000)
y, _t = step(sys_ss, t)

[ilayn]

Is this zero order hold discretized? What is the array that goes into expm? Because this is an overflow issue and should be handled first, not a NaN conversion problem.

[andreas-schwab]

normA is 0 and log2 is called with NaN

[ilayn]

The array is nonzero so that should not happen and if it is 0, that should return -inf. Let me check control library first how this model is discretized. We are treating the symptom we need the array that goes into expm.

[ilayn]

This model does go into the expm route because U is nonzero . What do you get out of this example? Could you please give a full account of what the issue is?

So something is happening here

https://github.com/python-control/python-control/blob/220ace34c8568cd14d71b9d2260336ba78235565/control/timeresp.py#L1159-L1191

[andreas-schwab]

(Pdb) p A
array([[-9.09111111e+03, -5.05000000e+03, 1.01000000e+03,
0.00000000e+00],
[ 1.11111111e+03, -1.11111111e+02, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00]])

[ilayn]

If that is the array then on main branch I get

>>> la.expm(np.array([[-9.09111111e+03, -5.05000000e+03, 1.01000000e+03,0.00000000e+00],
            [ 1.11111111e+03, -1.11111111e+02, 0.00000000e+00,0.00000000e+00],
            [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,1.00000000e+00],
            [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,0.00000000e+00]])
array([[0.        , 0.        , 0.01694884, 0.01707782],
       [0.        , 0.        , 0.16948837, 0.16925281],
       [0.        , 0.        , 1.        , 1.        ],
       [0.        , 0.        , 0.        , 1.        ]])

[andreas-schwab]

pick_pade_structure(Am) is returning 13, -2147483638

[ilayn]

This is what I get on main

In [8]: from scipy.linalg._matfuncs_expm import pick_pade_structure

In [9]: Am = np.empty((5, 4, 4), dtype=float)

In [10]: Am[0] = np.array([[-9.09111111e+03, -5.05000000e+03, 1.01000000e+03,
    ...:    ...: 0.00000000e+00],
    ...:    ...: [ 1.11111111e+03, -1.11111111e+02, 0.00000000e+00,
    ...:    ...: 0.00000000e+00],
    ...:    ...: [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
    ...:    ...: 1.00000000e+00],
    ...:    ...: [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
    ...:    ...: 0.00000000e+00]])

In [11]: pick_pade_structure(Am)
Out[11]: (13, 11)

Could you please check if you are on latest commit on main?

[andreas-schwab]

11 is the result of UNDEFINED BEHAVIOUR. You cannot trust it.

[ilayn]

Again, we are trying to find the root cause not the solution, I know where to put the guards if there is an overflow/undeflow. That is not the issue. We need to find where things go wrong. You are still fixing the resulting issue. Something is causing an overflow/underflow. So let's stop discussing about undefined behavior and find why it is causing a norm 0 on your machine. And it should return -inf not nan in c math.h so let's skip that part.

[andreas-schwab]

lm = ceil(log2(temp/normA/coeff[0])/6)

temp/normA is 0/0 is NaN

[ilayn]

You array is nonzero so normA can't be 0 because if these lines

scipy/scipy/linalg/_matfuncs_expm.c

Lines 616 to 618 in 860a46b

	normA = 0.0;
	for (i = 0; i < n; i++) { if (work_arr[n+i] > normA) { normA = work_arr[n+i]; } }

Next time it is modified is here

scipy/scipy/linalg/_matfuncs_expm.c

Line 753 in 860a46b

normA *= pow(2.0, -(*s));

pow result should be 0 I don't see how.

[andreas-schwab]

The debugger is always right.

[ilayn]

Then please use it properly and check the value of normA in those locations and report back what the values are. Because on three machines I have the result is 11. Hence, this seems like a won't fix problem so far.

[andreas-schwab]

11 is the result of UNDEFINED BEHAVIOUR. This the ULTIMATE SIN in programming!