{"id":1759,"date":"2024-05-24T23:30:00","date_gmt":"2024-05-25T03:30:00","guid":{"rendered":"http:\/\/aristotle2digital.blogwyrm.com\/?p=1759"},"modified":"2024-04-30T05:28:03","modified_gmt":"2024-04-30T09:28:03","slug":"time-series-6-recursive-least-squares","status":"publish","type":"post","link":"https:\/\/aristotle2digital.blogwyrm.com\/?p=1759","title":{"rendered":"Time Series 6 – Recursive Least Squares"},"content":{"rendered":"\n

The last post introduced the notion of the recursive least squares estimation approach as a natural precursor to the Kalman filter and derived the Kalman-like gain that governed the state estimation update equation.  Before moving on to a simple example of the Kalman filter, it seemed prudent to show the recursive least squares filter in action with a numerical example.<\/p>\n\n\n\n

The example involves a mass falling in an unknown uniform gravitational field with unknown initial conditions.  At regular time intervals, a measurement device measures the height of the mass above the ground but only to within certain precision due to inherent noise in it measuring processes.  The challenge is then to deliver the best estimate of the unknowns (initial height, initial speed, and the uniform gravitational acceleration).   We will assume that, beyond a measurement of the initial height, the other initial conditions are inaccessible to the measurement device since it only measures range.  This example is adapted from the excellent article Introduction to Kalman Filter: Derivation of the Recursive Least Squares Method<\/a><\/em> by Aleksander Haber.<\/p>\n\n\n\n

\"\"<\/figure><\/div>\n\n\n\n

To prove that the recursive least squares method does generate a real least squares estimation, let\u2019s first look at how a batch least squares approach to this problem would be performed.<\/p>\n\n\n\n

The dynamics of the falling mass are given by<\/p>\n\n\n\n

\\[ d(t) = d_0 + v_0 t + \\frac{1}{2} a t^2 \\; , \\]<\/p>\n\n\n\n

where $d(t)$ is the height above the ground at a given instant.<\/p>\n\n\n\n

We assume that the above equation is exact; in other words, there are no unknown forces in the problem.  In the language of Kalman-filtering and related optimal estimation, there is said to be no process noise, with it understood that the word \u2018process\u2019 is a genetic umbrella meaning the dynamics by which the state evolves in time, for which Newton\u2019s laws are an important subset.<\/p>\n\n\n\n

By design, measurements come in every $\\Delta t$ seconds and so the time evolution of the fall of the mass is sampled at times $t_k = k \\Delta t$.  At measurement time $t_k$ the amount the mass has dropped is<\/p>\n\n\n\n

\\[ x_k = x_0 + v_0 k \\Delta t + \\frac{1}{2} k^2 \\Delta t^2  a \\; \\]<\/p>\n\n\n\n

and the corresponding measurement is<\/p>\n\n\n\n

\\[ z_k = x_k + \\eta_k \\; , \\]<\/p>\n\n\n\n

where $\\eta_k$ is the noise the measuring device introduces in performing its job.<\/p>\n\n\n\n

\"\"<\/figure><\/div>\n\n\n\n

The batch least squares estimation begins by writing the evolution equation in a matrix form<\/p>\n\n\n\n

\\[ z_k = \\left[ \\begin{array}{ccc} 1 & k \\Delta t & \\frac{1}{2} k^2 \\Delta t^2 \\end{array} \\right] \\left[ \\begin{array}{c} x_0 \\\\ v_0 \\\\ a \\end{array} \\right]  + \\eta_k \\equiv H_{k,p} x_p + \\eta_k \\; .  \\]<\/p>\n\n\n\n

The array $x_p$ contains the constant initial conditions and unknown acceleration that is being estimated.<\/p>\n\n\n\n

For this numerical experiment, $x_0 = -4m$, $v_0 = 2 m\/s$ (positive being downward), and $a = 9.8 m\/s^2$ the usual value for surface of the Earth (give or take regional variations).  Measurements are taken every $\\Delta t = 1\/150 s$ over a time span of 14.2 seconds.  To illustrate the meaning of $H_{k,p}$, the value of the matrix at $k=0$ is<\/p>\n\n\n\n

\\[ H_{k=0,p} = \\left[ \\begin{array}{ccc} 1 & 0 & 0 \\end{array} \\right] \\; \\]<\/p>\n\n\n\n

corresponding to an elapsed time of $t_k = 0$ and at $k = 75$<\/p>\n\n\n\n

\\[ H_{k=0,p} = \\left[ \\begin{array}{ccc} 1 & 0.5 & 0.125 \\end{array} \\right] \\; \\]<\/p>\n\n\n\n

corresponding to an elapsed time of $t_k = 0.5$. <\/p>\n\n\n\n

The merit of this matrix approach is that by putting every measurement into a batch we can write the entire measured dynamics in a single equation as:<\/p>\n\n\n\n

\\[ z = H x + \\eta \\; . \\]<\/p>\n\n\n\n

The array $z$ has $M$ elements (for the number of measurement), where $M$ is the maximum value of the time index $k$ (for this specific example $M$ = 2131).  The matrix $H$ is a $M \\times n$ where $n$ is the number of estimate parameters ($n$ = 3, in this case) and $x$ is the array of parameters to be estimated).  The array $\\eta$ represents the $M$ values of random noise (assumed to be stationary) that were present in the measurements.  The values of $\\eta$ used in this numerical example were drawn from a normal distribution with a zero mean and a standard deviation of 1; the latter value determines the value of $R_k = 1$.<\/p>\n\n\n\n

Since, in the matrix equation above, the matrix $H$ isn\u2019t square, it may not be obvious how to solve for $x$.  However, there is a well-known technique using what is known as the Moore-Penrose pseudoinverse<\/a>. <\/p>\n\n\n\n

In short, to solve for $x$ in a least squares sense, we set $\\eta$ to zero (i.e., we ignore the noise assuming the classical argument that with enough measurements the influence of the noise will average to zero).  We then left-multiply by the transpose of $H$ to get<\/p>\n\n\n\n

\\[ H^T z = (H^T H) x \\; .\\]<\/p>\n\n\n\n

The quantity $H^T H$ is a $(n \\times K) \\cdot (K \\times n) = n \\times n $ matrix (in this case $3\\times 3$, which is invertible.  We then arrive at the compact form of the estimate for $x$ as<\/p>\n\n\n\n

\\[ x_{est} = (H^T H)^{-1} H^T z \\; .\\]<\/p>\n\n\n\n

For this numerical experiment, the resulting estimate is<\/p>\n\n\n\n

\\[ x_{est} = \\left[ \\begin{array}{c} -3.95187658 \\\\ 1.99300377 \\\\ 9.80040125 \u00a0\\end{array}\u00a0\u00a0\u00a0 \\right] \\; .\\]<\/p>\n\n\n\n

This is the same value one gets by using a package like numpy\u2019s linalg.lstsq.<\/p>\n\n\n\n

In the recursive method, all of the data are not stored and handled at once but rather each point is taken and processed as it comes in. <\/p>\n\n\n\n

The algorithm developed in the last post is summarized as:<\/p>\n\n\n\n