Recursive least squares (RLS) algorithm is used in adaptive filters to find the filter coefficients that relate to recursively producing the least squares (minimum of the sum of the absolute squared) of the error signal (difference between the desired and the actual signal). This is contrast to other algorithms that aim to reduce the mean square error. The difference is that RLS filters are dependent on the signals themselves, whereas MSE filters are dependent on their statistics (specifically, the autocorrelation of the input and the cross-correlation of the input and desired signals). If these statistics are known, an MSE filter with fixed co-efficients (i.e., independent of the incoming data) can be built.
Contents |
Motivation
Suppose that a signal <math>\mathbf{d}</math> is transmitted over a echoey, noisy channel that causes it to be received as
- <math>x(n)=\sum_{k=0}^q b_n(k) d(n-k)+v(n)</math>
where <math>v(n)</math> represents white noise. We will attempt to recover the desired signal <math>\mathbf{d}</math> by use of an FIR filter, <math>\mathbf{w}</math>:
- <math>\hat{d}(n) = y(n) = \mathbf{w}_n^\mathit{T} \mathbf{x}(n)</math>
Our goal is to estimate the parameters of the filter <math>\mathbf{w}</math>, and at each time n we refer to the new least squares estimate by <math>\mathbf{w_n}</math>. As time evolves, we would like to avoid completely redoing the least squares algorithm to find the new estimate for <math>\mathbf{w}_{n+1}</math>, in terms of <math>\mathbf{w}_n</math>. The benefit of the RLS algorithm is that there is no need to invert matrices, thereby saving computational power. Another advantage is that it provides intuition behind such results as the Kalman filter.
Discussion
The idea behind RLS filters is to minimize a weighted least squares error function. To stay within the adaptive filter terminology, the weighted least squares error function is the cost function defined as
- <math>C(\mathbf{\mathbf{w}_n})=\sum_{i=0}^{n}\lambda^{n-i}|e(i)|^{2}=\sum_{i=0}^{n}\lambda^{n-i}e(i)\,e^{*}(i)</math>
where <math>0<\lambda\le 1</math> is an exponential weighting factor which effectively limits the number of input samples based on which the cost function is minimized. The error signal <math>e(n)</math> and desired signal <math>d(n)</math> are defined in the block diagram below: 
The cost function is minimized by taking the partial derivatives for all entries <math>k</math> of the coefficient vector <math>\mathbf{w}_{n}</math> and setting the results to zero
- <math>\frac{\partial C(\mathbf{w}_{n})}{\partial w^{*}_{n}(k)}=\sum_{i=0}^{n}\lambda^{n-i}e(i)\,\frac{\partial e^{*}(i)}{\partial w^{*}_{n}(k)}=\sum_{i=0}^{n}\lambda^{n-i}e(i)\,x^{*}(i-k)=0</math>
Next, replace <math>e(n)</math> with the definition of the error signal
- <math>\sum_{i=0}^{n}\lambda^{n-i}\left[d(i)-\sum_{l=0}^{p}w_{n}(l)x(i-l)\right]x^{*}(i-k)= 0</math>
Rearranging the equation yields
- <math>\sum_{l=0}^{p}w_{n}(l)\left[\sum_{i=0}^{n}\lambda^{n-i}\,x(i-l)x^{*}(i-k)\right]= \sum_{i=0}^{n}\lambda^{n-i}d(i)x^{*}(i-k)</math>
This form can be expressed in terms of matrices
- <math>\mathbf{R}_{x}(n)\,\mathbf{w}_{n}=\mathbf{r}_{dx}(n)</math>
where <math>\mathbf{R}_{x}(n)</math> is the weighted autocorrelation matrix for <math>x(n)</math> and <math>\mathbf{r}_{dx}(n)</math> is the cross-correlation between <math>d(n)</math> and <math>x(n)</math>. Based on this expression we find the coefficients which minimize the cost function as
- <math>\mathbf{w}_{n}=\mathbf{R}_{x}^{-1}(n)\,\mathbf{r}_{dx}(n)</math>
This is the main result of the discussion.
Choosing <math>\lambda</math>
The smaller <math>\lambda</math> is, the smaller contribution of previous samples. This makes the filter more sensitive the recent samples, which means more fluctuations in the filter co-efficients. The <math>\lambda=1</math> case is referred to as the growing window RLS algorithm.
Recursive algorithm
The discussion resulted in a single equation to determine a coefficient vector which minimizes the cost function. In this section we want to derive a recursive solution of the form
- <math>\mathbf{w}_{n}=\mathbf{w}_{n-1}+\Delta\mathbf{w}_{n-1}</math>
where <math>\Delta\mathbf{w}_{n-1}</math> is a correction factor at time n-1. We start the derivation of the recursive algorithm by expressing the cross correlation <math>\mathbf{r}_{dx}(n)</math> in terms of <math>\mathbf{r}_{dx}(n-1)</math>
-
<math>\mathbf{r}_{dx}(n)</math> <math>=\sum_{i=0}^{n}\lambda^{n-i}d(i)\mathbf{x}^{*}(i)</math> <math>=\sum_{i=0}^{n-1}\lambda^{n-i}d(i)\mathbf{x}^{*}(i)+\lambda^{0}d(n)\mathbf{x}^{*}(n)</math> <math>=\lambda\mathbf{r}_{dx}(n-1)+d(n)\mathbf{x}^{*}(n)</math>
where <math>\mathbf{x}(i)</math> is the p+1 dimensional data vector
- <math>\mathbf{x}(i)=[x(i), x(i-1), \dots , x(i-p) ]^{T}</math>
Similarly we express <math>\mathbf{R}_{x}(n)</math> in terms of <math>\mathbf{R}_{x}(n-1)</math> by
-
<math>\mathbf{R}_{x}(n)</math> <math>=\sum_{i=0}^{n}\lambda^{n-i}\mathbf{x}^{*}(i)\mathbf{x}^{T}(i)</math> <math>=\lambda\mathbf{R}_{x}(n-1)+\mathbf{x}^{*}(n)\mathbf{x}^{T}(n)</math>
In order to generate the coefficient vector we are interested in the inverse of the deterministic autocorrelation matrix. For that task the Woodbury matrix identity comes in handy. With
-
<math>A</math> <math>=\lambda\mathbf{R}_{x}(n-1)</math> is (p+1)-by-(p+1) <math>U</math> <math>=\mathbf{x}^{*}(n)</math> is (p+1)-by-1 <math>V</math> <math>=\mathbf{x}^{T}(n)</math> is 1-by-(p+1) <math>C</math> <math>=I_1</math> is the 1-by-1 identity matrix
The Woodbury matrix identity follows
-
<math>\mathbf{R}_{x}^{-1}(n)</math> <math>=</math> <math>\left[\lambda\mathbf{R}_{x}(n-1)+\mathbf{x}^{*}(n)\mathbf{x}^{T}(n)\right]^{-1}</math> <math>=</math> <math>\lambda^{-1}\mathbf{R}_{x}^{-1}(n-1)</math> <math>-\lambda^{-1}\mathbf{R}_{x}^{-1}(n-1)\mathbf{x}^{*}(n)</math> <math>\left\{1+\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{R}_{x}^{-1}(n-1)\mathbf{x}^{*}(n)\right\}^{-1} \mathbf{x}^{T}(n)\lambda^{-1}\mathbf{R}_{x}^{-1}(n-1)</math>
To come in line with the standard literature, we define
-
<math>\mathbf{P}(n)</math> <math>=\mathbf{R}_{x}^{-1}(n)</math> <math>=\lambda^{-1}\mathbf{P}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)</math>
where the gain vector <math>g(n)</math> is
-
<math>\mathbf{g}(n)</math> <math>=\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)\left\{1+\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)\right\}^{-1}</math> <math>=\mathbf{P}(n-1)\mathbf{x}^{*}(n)\left\{\lambda+\mathbf{x}^{T}(n)\mathbf{P}(n-1)\mathbf{x}^{*}(n)\right\}^{-1}</math>
Before we move on, it is necessary to bring <math>\mathbf{g}(n)</math> into another form
-
<math>\mathbf{g}(n)\left\{1+\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)\right\}</math> <math>=\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)</math> <math>\mathbf{g}(n)+\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)</math> <math>=\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)</math>
Subtracting the second term on the left side yields
-
<math>\mathbf{g}(n)</math> <math>=\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)-\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}^{*}(n)</math> <math>=\lambda^{-1}\left[\mathbf{P}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\mathbf{P}(n-1)\right]\mathbf{x}^{*}(n)</math>
With the recursive definition of <math>\mathbf{P}(n)</math> the desired form follows
- <math>\mathbf{g}(n)=\mathbf{P}(n)\mathbf{x}^{*}(n)</math>
Now we are ready to complete the recursion. As discussed
-
<math>\mathbf{w}_{n}</math> <math>=\mathbf{P}(n)\,\mathbf{r}_{dx}(n)</math> <math>=\lambda\mathbf{P}(n)\,\mathbf{r}_{dx}(n-1)+d(n)\mathbf{P}(n)\,\mathbf{x}^{*}(n)</math>
The second step follows from the recursive definition of <math>\mathbf{r}_{dx}(n )</math>. Next we incorporate the recursive definition of <math>\mathbf{P}(n)</math> together with the alternate form of <math>\mathbf{g}(n)</math> and get
-
<math>\mathbf{w}_{n}</math> <math>=\lambda\left[\lambda^{-1}\mathbf{P}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)\right]\mathbf{r}_{dx}(n-1)+d(n)\mathbf{g}(n)</math> <math>=\mathbf{P}(n-1)\mathbf{r}_{dx}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\mathbf{P}(n-1)\mathbf{r}_{dx}(n-1)+d(n)\mathbf{g}(n)</math> <math>=\mathbf{P}(n-1)\mathbf{r}_{dx}(n-1)+\mathbf{g}(n)\left[d(n)-\mathbf{x}^{T}(n)\mathbf{P}(n-1)\mathbf{r}_{dx}(n-1)\right]</math>
With <math>\mathbf{w}_{n-1}=\mathbf{P}(n-1)\mathbf{r}_{dx}(n-1)</math> we arrive at the update equation
-
<math>\mathbf{w}_{n}</math> <math>=\mathbf{w}_{n-1}+\mathbf{g}(n)\left[d(n)-\mathbf{x}^{T}(n)\mathbf{w}_{n-1}\right]</math> <math>=\mathbf{w}_{n-1}+\mathbf{g}(n)\alpha(n)</math>
where <math>\alpha(n)=d(n)-\mathbf{x}^{T}(n)\mathbf{w}_{n-1}</math> is the a priori error. Compare this with the a posteriori error; the error calculated after the filter is updated:
- <math>e(n)=d(n)-\mathbf{x}^{T}(n)\mathbf{w}_n</math>
That means we found the correction factor
- <math>\Delta\mathbf{w}_{n-1}=\mathbf{g}(n)\alpha(n)</math>
This intuitively satisfying result indicates that the correction factor is directly proportional to both the error and the gain vector, which controls how much sensitivity is desired, through the weighting factor, <math>\lambda</math>.
RLS algorithm summary
The RLS algorithm for a p-th order algorithm can be summarized as
| Parameters: | <math>p=</math> filter order |
| <math>\lambda=</math> forgetting factor | |
| <math>\delta=</math> value to initialize <math>\mathbf{P}(0)</math> | |
| Initialisation: | <math>\mathbf{w}_{n}=0</math> |
| <math>\mathbf{P}(0)=\delta^{-1}I</math> where <math>I</math> is the <math>(p+1)</math>-by-<math>(p+1)</math> identity matrix | |
| Computation: | For <math>n=0,1,2,\dots </math> |
|
<math> \mathbf{x}(n) = \left[ \begin{matrix} x(n)\\ x(n-1)\\ \vdots\\ x(n-p) \end{matrix} \right] </math> |
|
| <math> \alpha(n) = d(n)-\mathbf{w}_{n-1}^{T}\mathbf{x}(n)</math> | |
| <math>\mathbf{g}(n)=\mathbf{P}(n-1)\mathbf{x}^{*}(n)\left\{\lambda+\mathbf{x}^{T}(n)\mathbf{P}(n-1)\mathbf{x}^{*}(n)\right\}^{-1}</math> | |
| <math>\mathbf{P}(n)=\lambda^{-1}\mathbf{P}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)</math> | |
| <math> \mathbf{w}_{n} = \mathbf{w}_{n-1}+\,\alpha(n)\mathbf{g}(n)</math>. |
Note that the recursion for <math>P</math> follows a Riccati equation and thus draws parallels to the Kalman filter.
See also
References
- Hayes, Monson H. (1996). "9.4: Recursive Least Squares", Statistical Digital Signal Processing and Modeling. Wiley, pg.541. ISBN 0-471-59431-8.
- Simon Haykin, Adaptive Filter Theory, Prentice Hall, 2002, ISBN 0-13-048434-2
- M.H.A Davis, R.B. Vinter, Stochastic Modelling and Control, Springer, 1985, ISBN 0-412-16200-8
