Recursive least squares filter

Short description: Adaptive filter algorithm for digital signal processing

Recursive least squares (RLS) is an adaptive filter algorithm that recursively finds the coefficients that minimize a weighted linear least squares cost function relating to the input signals. This approach is in contrast to other algorithms such as the least mean squares (LMS) that aim to reduce the mean square error. In the derivation of the RLS, the input signals are considered deterministic, while for the LMS and similar algorithms they are considered stochastic. Compared to most of its competitors, the RLS exhibits extremely fast convergence. However, this benefit comes at the cost of high computational complexity.

## Contents

* 1 Motivation
* 2 Discussion
* 2.1 Choosing λ
* 3 Recursive algorithm
* 4 RLS algorithm summary
* 5 Lattice recursive least squares filter (LRLS)
* 5.1 Parameter summary
* 5.2 LRLS algorithm summary
* 6 Normalized lattice recursive least squares filter (NLRLS)
* 6.1 NLRLS algorithm summary
* 7 See also
* 8 References
* 9 Notes

## Motivation

RLS was discovered by Gauss but lay unused or ignored until 1950 when Plackett rediscovered the original work of Gauss from 1821. In general, the RLS can be used to solve any problem that can be solved by adaptive filters. For example, suppose that a signal d(n) is transmitted over an echoey, noisy channel that causes it to be received as

x(n)=∑k=0qbn(k)d(n−k)+v(n)

where v(n) represents additive noise. The intent of the RLS filter is to recover the desired signal d(n) by use of a p+1-tap FIR filter, 𝐰:

d(n)≈∑k=0pw(k)x(n−k)=𝐰T𝐱n

where 𝐱n=[x(n)x(n−1)…x(n−p)]T is the column vector containing the p+1 most recent samples of x(n). The estimate of the recovered desired signal is

d^(n)=∑k=0pwn(k)x(n−k)=𝐰nT𝐱n

The goal is to estimate the parameters of the filter 𝐰, and at each time n we refer to the current estimate as 𝐰n and the adapted least-squares estimate by 𝐰n+1. 𝐰n is also a column vector, as shown below, and the transpose, 𝐰nT, is a row vector. The matrix product 𝐰nT𝐱n (which is the dot product of 𝐰n and 𝐱n) is d^(n), a scalar. The estimate is "good" if d^(n)−d(n) is small in magnitude in some least squares sense.

As time evolves, it is desired to avoid completely redoing the least squares algorithm to find the new estimate for 𝐰n+1, in terms of 𝐰n.

The benefit of the RLS algorithm is that there is no need to invert matrices, thereby saving computational cost. Another advantage is that it provides intuition behind such results as the Kalman filter.

## Discussion

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

500px

The error implicitly depends on the filter coefficients through the estimate d^(n):

e(n)=d(n)−d^(n)

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

C(𝐰n)=∑i=0nλn−ie2(i)

where 0<λ≤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 k of the coefficient vector 𝐰n and setting the results to zero

∂C(𝐰n)∂wn(k)=∑i=0n2λn−ie(i)⋅∂e(i)∂wn(k)=−∑i=0n2λn−ie(i)x(i−k)=0k=0,1,…,p

Next, replace e(n) with the definition of the error signal

∑i=0nλn−i[d(i)−∑ℓ=0pwn(ℓ)x(i−ℓ)]x(i−k)=0k=0,1,…,p

Rearranging the equation yields

∑ℓ=0pwn(ℓ)[∑i=0nλn−ix(i−ℓ)x(i−k)]=∑i=0nλn−id(i)x(i−k)k=0,1,…,p

This form can be expressed in terms of matrices

𝐑x(n)𝐰n=𝐫dx(n)

where 𝐑x(n) is the weighted sample covariance matrix for x(n), and 𝐫dx(n) is the equivalent estimate for the cross-covariance between d(n) and x(n). Based on this expression we find the coefficients which minimize the cost function as

𝐰n=𝐑x−1(n)𝐫dx(n)

This is the main result of the discussion.

### Choosing λ

The smaller λ 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 λ=1 case is referred to as the growing window RLS algorithm. In practice, λ is usually chosen between 0.98 and 1.[1] By using type-II maximum likelihood estimation the optimal λ can be estimated from a set of data.[2]

## 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

𝐰n=𝐰n−1+Δ𝐰n−1

where Δ𝐰n−1 is a correction factor at time n−1. We start the derivation of the recursive algorithm by expressing the cross covariance 𝐫dx(n) in terms of 𝐫dx(n−1)

𝐫dx(n) | =∑i=0nλn−id(i)𝐱(i)

|

| =∑i=0n−1λn−id(i)𝐱(i)+λ0d(n)𝐱(n)
| =λ𝐫dx(n−1)+d(n)𝐱(n)

where 𝐱(i) is the p+1 dimensional data vector

𝐱(i)=[x(i),x(i−1),…,x(i−p)]T

Similarly we express 𝐑x(n) in terms of 𝐑x(n−1) by

𝐑x(n) | =∑i=0nλn−i𝐱(i)𝐱T(i)

|

| =λ𝐑x(n−1)+𝐱(n)𝐱T(n)

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

A | =λ𝐑x(n−1) is (p+1)-by-(p+1)

|

U | =𝐱(n) is (p+1)-by-1 (column vector)
V | =𝐱T(n) is 1-by-(p+1) (row vector)
C | =𝐈1 is the 1-by-1 identity matrix

The Woodbury matrix identity follows

𝐑x−1(n) | = | [λ𝐑x(n−1)+𝐱(n)𝐱T(n)]−1

|
|

| = | 1λ{𝐑x−1(n−1)−𝐑x−1(n−1)𝐱(n)𝐱T(n)𝐑x−1(n−1)λ+𝐱T(n)𝐑x−1(n−1)𝐱(n)}
|  |

To come in line with the standard literature, we define

𝐏(n) | =𝐑x−1(n)

|

| =λ−1𝐏(n−1)−𝐠(n)𝐱T(n)λ−1𝐏(n−1)

where the gain vector g(n) is

𝐠(n) | =λ−1𝐏(n−1)𝐱(n){1+𝐱T(n)λ−1𝐏(n−1)𝐱(n)}−1

|

| =𝐏(n−1)𝐱(n){λ+𝐱T(n)𝐏(n−1)𝐱(n)}−1

Before we move on, it is necessary to bring 𝐠(n) into another form

𝐠(n){1+𝐱T(n)λ−1𝐏(n−1)𝐱(n)} | =λ−1𝐏(n−1)𝐱(n)

|

𝐠(n)+𝐠(n)𝐱T(n)λ−1𝐏(n−1)𝐱(n) | =λ−1𝐏(n−1)𝐱(n)

Subtracting the second term on the left side yields

𝐠(n) | =λ−1𝐏(n−1)𝐱(n)−𝐠(n)𝐱T(n)λ−1𝐏(n−1)𝐱(n)

|

| =λ−1[𝐏(n−1)−𝐠(n)𝐱T(n)𝐏(n−1)]𝐱(n)

With the recursive definition of 𝐏(n) the desired form follows

𝐠(n)=𝐏(n)𝐱(n)

Now we are ready to complete the recursion. As discussed

𝐰n | =𝐏(n)𝐫dx(n)

|

| =λ𝐏(n)𝐫dx(n−1)+d(n)𝐏(n)𝐱(n)

The second step follows from the recursive definition of 𝐫dx(n). Next we incorporate the recursive definition of 𝐏(n) together with the alternate form of 𝐠(n) and get

𝐰n | =λ[λ−1𝐏(n−1)−𝐠(n)𝐱T(n)λ−1𝐏(n−1)]𝐫dx(n−1)+d(n)𝐠(n)

|

| =𝐏(n−1)𝐫dx(n−1)−𝐠(n)𝐱T(n)𝐏(n−1)𝐫dx(n−1)+d(n)𝐠(n)
| =𝐏(n−1)𝐫dx(n−1)+𝐠(n)[d(n)−𝐱T(n)𝐏(n−1)𝐫dx(n−1)]

With 𝐰n−1=𝐏(n−1)𝐫dx(n−1) we arrive at the update equation

𝐰n | =𝐰n−1+𝐠(n)[d(n)−𝐱T(n)𝐰n−1]

|

| =𝐰n−1+𝐠(n)α(n)

where α(n)=d(n)−𝐱T(n)𝐰n−1 is the a priori error. Compare this with the a posteriori error; the error calculated after the filter is updated:

e(n)=d(n)−𝐱T(n)𝐰n

That means we found the correction factor

Δ𝐰n−1=𝐠(n)α(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, λ.

## RLS algorithm summary

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

Parameters: | p= filter order

|

| λ= forgetting factor
| δ= value to initialize 𝐏(0)
Initialization: | 𝐰(0)=0,
| x(k)=0,k=−p,…,−1,
| d(k)=0,k=−p,…,−1
| 𝐏(0)=δI where I is the identity matrix of rank p+1
Computation: | For n=1,2,…

𝐱(n)=[x(n)x(n−1)⋮x(n−p)]

| α(n)=d(n)−𝐱T(n)𝐰(n−1)
| 𝐠(n)=𝐏(n−1)𝐱(n){λ+𝐱T(n)𝐏(n−1)𝐱(n)}−1
| 𝐏(n)=λ−1𝐏(n−1)−𝐠(n)𝐱T(n)λ−1𝐏(n−1)
| 𝐰(n)=𝐰(n−1)+α(n)𝐠(n).

The recursion for P follows an algebraic Riccati equation and thus draws parallels to the Kalman filter.[3]

## 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).[4] 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 d(k). In the forward prediction case, we have d(k)=x(k) with the input signal x(k−1) as the most up to date sample. The backward prediction case is 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 x(k) is the most recent sample.[5]

### Parameter summary

κf(k,i) is the forward reflection coefficient

κb(k,i) is the backward reflection coefficient

ef(k,i) represents the instantaneous a posteriori forward prediction error

eb(k,i) represents the instantaneous a posteriori backward prediction error

ξbmind(k,i) is the minimum least-squares backward prediction error

ξfmind(k,i) is the minimum least-squares forward prediction error

γ(k,i) is a conversion factor between a priori and a posteriori errors

vi(k) are the feedforward multiplier coefficients.

ε is a small positive constant that can be 0.01

### LRLS algorithm summary

The algorithm for a LRLS filter can be summarized as

Initialization:

| For i=0,1,…,N
|  δ(−1,i)=δD(−1,i)=0 (if x(k)=0 for k<0)
|  ξbmind(−1,i)=ξfmind(−1,i)=ε
|  γ(−1,i)=1
|  eb(−1,i)=0
| End
Computation:
| For k≥0
|  γ(k,0)=1
|  eb(k,0)=ef(k,0)=x(k)
|  ξbmind(k,0)=ξfmind(k,0)=x2(k)+λξfmind(k−1,0)
|  e(k,0)=d(k)
|  For i=0,1,…,N
|  δ(k,i)=λδ(k−1,i)+eb(k−1,i)ef(k,i)γ(k−1,i)
|  γ(k,i+1)=γ(k,i)−eb2(k,i)ξbmind(k,i)
|  κb(k,i)=δ(k,i)ξfmind(k,i)
|  κf(k,i)=δ(k,i)ξbmind(k−1,i)
|  eb(k,i+1)=eb(k−1,i)−κb(k,i)ef(k,i)
|  ef(k,i+1)=ef(k,i)−κf(k,i)eb(k−1,i)
|  ξbmind(k,i+1)=ξbmind(k−1,i)−δ(k,i)κb(k,i)
|  ξfmind(k,i+1)=ξfmind(k,i)−δ(k,i)κf(k,i)
|  Feedforward filtering
|  δD(k,i)=λδD(k−1,i)+e(k,i)eb(k,i)γ(k,i)
|  vi(k)=δD(k,i)ξbmind(k,i)
|  e(k,i+1)=e(k,i)−vi(k)eb(k,i)
|  End
| End

## 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

Initialization:

| For i=0,1,…,N.
|  δ‾(−1,i)=0 (if x(k)=d(k)=0 for k<0)
|  δ‾D(−1,i)=0
|  e‾b(−1,i)=0
| End
|  σx2(−1)=λσd2(−1)=ε
Computation:
| For k≥0
|  σx2(k)=λσx2(k−1)+x2(k) (Input signal energy)
|  σd2(k)=λσd2(k−1)+d2(k) (Reference signal energy)
|  e‾b(k,0)=e‾f(k,0)=x(k)σx(k)
|  e‾(k,0)=d(k)σd(k)
|  For i=0,1,…,N
|  δ‾(k,i)=δ(k−1,i)(1−e‾b2(k−1,i))(1−e‾f2(k,i))+e‾b(k−1,i)e‾f(k,i)
|  e‾b(k,i+1)=e‾b(k−1,i)−δ‾(k,i)e‾f(k,i)(1−δ‾2(k,i))(1−e‾f2(k,i))
|  e‾f(k,i+1)=e‾f(k,i)−δ‾(k,i)e‾b(k−1,i)(1−δ‾2(k,i))(1−e‾b2(k−1,i))
|  Feedforward filter
|  δ‾D(k,i)=δ‾D(k−1,i)(1−e‾b2(k,i))(1−e‾2(k,i))+e‾(k,i)e‾b(k,i)
|  e‾(k,i+1)=1(1−e‾b2(k,i))(1−δ‾D2(k,i))[e‾(k,i)−δ‾D(k,i)e‾b(k,i)]
|  End
| End

## See also

* Adaptive filter
* Kernel adaptive filter
* Least mean squares filter
* Zero-forcing equalizer

## References

* Hayes, Monson H. (1996). "9.4: Recursive Least Squares". Statistical Digital Signal Processing and Modeling. Wiley. p. 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
* Weifeng Liu, Jose Principe and Simon Haykin, Kernel Adaptive Filtering: A Comprehensive Introduction, John Wiley, 2010, ISBN 0-470-44753-2
* R.L.Plackett, Some Theorems in Least Squares, Biometrika, 1950, 37, 149–157, ISSN 0006-3444
* C.F.Gauss, Theoria combinationis observationum erroribus minimis obnoxiae, 1821, Werke, 4. Gottinge

## Notes

1. ↑ Emannual C. Ifeacor, Barrie W. Jervis. Digital signal processing: a practical approach, second edition. Indianapolis: Pearson Education Limited, 2002, p. 718
2. ↑ Steven Van Vaerenbergh, Ignacio Santamaría, Miguel Lázaro-Gredilla "Estimation of the forgetting factor in kernel recursive least squares", 2012 IEEE International Workshop on Machine Learning for Signal Processing, 2012, accessed June 23, 2016.
3. ↑ Welch, Greg and Bishop, Gary "An Introduction to the Kalman Filter", Department of Computer Science, University of North Carolina at Chapel Hill, September 17, 1997, accessed July 19, 2011.
4. ↑ Diniz, Paulo S.R., "Adaptive Filtering: Algorithms and Practical Implementation", Springer Nature Switzerland AG 2020, Chapter 7: Adaptive Lattice-Based RLS Algorithms. https://doi.org/10.1007/978-3-030-29057-3_7
5. ↑ Albu, Kadlec, Softley, Matousek, Hermanek, Coleman, Fagan "Implementation of (Normalised) RLS Lattice on Virtex" , Digital Signal Processing, 2001, accessed December 24, 2011.

0.00

(0 votes)

Original source: https://en.wikipedia.org/wiki/Recursive least squares filter. Read more

|

Retrieved from "https://handwiki.org/wiki/index.php?title=Recursive_least_squares_filter&oldid=4227763"