Parameter space for Kalman is q/r rather than (q,r)

[pavelkomarov]

I've been playing with this script to plot the frequency response of a Kalman filter

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import dlti, freqz

dt = 1e-2
q = 1000
r = 0.01

I = np.eye(3)
A = np.array([[1, dt, (dt**2)/2], # states are x, x', x"
              [0, 1, dt],
              [0, 0,  1]])
C = np.array([[1, 0, 0]]) # we measure only y = noisy x
R = np.array([[r]])
Q = np.array([[0, 0, 0],
              [0, 0, 0],
              [0, 0, q]]) # uncertainty is around the acceleration
P0 = np.array(100*np.eye(3)) # See #110 for why this choice of P0

P = P0
for i in range(100):
	P_ = A @ P @ A.T + Q
	K = P_ @ C.T @ np.linalg.inv(C @ P_ @ C.T + R)
	P = (I - K @ C) @ P_

system = dlti((I - K @ C) @ A, K, C, np.array([[0]]), dt=dt) # Discrete LTI system, using closed-loop system matrices
transfer_fun = system.to_tf() # make transfer function

b = transfer_fun.num
a = transfer_fun.den

w, h = freqz(b, a, worN=1024)

plt.figure(figsize=(10, 3))
plt.plot(w / np.pi, 20 * np.log10(abs(h)))
plt.title(rf'Frequency response of constant acceleration Kalman Filter with $\Delta t = 0.01$, $q = 100$, $r=100$')
plt.xlabel("Frequency as a fraction of the Nyquist rate")
plt.ylabel("Magnitude (dB)")
plt.grid(True)
plt.tight_layout()
plt.show()

You get something sort of like this:

The results appear to be the same for certain ratios of $q:r$, irrespective of the scale of $q$ and $r$. This is odd, but I've been playing with filtering on synthetic examples just now and am seeing very similar results when I raise $q$ by an order of magnitude as when I lower $r$ by an order of magnitude:

versus:

I suspect, then, that the optimization need only consider a single parameter, and maybe we should only surface a qr_ratio parameter to users.

[pavelkomarov]

Some fun results from today: I worked out the frequency response of the RTS smoothing step. Here's the code

L = P @ A.T @ np.linalg.inv(A @ P @ A.T + Q)

system2 = dlti(L, (I - L @ A), I, np.zeros((3,3)), dt=dt)
transfer_fun2 = system2.to_tf()

b2 = transfer_fun2.num
a2 = transfer_fun2.den

w2, h2 = freqz(b2[0], a2, worN=1024)

When combined with the frequency response of the filter, you get this:

[pavelkomarov]

I think I see what's going on. If I increase $q$ and $r$ by an order of magnitude, the entries of the steady-state $P$ matrix also increase by an order of magnitude. Glad to see that has an effect. But all entries are evenly scaled up; its shape doesn't change. If I change the ratio, the shape of $P$ changes.

What is our best estimate in a Kalman filter? It's the central point from combining two Gaussian distributions with covariance $R$ and $C(APA^T + Q)C^T$ (which are in measurement space, technically, and we infer the central point back in state space). If I scale up $R = [r]$ and scale up $P$ and $Q$ by the same factor, then I don't change the shape of either of those measurement-space distributions, and the combination will have the same mean as before, although with a greater variance. But if I change the shape of one of the distributions (to be more oblong, for instance), the central point of the combination will move.

Since the filter's output in the plot above is coming from that mean point, and we're blind to the covariance around it, the estimate ends up looking the same for the same $q:r$ ratios, regardless of their scale. Mystery solved.

Should we now do anything about this in terms of searching parameter space? In theory if the values of $q$ and $r$ get too large or too small, you could end up with numerical issues, but the optimization limits them both to the range [1e-10, 1e10], so their ratio is in the range [1-20, 1e20]. I could see that starting to cause numerical issues if one of them is held constant at 1, so you'd probably want to set them by splitting the difference across the $0^\text{th}$ power.

I like being able to give $q$ and $r$ separately from a user-interface standpoint, though, because in modeling we usually make different process noise level and measurement noise level guesses, and we tend not to then go find their ratio. I suppose we could force users to do so and direct them to this issue as an explanation. Alternatively, we could give the Kalman optimization special treatment and find a way to break two parameters out from one. There's not presently a clean way to do this in the optimization code, and inventing one may not be worth it, since optimizing over the two parameters doesn't take very long or pose too much of a challenge as is.

[pavelkomarov]

One way of addressing this would be to go around the line:

cartesian_product = product(*[np.array(params[k]).astype(float) for k in search_space_types])

Instead of taking the cartesian product of {'forwardbackward': [True, False], 'q': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8], 'r': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8]}, we'd put a line in there to say "If the method is Kalman, actually do this simpler combination where you don't need quite as many elements." We can find the order of $q:r$ for $q \times r$ part of the Cartesian product as:

$$\begin{bmatrix} 0 & -4 & -7 & -9 & -12 & -16\\ 4 & 0 & -3 & -5 & -8 & -12\\ 7 & 3 & 0 & -2 & -5 & -9\\ 9 & 5 & 2 & 0 & -3 & -7\\ 12 & 8 & 5 & 3 & 0 & -4\\ 16 & 12 & 9 & 7 & 4 & 0\\ \end{bmatrix}$$

19 of the entries in this matrix are unique, with 17 repeats, so we're doing about a factor of 2 extra work, which doesn't seem terrible.

However, then within each Nelder-Mead, higher dimension does create longer running times as well, so the truly best way to address this would be to collapse the two into a single parameter and unpack before the parameters are passed to the func in _objective_function. That's a little more challenging, and the unpacking-from-dictionary code is already pretty complicated, so I'm not super keen to treat special cases.

My vote is this one doesn't really need to be addressed at the cost of increased code complexity, although I'd like it to be easy to address. I'll ponder. Worst case we can sit on this analysis for the future.