RLS and the Kalman filter

Two adaptive filters learn the same 8-tap system from a strongly coloured input: the LMS of 26.3, and RLS. The curves show how far each one’s taps still are from the answer. Watch how early the RLS curve drops: after 50 samples it is at −40.7 dB, while LMS has barely moved.

RLS arrives in a few steps

In Adaptive filters: LMS (26.3), LMS learned the taps of an unknown system one small step at a time. Each step used the slope from a single sample. A coloured input made it crawl along the flat directions of the error bowl.

Think of a set of equations. You can solve them outright, or you can improve a guess a little at a time. LMS improves a guess. Recursive least squares, or RLS, solves outright at every new sample, and it does so cheaply.

This page is marked advanced because it works with matrices. From linear algebra I use only three things. They are a matrix times a vector, the transpose $^\top$ that turns a column into a row, and the inverse of a matrix, which undoes it.

Least squares over the past

The setting is 26.3’s system identification. An unknown FIR system is driven by an input $x[n]$, and its output plus a little white noise $v[n]$ is the desired signal $d_\text{des}[n]$. My filter has $N_h$ taps in a column $\mathbf{h}$. At sample $n$ it sees the tap-input vector $\mathbf{x}_n=(x[n],\dots,x[n-N_h+1])^\top$, newest first, as in 26.3.

In “Least squares or least worst” of Optimal FIR design (19.3), the least-squares design made the mean square of the error as small as possible. RLS does the same with the errors of every sample so far, added up. I count samples from $n=1$ on this page, and I give old errors less weight. After $n$ samples, the taps should make this sum as small as possible:

The number $\gamma$ is the forgetting factor, just below 1. It has no argument, so it is not 19.3’s weight $\gamma(\Omega)$. An error $k$ samples old counts $\gamma^k$ times as much as the newest one. With $\gamma=0.999$, an error 1000 samples old counts 0.368 times, so RLS remembers roughly the last $1/(1-\gamma)=1000$ samples.

Why forget at all? If the unknown system changes slowly, old samples describe a system that is gone. Forgetting lets RLS follow it.

In “The best FIR filter, from correlations” of The Wiener filter (26.2), the best taps solved the Wiener–Hopf equations $\mathbf{R}\mathbf{h}=\mathbf{r}$, built from expected values. Setting the slope of the sum above to zero along each tap gives the same equations, with weighted sums of the data in place of expected values:

The best taps after $n$ samples, $\mathbf{h}_n$, solve $\mathbf{R}_n\mathbf{h}_n=\mathbf{r}_n$.

Solving these afresh at every sample would cost on the order of $N_h^3$ operations each time. But the new matrix is the old one, scaled, plus one piece: $\mathbf{R}_n=\gamma\mathbf{R}_{n-1}+\mathbf{x}_n\mathbf{x}_n^\top$. So RLS keeps the inverse, $\mathbf{P}_n=\mathbf{R}_n^{-1}$, and updates it with that one piece. Here are the four lines it runs at each sample:

The first line is the gain $\mathbf{g}_n$, a column of $N_h$ numbers. The second is the error made with the old taps. The third moves the taps by the gain times that error. The last updates the inverse without inverting anything, and the linear algebra callout further on names the trick.

Before the first sample I start from $\mathbf{h}_0=\mathbf{0}$ and $\mathbf{P}_0=100\,\mathbf{I}$. That $\mathbf{P}_0$ is the inverse of $0.01\,\mathbf{I}$, a tiny made-up correlation that stands in until data arrive. It adds a small extra term to the sum, so the first solutions exist, and that term shrinks as $\gamma^n$. Apart from it, $\mathbf{h}_n$ is the exact least-squares answer at every $n$.

Now compare the steps. LMS moves its taps along $\mu\mathbf{x}_n$. RLS moves them along $\mathbf{g}_n$, which is $\mathbf{P}_{n-1}\mathbf{x}_n$ divided by a number.

Because $\mathbf{P}$ is the inverse of the correlation matrix, it stretches the step along the bowl’s flat directions by as much as the input’s colour squeezed them. So every direction is learned at one pace, whatever the colour.

The task

Let’s give both the same job. The unknown system has 8 taps, $h_\text{true}[k]=0.8^k\cos(0.4\pi k)$ for $k=0$ to 7. They are 1.000, 0.247, −0.518, −0.414, 0.127, 0.328, 0.081 and −0.170.

The input is AR(1) noise, as in Random processes (24.2): each value is 0.9 times the last plus a new white Gaussian draw. The draws have variance 0.19 (seed 264), and the first 300 outputs are dropped. That variance makes the input’s power $0.19/(1-0.9^2)=1$; this draw measures 0.990.

The noise $v[n]$ has standard deviation $\sigma_v=0.01$ (seed 2640) and measures RMS 0.0100.

This input is strongly coloured. Its correlation matrix $\mathbf{R}$ has the entries $0.9^{\lvert i-k\rvert}$, and its eigenvalues run from 0.055 to 6.20, a spread of 113. In 26.3’s bowl the spread was 4.0.

LMS gets the step size $\mu=0.02$, the $\mu$ of 26.3. Here $\operatorname{tr}\mathbf{R}$, the trace, adds up the diagonal, so it is 8. A common safe bound for LMS is $2/\operatorname{tr}\mathbf{R}=0.250$, below 26.3’s $2/\lambda_\text{max}=0.322$, and $\mu$ sits well under it. RLS gets $\gamma=0.999$.

To measure how far the taps are from the answer, I use the misalignment, in decibels:

Here $\lVert\mathbf{h}\rVert^2$ is the sum of the squares of the entries. At 0 dB the error is as large as the answer itself. At −30 dB its squared size is a thousandth of the answer’s.

The picture at the top of the page plots this misalignment for both filters. After 2000 samples LMS has crawled to −47.1 dB, and RLS sits at −56.1 dB, a floor set by the noise.

Notice how early RLS drops. After 8 samples it has 8 equations for its 8 unknowns, and it is at −33.3 dB. After 16 samples it is at −33.7 dB, while LMS is at −0.8 dB.

LMS is slow for the reason 26.3 gave. Each direction of the bowl shrinks by $1-\mu\lambda$ per step, so the flattest one, with $\lambda_\text{min}=0.055$, has the time constant $1/(\mu\lambda_\text{min})$, about 914 samples. LMS first gets below −30 dB at sample 917.

What RLS costs

RLS is not free. Count the multiplications in its four lines, treating a division as a multiplication. $\mathbf{P}_{n-1}\mathbf{x}_n$ takes $N_h^2$. $\mathbf{P}$ is symmetric, so $\mathbf{x}_n^\top\mathbf{P}_{n-1}$ is the same numbers as a row, and $\mathbf{x}_n^\top\mathbf{P}_{n-1}\mathbf{x}_n$ needs $N_h$ more.

The gain takes $N_h$ divisions, the error $N_h$ products and the taps $N_h$. The new $\mathbf{P}_n$ takes $N_h^2$ products and $N_h^2$ divisions by $\gamma$. That is $3N_h^2+4N_h$, which is 224 for 8 taps.

LMS needs $N_h$ for its output, 1 for $\mu e[n]$ and $N_h$ for the taps: $2N_h+1=17$. So RLS costs about 13 times as much per sample here, and the ratio grows with $N_h$.

Is the colour really the culprit? Give LMS a white input of the same nominal power (seed 2641, measured 1.017) and the same noise. It gets below −30 dB after 150 samples instead of 917, and RLS needs 10 there, against 8 with the coloured input. The colour, not the task, slows LMS.

The forgetting factor sets RLS’s floor. With a memory of about 1000 samples, the noise never quite averages away. Theory puts the floor at $\tfrac{1-\gamma}{2}\,\sigma_v^2\operatorname{tr}\mathbf{R}^{-1}/\lVert\mathbf{h}_\text{true}\rVert^2$, which is −56.9 dB here, and the clip ends at −56.1 dB. A smaller $\gamma$ follows a changing system sooner, but sits on a higher floor.

Predict, then correct

Now a different problem in the same spirit. A target moves in a plane, and every 0.1 s a sensor reports its position with an error of a few metres. I want its true path.

A ship’s navigator works this way. Between fixes, they work out the position from the last one, the speed and the heading; this is dead reckoning. Each new fix nudges that estimate: more if the fix is reliable, less if it is rough.

A model of the motion

In “What is stored, and what is pushed in” of Difference equations (6.1), the stored value was what a system carries from one sample to the next. For a moving target I carry the position and the velocity, and I call them the state. I treat the two axes separately. Along each axis the state is the column $\mathbf{s}_n=(\text{position},\text{velocity})^\top$ at time $nT_s$, with $T_s=0.1$ s.

In one step, the position moves on by $T_s$ times the velocity, and the velocity stays the same. In matrix form that is $\mathbf{A}\mathbf{s}_{n-1}$, with

Nothing steers the target, but a random acceleration jostles it, new at each step, with standard deviation $\sigma_a=0.3$ m/s². Over one step it adds $T_s^2/2$ times that acceleration to the position and $T_s$ times it to the velocity. This is the constant-velocity model: the velocity stays constant apart from those random pushes.

How big is a push? In “A cloud that leans: correlation” of Random variables for signals (24.1), the top of the correlation coefficient was $\mathbb{E}\{(x-\mu_x)(y-\mu_y)\}$. That average is the covariance of the two values.

Put the variances on the diagonal of a matrix and the covariances off it, and you have a covariance matrix. The push’s is

because each entry is a product of two of the push’s parts, $T_s^2/2$ and $T_s$, times the acceleration’s variance.

The sensor measures the position only, with white noise $v[n]$ of standard deviation $\sigma_v$:

The row $\mathbf{c}^\top$ picks the position out of the state. Along the x axis $\sigma_v$ is 3 m, and along the y axis 1.5 m.

This is a state-space model: one equation for how the state moves, and one for what the sensor sees. The page assumes the model is right, and that the push and the sensor noise are Gaussian with known sizes.

Predict, then update

The Kalman filter keeps an estimate of the state, which I also write $\mathbf{s}_n$, and the covariance matrix $\mathbf{P}_n$ of that estimate’s error: its own uncertainty. RLS’s $\mathbf{P}_n$ measured the same kind of thing, how unsure the taps were.

At each new time it does two things. Predict: run the model forward without a push, and grow the uncertainty by the push’s covariance. A minus sign up high marks a prediction:

Update: compare the measurement with the predicted position. The difference $y[n]-\mathbf{c}^\top\mathbf{s}_n^-$ is the innovation, the part of the measurement the model did not expect. Move the estimate by a fraction of it, set by the Kalman gain $\mathbf{g}_n$:

Why that gain? Write $P^-_{11}$ for the top-left entry of $\mathbf{P}_n^-$, the variance of the predicted position. The gain’s first entry, for the position, is $P^-_{11}/(P^-_{11}+\sigma_v^2)$.

If the prediction is much surer than the sensor, it is near 0, and the filter keeps its prediction. If the sensor is much surer, it is near 1, and the filter jumps to the measurement. In between, each is weighted by how sure it is.

The gain’s second entry moves the velocity, which no sensor measures: the filter learns it from how the positions drift.

Under the page’s assumptions, no other estimate made from the measurements so far has a smaller mean-square error. That is the sense in which R. E. Kalman’s filter, from his paper of 1960, is the optimal recursive estimator.

I start from the first measurement with zero velocity, $\mathbf{s}_0=(y[0],\ 0)^\top$, and from $\mathbf{P}_0=\mathrm{diag}(\sigma_v^2,100)$, where $\mathrm{diag}$ puts its numbers on the diagonal and zeros elsewhere. The 100 says I have little idea of the velocity: a standard deviation of 10 m/s.

The uncertainty ellipse

Along one axis, write $P_{11}$ for the top-left entry of $\mathbf{P}_n$. Then $\sqrt{P_{11}}$ is the standard deviation of the position error. With both axes, I draw an uncertainty ellipse around the estimate, reaching the x filter’s $2\sqrt{P_{11}}$ along x and the y filter’s along y. Its axes lie along x and y because the two filters run independently.

Picture 24.1’s cloud of Gaussian pairs with independent coordinates, and draw around it the ellipse two standard deviations out along each axis: it holds 86.5 % of the draws. So if the filter is honest, the true position should lie inside its ellipse about 86.5 % of the time.

Here is the run. The target starts at the origin with velocity (1.0, 0.5) m/s and follows the model above. Its random accelerations come from seed 2642, with $\sigma_a=0.3$ m/s²; their unit draws measure standard deviation 1.024.

The sensor’s noise is 3 m along x and 1.5 m along y (seed 2643), and measures 3.095 m and 1.518 m. The filter knows the model and those sizes, and runs once per axis for 30 s.

Each step plays in two halves, predict then correct: the ellipse grows, then shrinks as the estimate moves toward the new dot. By 29.9 s the estimate’s RMS error is 0.96 m, against the measurements’ 3.45 m.

Watch the ellipse shrink and then hold. It settles at ±1.25 m along x and ±0.74 m along y. Two measurement standard deviations are ±6 m and ±3 m, so the estimate is about four to five times surer than any single dot.

From 5 s on, once the start has been forgotten, the true position lay inside the ellipse 84.4 % of the time: 211 of 250 times. That is close to the Gaussian’s 86.5 %. The ellipse means what it says because the model is right. Over the same span the dots miss by RMS 3.45 m, and the estimate by 0.88 m.

The gains settle too: 0.044 along x and 0.061 along y. So each new dot moves the predicted position 4.4 % of the way toward it along x, and 6.1 % along y. The y sensor is twice as precise, so the filter believes it more.

RLS is a Kalman filter

Look back at RLS with this in mind. Let the state be the taps, with the model “the taps do not move”: $\mathbf{A}=\mathbf{I}$ and no push.

Let the measurement be $d_\text{des}[n]$, with the row $\mathbf{x}_n^\top$ in place of $\mathbf{c}^\top$ and a measurement-noise variance of 1. Let each prediction divide $\mathbf{P}$ by $\gamma$. Then the Kalman filter’s update is RLS’s four lines.

The forgetting factor plays the part of the push: it makes the filter a little less sure at every step, so the taps can still move.

Worked example

1. One Kalman step by hand. Take the tracker’s x axis, so $\sigma_v=3$ m and $\mathbf{P}_0=\mathrm{diag}(9,100)$. The first two measurements are $y[0]=3.86$ m and $y[1]=-2.79$ m; the true positions are 0 and 0.10 m. So $\mathbf{s}_0=(3.86,\ 0)^\top$.

Predict, ignoring $\mathbf{Q}$, which is tiny. The position moves on by $0.1\times0$, so $\mathbf{s}_1^-=(3.86,\ 0)^\top$. The top-left entry of $\mathbf{A}\mathbf{P}_0\mathbf{A}^\top$ is $9+2\cdot0.1\cdot0+0.01\cdot100=10$, and the entry below it is $0.1\cdot100=10$. With $\mathbf{Q}$ the first is 10.00000225.

The gain is $(10,\ 10)^\top/(10+9)$, so both entries are 0.526: the new estimate moves about half of the way to the new measurement. The innovation is $-2.79-3.86=-6.65$ m. The position becomes $3.86+0.526\times(-6.65)=0.36$ m, and the velocity $0+0.526\times(-6.65)=-3.50$ m/s.

The position variance falls to $10-0.526\times10=4.74$, so the ellipse reaches $2\sqrt{4.74}=4.35$ m along x, down from 6.00 m. The velocity is far off, since the truth is about 1.0 m/s. Its variance is still 94.7, so the filter knows it is unsure, and the next measurements fix it.

2. RLS against LMS, counted in work. RLS gets below −30 dB after 8 samples at 224 multiplications each, $8\times224=1792$ in all, and LMS after 917 samples at 17 each, $917\times17=15\,589$. So even in total work, RLS gets there about 8.7 times cheaper here. With the white input, LMS needs $150\times17=2550$, still more than RLS.

Where you’ll meet this

GPS receivers and inertial navigation blend a motion model with noisy fixes in a Kalman filter. Radar trackers do the same for each target they follow, as in Radar signal processing (33.1). Robots use Kalman filters to estimate where they are, and economists use state-space models for quantities they cannot measure directly. A Kalman filter helped navigate the Apollo spacecraft to the Moon.

Real targets turn and brake, and many sensors are not linear in the state. The extended and unscented Kalman filters handle such nonlinear models; this page only names them.

RLS fits where the filter is short and must learn fast, such as some adaptive equalisers in Channels and equalisation (33.2). For long filters, such as echo cancellers with hundreds of taps, its $N_h^2$ cost is why LMS and NLMS remain the usual choice.

Reference card