Least Squares and Filtering

We all know least squares fit given a set of measurement data \(x_1^*, x_2^*, \ldots, x_N^*\) which are actually sampling the actual values \(x_1, x_2, \ldots, x_N\) with some noise. Usually this noise is assumed to be Gaussian with zero mean and some variance.

If I have no knowledge of the underlying model, I might simply choose a zeroth order polynomial to fit the data, say \(P(x) = a_0\) where \(a_0\) is a constant. Our task is to find the best estimate of \(a_0\). I will denote with hat all the estimates. So \(\hat{x}_k\) is the best estimate of \(a_0\) given \(k\) observations of \(x\).

Clearly, given \(k\) data points, my best estimate is going to be

\[ \hat{x}_k = \frac{1}{k} \sum_{i=1}^k x_i^* \]

So lets simulate some data and see how it looks. We will randomly generate \(10\) data points using the actual model \(x = 2t + 1 + \epsilon\) where \(\epsilon \sim \mathcal{N}(0, 0.5)\). We will then fit a zeroth order polynomial to the data and see how well it fits the data. To quantify the goodness of fit, we will also calculate the mean squared error between the fitted polynomial and the data.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
t = np.linspace(0, 1, 10)

# Generate data
x = 2*t + 1 + np.random.normal(0, 0.1, 10)

# Fit a zeroth order polynomial
p0 = np.mean(x)

# Plot the data and the fitted polynomial
plt.scatter(t, x)
plt.plot(t, np.ones(10)*p0, label='Zeroth order polynomial')
plt.xlabel('t')
plt.ylabel('x')
plt.legend()
plt.show()

# Print the fitted polynomial
print('Fitted zeroth order polynomial:', p0)

# Calculate the mean squared error between the fitted polynomial and the data
mse_p0 = np.mean((np.ones(10)*p0 - x)**2)
print('Mean squared error between the zeroth order polynomial and the data:', mse_p0)

Fitted zeroth order polynomial: 2.073802317072883
Mean squared error between the zeroth order polynomial and the data: 0.3543766223136096

Now clearly this is a very bad fit. But we can blame only ourselves. So lets try something different.

Since,

\[ \hat{x}_k = \frac{1}{k} \sum_{i=1}^k x_i^* \]

therefore

\[ \begin{align*} \hat{x}_{k+1} & = \frac{1}{k+1} \sum_{i=1}^{k+1} x_i^*\\ & = \frac{1}{k+1} \left( x_{k+1}^* + \sum_{i=1}^{k} x_i^*\right)\\ & = \frac{1}{k+1} \left( x_{k+1}^* + k \hat{x}_k \right)\\ & = \frac{1}{k+1} \left( x_{k+1}^* + k \hat{x}_k + \hat{x}_k - \hat{x}_k \right)\\ & = \hat{x}_k + \frac{1}{k+1} \left( x_{k+1}^* - \hat{x}_k \right) \\ \end{align*}\]

This is a recursive formula for updating the estimate of \(a_0\) given a new observation \(x_{k+1}^*\). This is called the recursive least squares algorithm.

I can rewrite the same formula as

\[ \hat{x}_{k} = \hat{x}_{k-1} + \frac{1}{k} \left( x_{k}^* - \hat{x}_{k-1} \right) \]

Now let us implement this algorithm for the same data and see how well it fits the data. But this time we will plot all the estimates of the zeroth order polynomial as we get more and more data.

# Fit a zeroth order polynomial using recursive least squares
p0 = np.zeros(10)
p0[0] = x[0]
for i in range(1, 10):
    p0[i] = p0[i-1] + (x[i] - p0[i-1])/(i+1)

# Plot the data and the fitted polynomial
plt.scatter(t, x)
for i in range(10):
    plt.plot(t, p0, label='Zeroth order polynomial')
plt.xlabel('t')
plt.ylabel('x')
plt.show()

# Print the last estimate of the zeroth order polynomial

print('Last estimate of the fitted zeroth order polynomial:', p0[-1])


# Calculate the mean squared error between the fitted polynomial and the data
mse_p0 = np.mean((p0 - x)**2)
print('Mean squared error between the zeroth order polynomial and the data:', mse_p0)

Last estimate of the fitted zeroth order polynomial: 2.073802317072883
Mean squared error between the zeroth order polynomial and the data: 0.30537121854051574

If you compare the final estimate of \(a_0\) using the recursive least squares algorithm with the estimate obtained using all the data points (regular least squares), you will see that they match exactly to the machine precision. This is to be expected since the recursive least squares algorithm is nothing but the algebraic reordering of the regular least squares algorithm.

This second approach is called filtering. We improve on the estimate as more and more data comes in. This is very useful in real time applications where we need to make decisions based on the incoming data.

Lastly, I note that in this case, the mean squared error between the fitted polynomial and the data is smaller for the recursive least squares algorithm compared to the regular least squares algorithm.

I do not know if this holds for all cases.

You can read more about this in the excellent book Zarchan (2009).

References

Zarchan, P. 2009. Fundamentals of Kalman Filtering: A Practical Approach. AIAA.