RLS algorithm

The RLS algorithm ( R ecursive- L EAST S quares- algorithm ) based on the method of least squares . It is used to solve overdetermined linear systems of equations and especially to estimate model parameters when identifying linear systems or in neuro-informatics . The recursivity allows online use with the latest data while maintaining the same complexity in each recursion step.

Derivation

The basis for the recursive algorithm is first of all the formal quadratic minimization problem. In the example of the system identification there are pairs of input and output data in the matrix , the relationship of which is to be mapped using a linear model with the parameters to be estimated . With the additional measured output values in the vector , the system can be described by: ${\displaystyle \mathbf {X} }$ ${\displaystyle {\underline {\hat {\Theta }}}}$ ${\displaystyle {\underline {y}}}$

${\displaystyle {\underline {y}}=\mathbf {X} \cdot {\underline {\hat {\Theta }}}}$

The corresponding optimization problem is formulated as follows:

${\displaystyle {\underset {\underline {\hat {\Theta }}}{\min }}f({\underline {\hat {\Theta }}})={\underset {\underline {\hat {\Theta }}}{\min }}{\frac {1}{2}}||{\underline {y}}-\mathbf {X} \cdot {\underline {\hat {\Theta }}}||^{2}={\underset {\underline {\hat {\Theta }}}{\min }}\left[{\frac {1}{2}}\left({\underline {y}}-\mathbf {X} \cdot {\underline {\hat {\Theta }}}\right)^{T}\left({\underline {y}}-\mathbf {X} \cdot {\underline {\hat {\Theta }}}\right)\right]}$

The square of the difference between the measurement and model data should be minimized. This procedure has the advantage that a quadratic function has exactly one global extreme at its vertex. At this point the first derivative becomes zero, which leads to the solution of the minimization problem:

${\displaystyle {\frac {\partial f({\underline {\hat {\Theta }}})}{\partial {\underline {\hat {\Theta }}}}}=-\mathbf {X} ^{T}\cdot \left({\underline {y}}-\mathbf {X} \cdot {\underline {\hat {\Theta }}}\right){\overset {!}{=}}0}$

Changing provides the parameter vector you are looking for:

${\displaystyle {\underline {\hat {\Theta }}}=\left(\mathbf {X} ^{T}\cdot \mathbf {X} \right)^{-1}\cdot \mathbf {X} ^{T}\cdot {\underline {y}}}$

To solve this, the number of data pairs must at least correspond to the number of parameters sought. The more data pairs there are, the more entries the matrix has and the more difficult it is to calculate. ${\displaystyle \mathbf {X} }$

Recursion

When switching to recursive calculation, the calculation effort remains the same even if new data is added in every step, since the previous result is used as the starting point and thus only one new data pair is involved in the calculation.

The system of equations to be solved is in the recursion step , where at least corresponds to the number of parameters ${\displaystyle k}$ ${\displaystyle k}$

${\displaystyle {\underline {y}}(k)=\mathbf {X} (k)\cdot {\underline {\hat {\Theta }}}(k)}$

with the corresponding solution for the parameter vector : ${\displaystyle {\underline {\hat {\Theta }}}(k)}$

${\displaystyle {\underline {\hat {\Theta }}}(k)=\left(\mathbf {X} ^{T}(k)\cdot \mathbf {X} (k)\right)^{-1}\cdot \mathbf {X} ^{T}(k)\cdot {\underline {y}}(k)}$

For the system of equations expands to: ${\displaystyle k+1}$

${\displaystyle {\underline {y}}(k+1)=\mathbf {X} (k+1)\cdot {\underline {\hat {\Theta }}}(k+1)}$

and for the solution applies accordingly:

${\displaystyle {\underline {\hat {\Theta }}}(k+1)=\left(\mathbf {X} ^{T}(k+1)\cdot \mathbf {X} (k+1)\right)^{-1}\cdot \mathbf {X} ^{T}(k+1)\cdot y(k+1)}$

The data in the step can be divided into existing and newly added data as follows: ${\displaystyle k+1}$

${\displaystyle {\underline {y}}(k+1)={\begin{bmatrix}{\underline {y}}(k)\\y(k+1)\end{bmatrix}},\qquad \mathbf {X} (k+1)={\begin{bmatrix}\mathbf {X} (k)\\{\underline {x}}^{T}(k+1)\end{bmatrix}}}$

Insertion into the equation for the parameter vector yields:

${\displaystyle {\underline {\hat {\Theta }}}(k+1)=\left({\begin{bmatrix}\mathbf {X} ^{T}(k)&{\underline {x}}(k+1)\end{bmatrix}}\cdot {\begin{bmatrix}\mathbf {X} (k)\\{\underline {x}}^{T}(k+1)\end{bmatrix}}\right)^{-1}\cdot {\begin{bmatrix}\mathbf {X} ^{T}(k)&{\underline {x}}(k+1)\end{bmatrix}}\cdot {\begin{bmatrix}{\underline {y}}(k)\\y(k+1)\end{bmatrix}}}$

The expression can be further classified into:

${\displaystyle {\underline {\hat {\Theta }}}(k+1)=\left(\mathbf {X} ^{T}(k)\cdot \mathbf {X} (k)+{\underline {x}}(k+1)\cdot {\underline {x}}^{T}(k+1)\right)^{-1}\cdot \left(\mathbf {X} ^{T}(k)\cdot \mathbf {X} (k)\cdot {\underline {\hat {\Theta }}}(k)+{\underline {x}}(k+1)\cdot y(k+1)\right)}$

With the abbreviation and the use of the lemma for matrix inversion , the recursion equation results in: ${\displaystyle \mathbf {P} ^{-1}=\mathbf {X} ^{T}(k)\cdot \mathbf {X} (k)}$

${\displaystyle {\underline {\hat {\Theta }}}(k+1)=\left(\mathbf {P} (k)-{\frac {\mathbf {P} (k)\cdot {\underline {x}}(k+1)\cdot {\underline {x}}^{T}(k+1)\cdot \mathbf {P} (k)}{1+{\underline {x}}^{T}(k+1)\cdot \mathbf {P} (k)\cdot {\underline {x}}(k+1)}}\right)\cdot \left(\mathbf {P} ^{-1}(k)\cdot {\underline {\hat {\Theta }}}(k)+{\underline {x}}(k+1)\cdot y(k+1)\right)}$

Simplified:

${\displaystyle {\underline {\hat {\Theta }}}(k+1)={\underline {\hat {\Theta }}}(k)+\underbrace {\frac {\mathbf {P} (k)\cdot {\underline {x}}(k+1)}{1+{\underline {x}}^{T}(k+1)\cdot \mathbf {P} (k)\cdot {\underline {x}}(k+1)}} _{{\text{Verstärkungsterm }}{\underline {\gamma }}(k)}\cdot \underbrace {\left(y(k+1)-{\underline {x}}^{T}(k+1)\cdot {\underline {\hat {\Theta }}}(k)\right)} _{\text{Korrekturterm}}}$

The matrix can also be updated recursively: ${\displaystyle \mathbf {P} }$

${\displaystyle \mathbf {P} (k+1)=\mathbf {P} (k)-{\underline {\gamma }}(k)\cdot {\underline {x}}^{T}(k+1)\cdot \mathbf {P} (k)}$

In comparison to the non-recursive algorithm, it can be seen that an inversion of the matrix is no longer necessary, only a division by a scalar. ${\displaystyle \left(\mathbf {X} ^{T}\cdot \mathbf {X} \right)}$

Forgetting factor

By introducing a forgetting factor, the historical measurement data lose their importance for the optimization and the current ones are weighted more heavily. This is useful to compensate for changes in the operating point, disturbances or invariance of the system to be modeled. Usually lambda is chosen in the range . ${\displaystyle \lambda <1}$ ${\displaystyle 0{,}95<\lambda <1}$

${\displaystyle {\underline {\hat {\Theta }}}(k+1)={\underline {\hat {\Theta }}}(k)+\underbrace {\frac {\mathbf {P} (k)\cdot {\underline {x}}(k+1)}{\lambda +{\underline {x}}^{T}(k+1)\cdot \mathbf {P} (k)\cdot {\underline {x}}(k+1)}} _{{\text{Verstärkungsterm }}{\underline {\gamma }}(k)}\cdot \underbrace {\left(y(k+1)-{\underline {x}}^{T}(k+1)\cdot {\underline {\hat {\Theta }}}(k)\right)} _{\text{Korrekturterm}}}$

${\displaystyle \mathbf {P} (k+1)={\frac {1}{\lambda }}\cdot \left(\mathbf {P} (k)-{\underline {\gamma }}(k)\cdot {\underline {x}}^{T}(k+1)\cdot \mathbf {P} (k)\right)}$

application

For a good result, the RLS algorithm requires a maximum of as many recursion steps as parameters are to be identified. This speed and its simple calculability enable a wide range of online application options for system identification or as an adaptive filter on less potent microprocessor systems . Both and must be initialized at the beginning of the recursion . If information about the parameters is available a priori, it can be used here, otherwise the parameters are set to zero. As a rule, an initialization as a diagonal matrix with values ​​greater than 100 is suitable for the matrix . High values ​​mean little confidence in the parameters, which leads to greater parameter corrections, which are initially quite desirable. ${\displaystyle {\underline {\hat {\Theta }}}}$ ${\displaystyle \mathbf {P} }$ ${\displaystyle \mathbf {P} }$

In order for the process to remain stable, the matrix must remain positive and definite . However, numerical inaccuracies during the calculation can result in negative eigenvalues. Therefore, implementations of the RLS algorithm have been developed that prove to be numerically more stable. B. can be achieved by integrating a QR decomposition . However, this increases the computational effort. ${\displaystyle \mathbf {P} }$

example

For a system to be identified, a PT2 system was chosen as a model function, whose time-discrete implementation is available in the form of the following recursion equation:

${\displaystyle {\hat {y}}(k)=-a_{1}\cdot {\hat {y}}(k-1)-a_{2}\cdot {\hat {y}}(k-2)+a_{3}\cdot u(k-1)+a_{4}\cdot u(k-2)}$

With the combination of model input and output to and the parameter vector, the matrix notation follows: ${\displaystyle {\underline {h}}^{T}={\begin{bmatrix}-{\hat {y}}(k-1)&-{\hat {y}}(k-2)&u(k-1)&u(k-2)\end{bmatrix}}}$ ${\displaystyle {\underline {\hat {\Theta }}}^{T}={\begin{bmatrix}a_{1}&a_{2}&a_{3}&a_{4}\end{bmatrix}}}$

${\displaystyle {\hat {y}}(k)={\underline {h}}^{T}\cdot {\underline {\hat {\Theta }}}}$

The following implementation of the RLS algorithm now results for the selected model approach:

${\displaystyle {\underline {\hat {\Theta }}}(k+1)={\underline {\hat {\Theta }}}(k)+{\frac {\mathbf {P} (k)\cdot {\underline {h}}(k+1)}{\lambda +{\underline {h}}^{T}(k+1)\cdot \mathbf {P} (k)\cdot {\underline {h}}(k+1)}}\cdot \left(y(k+1)-{\underline {h}}^{T}(k+1)\cdot {\underline {\hat {\Theta }}}(k)\right)}$

${\displaystyle \mathbf {P} (k+1)={\frac {1}{\lambda }}\cdot \left(\mathbf {P} (k)-{\underline {\gamma }}(k)\cdot {\underline {h}}^{T}(k+1)\cdot \mathbf {P} (k)\right)}$

literature

• Dierk Schröder: Intelligent processes . Springer-Verlag, Berlin, Heidelberg 2010, ISBN 978-3-642-11397-0 .
• Martin Werner: Digital signal processing with MATLAB® internship . Vieweg & Sohn Verlag, Wiesbaden 2008, ISBN 978-3-8348-0393-1 .