RLS-Algorithmus

Der RLS-Algorithmus (Recursive-Least-Squares-Algorithmus) basiert auf der Methode der kleinsten Quadrate. Er wird zur Lösung überbestimmter linearer Gleichungssysteme und insbesondere zur Schätzung von Modellparametern bei der Identifikation linearer Systeme oder in der Neuroinformatik genutzt. Die Rekursivität erlaubt die Online-Nutzung mit aktuell anfallenden Daten bei gleichbleibender Komplexität in jedem Rekursionsschritt.

Herleitung

Die Basis für den rekursiven Algorithmus ist zunächst das formale quadratische Minimierungsproblem. Im Beispiel der Systemidentifikation liegen Paare aus Ein- und Ausgangsdaten in der Matrix $\mathbf {X}$ vor, deren Zusammenhang durch ein lineares Modell mit den zu schätzenden Parametern ${\underline {\hat {\Theta }}}$ abgebildet werden soll. Mit den zusätzlich im Vektor ${\underline {y}}$ vorliegenden gemessenen Ausgangswerten lässt sich das System beschreiben durch:

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

Das entsprechende Optimierungsproblem formuliert sich zu:

${\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]$

Es soll das Quadrat der Differenz aus Mess- und Modelldaten minimiert werden. Dieses Vorgehen hat den Vorteil, dass eine quadratische Funktion genau ein globales Extremum in ihrem Scheitelpunkt besitzt. An dieser Stelle wird die erste Ableitung zu Null, was zur Lösung des Minimierungsproblems führt:

${\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$

Umstellen liefert den gesuchten Parametervektor:

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

Zur Lösung muss die Anzahl der Datenpaare mindestens der Anzahl der gesuchten Parameter entsprechen. Je mehr Datenpaare vorliegen, desto mehr Einträge hat die Matrix $\mathbf {X}$ und desto schwieriger gestaltet sich die Berechnung.

Rekursion

Beim Übergang zur rekursiven Berechnung bleibt auch bei Hinzukommen neuer Daten in jedem Schritt der Rechenaufwand gleich, da das vorherige Ergebnis als Ausgangspunkt genutzt wird und so nur je ein neues Datenpaar in die Kalkulation involviert ist.

Im Rekursionsschritt $k$, wobei $k$ mindestens der Anzahl der Parameter entspricht, ist das zu lösende Gleichungssystem

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

mit der zugehörigen Lösung für den Parametervektor ${\underline {\hat {\Theta }}}(k)$:

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

Für $k+1$ erweitert sich das Gleichungssystem zu:

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

und für die Lösung gilt entsprechend:

${\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)$

Die Daten im Schritt $k+1$ lassen sich in bereits vorliegende und neu hinzugekommene Daten folgendermaßen aufteilen:

${\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}}$

Einsetzen in die Gleichung für den Parametervektor liefert:

${\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}}$

Der Ausdruck lässt sich noch weiter ordnen zu:

${\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)$

Mit der Abkürzung $\mathbf {P} ^{-1}=\mathbf {X} ^{T}(k)\cdot \mathbf {X} (k)$ und der Nutzung des Lemmas zur Matrixinversion ergibt sich die Rekursionsgleichung zu:

${\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)$

Vereinfacht:

${\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}}$

Die Aktualisierung der Matrix $\mathbf {P}$ kann ebenfalls rekursiv erfolgen:

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

Im Vergleich zum nichtrekursiven Algorithmus ist ersichtlich, dass keine Inversion der Matrix $\left(\mathbf {X} ^{T}\cdot \mathbf {X} \right)$ mehr notwendig ist, sondern lediglich eine Division durch ein Skalar.

Vergessensfaktor

Durch die Einführung eines Vergessensfaktors $\lambda <1$ verlieren die historischen Messdaten für die Optimierung an Bedeutung und die aktuellen werden stärker gewichtet. Dies ist sinnvoll, um Arbeitspunktwechsel, Störungen oder Invarianzen des zu modellierenden Systems auszugleichen. Üblicherweise wird Lambda im Bereich $0{,}95<\lambda <1$ gewählt.

${\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}}$

$\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)$

Anwendung

Der RLS-Algorithmus benötigt für ein gutes Ergebnis maximal so viele Rekursionsschritte wie Parameter identifiziert werden sollen. Diese Schnelligkeit und seine einfache Berechenbarkeit ermöglichen vielfältige Online-Anwendungsmöglichkeiten zur Systemidentifikation oder als adaptives Filter auch auf weniger potenten Mikroprozessorsystemen. Zu Beginn der Rekursion müssen sowohl ${\underline {\hat {\Theta }}}$ als auch $\mathbf {P}$ initialisiert werden. Liegen a priori Informationen zu den Parametern vor, so können diese hier verwendet werden, ansonsten werden die Parameter zu Null gesetzt. Für die Matrix $\mathbf {P}$ eignet sich in der Regel eine Initialisierung als Diagonalmatrix mit Werten größer 100. Hohe Werte bedeuten geringes Vertrauen in die Parameter, was zu stärkeren Parameterkorrekturen führt, die anfänglich durchaus erwünscht sind.

Damit das Verfahren stabil bleibt, muss die Matrix $\mathbf {P}$ positiv definit bleiben. Durch numerische Ungenauigkeiten während der Berechnung können jedoch negative Eigenwerte entstehen. Deshalb wurden Implementierungen des RLS-Algorithmus entwickelt, die sich als numerisch stabiler erweisen, was z. B. durch die Einbindung einer QR-Zerlegung erreicht werden kann. Allerdings steigt hierbei der Rechenaufwand.

Beispiel

Für ein zu identifizierendes System wurde ein PT2-System als Modellfunktion gewählt, dessen zeitdiskrete Realisierung in Form der folgenden Rekursionsgleichung vorliegt:

${\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)$

Mit der Zusammenfassung von Modellein- und -ausgang zu ${\underline {h}}^{T}={\begin{bmatrix}-{\hat {y}}(k-1)&-{\hat {y}}(k-2)&u(k-1)&u(k-2)\end{bmatrix}}$ und dem Parametervektor ${\underline {\hat {\Theta }}}^{T}={\begin{bmatrix}a_{1}&a_{2}&a_{3}&a_{4}\end{bmatrix}}$ folgt die matrizielle Schreibweise:

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

Für den gewählten Modellansatz ergibt sich nun folgende Realisierung des RLS-Algorithmus:

${\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)$

$\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)$

Literatur

• Dierk Schröder: Intelligente Verfahren. Springer-Verlag, Berlin, Heidelberg 2010, ISBN 978-3-642-11397-0. 
• Martin Werner: Digitale Signalverarbeitung mit MATLAB®-Praktikum. Vieweg & Sohn Verlag, Wiesbaden 2008, ISBN 978-3-8348-0393-1.