Recursive least squares filter

Discussion

The idea behind RLS filters is to minimize a cost function ${\displaystyle C}$ by appropriately selecting the filter coefficients ${\displaystyle {\mathbf{w}}_n}$, updating the filter as new data arrives. The error signal ${\displaystyle e(n)}$ and desired signal ${\displaystyle d(n)}$ are defined in the negative feedback diagram below:

The error implicitly depends on the filter coefficients through the estimate ${\displaystyle \hat{d}(n)}$:

${\displaystyle e(n)=d(n)-\hat{d}(n)}$

The weighted least squares error function ${\displaystyle C}$—the cost function we desire to minimize—being a function of ${\displaystyle e(n)}$ is therefore also dependent on the filter coefficients:

${\displaystyle C({\mathbf{w}}_n)=\sum _{i=0}^n{\lambda }^{n-i}{e}^2(i)}$

where ${\displaystyle 0<\lambda \le 1}$ is the "forgetting factor" which gives exponentially less weight to older error samples.

The cost function is minimized by taking the partial derivatives for all entries ${\displaystyle k}$ of the coefficient vector ${\displaystyle {\mathbf{w}}_n}$ and setting the results to zero

${\displaystyle \frac{\mathrm{\partial }C({\mathbf{w}}_n)}{\mathrm{\partial }{w}_n(k)}=\sum _{i=0}^{n}2{\lambda }^{n-i}e(i)\cdot \frac{\mathrm{\partial }e(i)}{\mathrm{\partial }{w}_n(k)}=-\sum _{i=0}^{n}2{\lambda }^{n-i}e(i)x(i-k)=0k=0,1,\dots ,p}$

Next, replace ${\displaystyle e(n)}$ with the definition of the error signal

${\displaystyle \sum _{i=0}^n{\lambda }^{n-i}\left[d(i)-\sum _{\ell =0}^p{w}_n(\ell )x(i-\ell )\right]x(i-k)=0k=0,1,\dots ,p}$

Rearranging the equation yields

${\displaystyle \sum _{\ell =0}^p{w}_n(\ell )\left[\sum _{i=0}^n{\lambda }^{n-i}x(i-\ell )x(i-k)\right]=\sum _{i=0}^n{\lambda }^{n-i}d(i)x(i-k)k=0,1,\dots ,p}$

This form can be expressed in terms of matrices

${\displaystyle {\mathbf{R}}_x(n){\mathbf{w}}_n={\mathbf{r}}_{dx}(n)}$

where ${\displaystyle {\mathbf{R}}_x(n)}$ is the weighted sample covariance matrix for ${\displaystyle x(n)}$, and ${\displaystyle {\mathbf{r}}_{dx}(n)}$ is the equivalent estimate for the cross-covariance between ${\displaystyle d(n)}$ and ${\displaystyle x(n)}$. Based on this expression we find the coefficients which minimize the cost function as

${\displaystyle {\mathbf{w}}_n={\mathbf{R}}_x^{-1}(n){\mathbf{r}}_{dx}(n)}$

This is the main result of the discussion.

Choosing λ

The smaller ${\displaystyle \lambda }$ is, the smaller is the contribution of previous samples to the covariance matrix. This makes the filter more sensitive to recent samples, which means more fluctuations in the filter co-efficients. The ${\displaystyle \lambda =1}$ case is referred to as the growing window RLS algorithm. In practice, ${\displaystyle \lambda }$ is usually chosen between 0.98 and 1. By using type-II maximum likelihood estimation the optimal ${\displaystyle \lambda }$ can be estimated from a set of data.

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

${\displaystyle {\mathbf{w}}_n={\mathbf{w}}_{n-1}+\mathrm{\Delta }{\mathbf{w}}_{n-1}}$

where ${\displaystyle \mathrm{\Delta }{\mathbf{w}}_{n-1}}$ is a correction factor at time ${\displaystyle n-1}$. We start the derivation of the recursive algorithm by expressing the cross covariance ${\displaystyle {\mathbf{r}}_{dx}(n)}$ in terms of ${\displaystyle {\mathbf{r}}_{dx}(n-1)}$

where ${\displaystyle \mathbf{x}(i)}$ is the ${\displaystyle p+1}$ dimensional data vector

${\displaystyle \mathbf{x}(i)=[x(i),x(i-1),\dots ,x(i-p){]}^T}$

Similarly we express ${\displaystyle {\mathbf{R}}_x(n)}$ in terms of ${\displaystyle {\mathbf{R}}_x(n-1)}$ by

In order to generate the coefficient vector we are interested in the inverse of the deterministic auto-covariance matrix. For that task the Woodbury matrix identity comes in handy. With

The Woodbury matrix identity follows

To come in line with the standard literature, we define

where the gain vector ${\displaystyle g(n)}$ is

Before we move on, it is necessary to bring ${\displaystyle \mathbf{g}(n)}$ into another form

Subtracting the second term on the left side yields

With the recursive definition of ${\displaystyle \mathbf{P}(n)}$ the desired form follows

${\displaystyle \mathbf{g}(n)=\mathbf{P}(n)\mathbf{x}(n)}$

Now we are ready to complete the recursion. As discussed

The second step follows from the recursive definition of ${\displaystyle {\mathbf{r}}_{dx}(n)}$. Next we incorporate the recursive definition of ${\displaystyle \mathbf{P}(n)}$ together with the alternate form of ${\displaystyle \mathbf{g}(n)}$ and get

With ${\displaystyle {\mathbf{w}}_{n-1}=\mathbf{P}(n-1){\mathbf{r}}_{dx}(n-1)}$ we arrive at the update equation

where ${\displaystyle \alpha (n)=d(n)-{\mathbf{x}}^T(n){\mathbf{w}}_{n-1}}$ is the a priori error. Compare this with the a posteriori error; the error calculated after the filter is updated:

${\displaystyle e(n)=d(n)-{\mathbf{x}}^T(n){\mathbf{w}}_n}$

That means we found the correction factor

${\displaystyle \mathrm{\Delta }{\mathbf{w}}_{n-1}=\mathbf{g}(n)\alpha (n)}$

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, ${\displaystyle \lambda }$.

RLS algorithm summary

The RLS algorithm for a p-th order RLS filter can be summarized as

The recursion for ${\displaystyle P}$ follows an algebraic Riccati equation and thus draws parallels to the Kalman filter.

Lattice recursive least squares filter (LRLS)

The lattice recursive least squares adaptive filter is related to the standard RLS except that it requires fewer arithmetic operations (order N). It offers additional advantages over conventional LMS algorithms such as faster convergence rates, modular structure, and insensitivity to variations in eigenvalue spread of the input correlation matrix. The LRLS algorithm described is based on a posteriori errors and includes the normalized form. The derivation is similar to the standard RLS algorithm and is based on the definition of ${\displaystyle d(k)}$. In the forward prediction case, we have ${\displaystyle d(k)=x(k)}$ with the input signal ${\displaystyle x(k-1)}$ as the most up to date sample. The backward prediction case is ${\displaystyle d(k)=x(k-i-1)}$, where i is the index of the sample in the past we want to predict, and the input signal ${\displaystyle x(k)}$ is the most recent sample.

Parameter summary

${\displaystyle {\kappa }_f(k,i)}$ is the forward reflection coefficient

${\displaystyle {\kappa }_b(k,i)}$ is the backward reflection coefficient

${\displaystyle e_f(k,i)}$ represents the instantaneous a posteriori forward prediction error

${\displaystyle e_b(k,i)}$ represents the instantaneous a posteriori backward prediction error

${\displaystyle {\xi }_{b_{min}}^d(k,i)}$ is the minimum least-squares backward prediction error

${\displaystyle {\xi }_{f_{min}}^d(k,i)}$ is the minimum least-squares forward prediction error

${\displaystyle \gamma (k,i)}$ is a conversion factor between a priori and a posteriori errors

${\displaystyle v_i(k)}$ are the feedforward multiplier coefficients.

${\displaystyle \epsilon }$ is a small positive constant that can be 0.01

LRLS algorithm summary

The algorithm for a LRLS filter can be summarized as

Normalized lattice recursive least squares filter (NLRLS)

The normalized form of the LRLS has fewer recursions and variables. It can be calculated by applying a normalization to the internal variables of the algorithm which will keep their magnitude bounded by one. This is generally not used in real-time applications because of the number of division and square-root operations which comes with a high computational load.

NLRLS algorithm summary

The algorithm for a NLRLS filter can be summarized as