Matrix-Riccati-Gleichung

Als Matrix-Riccati-Gleichungen oder algebraische Riccati-Gleichungen wird ein Typ von nichtlinearen Gleichungen für Matrizen bezeichnet, die sich bei Dimension 1 auf eine algebraische, quadratische Gleichung zurückführen lassen. Daher kommt auch die Bezeichnung des Problems in Anlehnung an die entsprechende Riccati-Differentialgleichung. Bei allgemeinen Dimensionen m , n ∈ N {\displaystyle m,n\in \mathbb {N} } ist in allgemeiner Form der Matrix-Riccati-Gleichung eine Matrix X ∈ R m × n {\displaystyle X\in \mathbb {R} ^{m\times n}} gesucht, welche die Gleichung

X B X + X A − D X − C = 0 ∈ R m × n {\displaystyle XBX+XA-DX-C=0\in \mathbb {R} ^{m\times n}}

erfüllt. Die anderen, vorgegebenen Matrizen haben die dazu passenden Dimensionen C , B T ∈ R m × n {\displaystyle C,B^{T}\in \mathbb {R} ^{m\times n}} , A ∈ R n × n {\displaystyle A\in \mathbb {R} ^{n\times n}} , D ∈ R m × m {\displaystyle D\in \mathbb {R} ^{m\times m}} . Ein Spezialfall dieser Gleichung ist X 2 = C {\displaystyle X^{2}=C} , welche als Lösungen die Quadratwurzel einer Matrix X = C 1 / 2 {\displaystyle X=C^{1/2}} hat, wenn solche existieren.

Außer bei der Quadratwurzel treten Matrix-Riccati-Gleichungen bei den folgenden Problemen auf.

Soll die ( m + n ) × ( m + 1 ) {\displaystyle (m+n)\times (m+1)} -Blockmatrix

M := ( A B C D )  mit  ( I n 0 X I m ) {\displaystyle M:={\begin{pmatrix}A&B\\C&D\end{pmatrix}}{\mbox{ mit }}{\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}}

auf obere Block-Dreieckform transformiert werden, bekommt man

( I n 0 − X I m ) ( A B C D ) ( I n 0 X I m ) = ( A + B X B 0 D − X B ) =: M ^ , {\displaystyle {\begin{pmatrix}I_{n}&0\\-X&I_{m}\end{pmatrix}}{\begin{pmatrix}A&B\\C&D\end{pmatrix}}{\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}={\begin{pmatrix}A+BX&B\\0&D-XB\end{pmatrix}}=:{\hat {M}},}

wenn X {\displaystyle X} Lösung der obigen Riccati-Gleichung ist, dann verschwindet der linke untere Block C + D X − X A − X B X = 0 {\displaystyle C+DX-XA-XBX=0} in der transformierten Matrix. Bei den beiden Einheitsmatrizen ist die Dimension als Index vermerkt, I k ∈ R k × k {\displaystyle I_{k}\in \mathbb {R} ^{k\times k}} . Die Multiplikation der 3 Matrizen stellt tatsächlich eine Ähnlichkeitstransformation dar, da der linke und der rechte Faktor zueinander invers sind. Daher ergeben sich die Eigenwerte der Gesamtmatrix M {\displaystyle M} aus der Vereinigung der Eigenwerte der beiden Hauptdiagonalblöcke A + B X {\displaystyle A+BX} und D − X B {\displaystyle D-XB} vom M ^ {\displaystyle {\hat {M}}} . Darüber hinaus bilden die ersten n {\displaystyle n} Spalten ( I n X ) {\displaystyle {\begin{pmatrix}I_{n}\\X\end{pmatrix}}} der Transformationsmatrix eine Basis für den zu A + B X {\displaystyle A+BX} gehörigen invarianten Unterraum (Summe von Eigenräumen) von M {\displaystyle M} , aus dem sich bei Bedarf die Eigenvektoren bestimmen lassen. Es gilt also

( A B C D ) ( I n X ) = ( I n X ) ( A + B X ) . {\displaystyle {\begin{pmatrix}A&B\\C&D\end{pmatrix}}{\begin{pmatrix}I_{n}\\X\end{pmatrix}}={\begin{pmatrix}I_{n}\\X\end{pmatrix}}(A+BX).}

Anwendung findet diese Eigenschaft z. B. bei der Nachbesserung von Eigenvektor-Basen: wenn M {\displaystyle M} durch Störungen aus einer Block-Dreieckmatrix hervorging, ist C {\displaystyle C} klein und unter geeigneten Voraussetzungen auch X {\displaystyle X} . Dann kann die Block-Dreieckform in der angegebenen Weise wiederhergestellt werden ([Stewart]).

Bei einem linearen System von Differentialgleichungen y ′ ( t ) = A y ( t ) + S u ( t ) {\displaystyle y'(t)=Ay(t)+Su(t)} für einen Zustand y ( t ) ∈ R n {\displaystyle y(t)\in \mathbb {R} ^{n}} mit konstanten Koeffizienten A ∈ R n × n {\displaystyle A\in \mathbb {R} ^{n\times n}} , S ∈ R n × p {\displaystyle S\in \mathbb {R} ^{n\times p}} soll diejenige optimale Steuerung u ( t ) ∈ R p {\displaystyle u(t)\in \mathbb {R} ^{p}} bestimmt werden, welche bei unendlichem Zeithorizont das Funktional

∫ 0 ∞ ( y ( t ) T Q y ( t ) + u ( t ) T R u ( t ) ) d t {\displaystyle \int _{0}^{\infty }{\big (}y(t)^{T}Qy(t)+u(t)^{T}Ru(t){\big )}dt}

minimiert. Darin ist R ∈ R p × p {\displaystyle R\in \mathbb {R} ^{p\times p}} symmetrisch und positiv definit, Q ∈ R n × n {\displaystyle Q\in \mathbb {R} ^{n\times n}} symmetrisch und positiv semi-definit. Verwendet man eine Steuerung durch Rückkopplung u ( t ) = − K y ( t ) {\displaystyle u(t)=-Ky(t)} , ist das Optimum bei unendlichem Zeithorizont gegeben durch u ( t ) = − R − 1 S T X y ( t ) {\displaystyle u(t)=-R^{-1}S^{T}Xy(t)} , wobei X = X T {\displaystyle X=X^{T}} die (maximale) symmetrische Lösung der Riccati-Gleichung

Q + X A + A T X − X S R − 1 S T X = 0 {\displaystyle Q+XA+A^{T}X-XSR^{-1}S^{T}X=0}

ist, für welche die Matrix A − S K = A − S R − 1 S T X {\displaystyle A-SK=A-SR^{-1}S^{T}X} asymptotisch stabil ist mit allen Eigenwerten in der linken komplexen Halbebene. Für mehr Hintergrund wird auf den Artikel LQ-Regler verwiesen. Diese Gleichung ist also ein Spezialfall der Gleichung aus der Einleitung mit m = n {\displaystyle m=n} , C = − Q {\displaystyle C=-Q} , D = − A T {\displaystyle D=-A^{T}} , B = − S R − 1 S T = B T {\displaystyle B=-SR^{-1}S^{T}=B^{T}} . Die hierzu gehörige Blockmatrix

L = ( A − S R − 1 S T − Q − A T ) {\displaystyle L={\begin{pmatrix}A&-SR^{-1}S^{T}\\-Q&-A^{T}\end{pmatrix}}}

ist eine hamiltonsche Matrix, da B {\displaystyle B} und C {\displaystyle C} hier symmetrisch sind. Bei dieser Matrix L {\displaystyle L} tritt mit jedem Eigenwert λ {\displaystyle \lambda } auch − λ {\displaystyle -\lambda } als Eigenwert auf.

Da die Matrix-Riccati-Gleichung eine algebraische Gleichung vom Grad 2 für die m ⋅ n {\displaystyle m\cdot n} Unbekannten in der Matrix X {\displaystyle X} ist, kann zur Lösung das Newton-Verfahren eingesetzt werden. Die Ableitung der Abbildung X ↦ X B X + X A − D X − C {\displaystyle X\mapsto XBX+XA-DX-C} an der Stelle X ∈ R m × n {\displaystyle X\in \mathbb {R} ^{m\times n}} ist die lineare Abbildung

H ↦ H ( A + B X ) + ( X B − D ) H  für  H ∈ R m × n . {\displaystyle H\mapsto H(A+BX)+(XB-D)H{\text{ für }}H\in \mathbb {R} ^{m\times n}.}

Mit einer aktuellen Näherung X k ∈ R m × n {\displaystyle X_{k}\in \mathbb {R} ^{m\times n}} bekommt man das Inkrement H k = X k + 1 − X k {\displaystyle H_{k}=X_{k+1}-X_{k}} zu einer verbesserten Näherung also aus dem linearen Gleichungssystem

H k ( A + B X k ) + ( X k B − D ) H k = C + D X k − X k A − X k B X k , k ≥ 0 , {\displaystyle H_{k}(A+BX_{k})+(X_{k}B-D)H_{k}=C+DX_{k}-X_{k}A-X_{k}BX_{k},\quad k\geq 0,}

wo auf der rechten Seite, wie gewohnt, das negative Residuum der Riccati-Gleichung steht. Das Ganze stellt eine Sylvester-Gleichung dar, im zugehörigen Artikel werden numerische Methoden zu ihrer Auflösung behandelt. Diese lineare Gleichung ist eindeutig lösbar, wenn die beiden Matrizen A + B X k {\displaystyle A+BX_{k}} und D − X k B {\displaystyle D-X_{k}B} keine gemeinsamen Eigenwerte besitzen, z. B. wenn die Realteile aller Eigenwerte von A + B X k {\displaystyle A+BX_{k}} oberhalb und die von D − X k B {\displaystyle D-X_{k}B} unterhalb eines geeigneten Wertes (etwa null) liegen.

Involutorische Matrizen V ∈ R N × N {\displaystyle V\in \mathbb {R} ^{N\times N}} sind Lösungen der einfachen Riccati-Gleichung V 2 = I {\displaystyle V^{2}=I} . Auch die Newton-Iteration für diese spezielle Gleichung ist einfach,

V k + 1 := 1 2 ( V k + V k − 1 ) ,   k = 0 , 1 , … , {\displaystyle V_{k+1}:={\frac {1}{2}}(V_{k}+V_{k}^{-1}),\ k=0,1,\ldots ,}

und man kann zeigen, dass diese Signum-Iteration immer und quadratisch konvergiert, sofern die Startmatrix V 0 := M ∈ R N × N {\displaystyle V_{0}:=M\in \mathbb {R} ^{N\times N}} keine rein imaginären Eigenwerte (einschließlich null) besitzt. Alle Matrizen V k , k ≥ 0 , {\displaystyle V_{k},\,k\geq 0,} kommutieren miteinander und besitzen daher die gleiche Jordan-Basis, und dies gilt auch für die Grenzwert-Matrix S ( M ) := lim k → ∞ V k {\displaystyle S(M):=\lim _{k\to \infty }V_{k}} . Die zugehörigen Eigenwerte der V k {\displaystyle V_{k}} konvergieren gegen 1 {\displaystyle 1} bzw. − 1 {\displaystyle -1} , wenn der Realteil im Eigenwert von V 0 = M {\displaystyle V_{0}=M} positiv bzw. negativ war. Daher besitzt S ( M ) {\displaystyle S(M)} nur die beiden Eigenwerte ± 1 {\displaystyle \pm 1} und wird als Signum-Funktion von M {\displaystyle M} bezeichnet, S ( M ) {\displaystyle S(M)} ist also eine Involution mit S 2 = I {\displaystyle S^{2}=I} . Da die Eigenwerte von S ( M ) {\displaystyle S(M)} bekannt sind, bekommt man Basen für die invarianten Unterräume zu + 1 {\displaystyle +1} bzw. − 1 {\displaystyle -1} , indem man Basen für die Kerne von S ( M ) − I {\displaystyle S(M)-I} bzw. S ( M ) + I {\displaystyle S(M)+I} bestimmt, etwa mit der QR-Zerlegung. Diese sind dann auch Basen für die invarianten Unterräume der Ausgangsmatrix M {\displaystyle M} zu den Eigenwerten mit positivem bzw. negativem Realteil.

Diesen Hintergrund kann man mit N = m + n {\displaystyle N=m+n} zur Lösung der ursprünglichen Riccati-Gleichung verwenden, wenn aufgrund der Struktur von M {\displaystyle M} bzw. H {\displaystyle H} die Zahl der Eigenwerte mit positivem und negativem Realteil klar ist. Das gilt für die Quadratwurzel und das Steuerungsproblem.

Für die Quadratwurzel verschwinden die Hauptdiagonalblöcke von M {\displaystyle M} und auch bei den V k {\displaystyle V_{k}} ist das so, sei also

M = ( 0 I C 0 ) = V 0 , V k = ( 0 R k P k 0 ) , k ≥ 0 , {\displaystyle M={\begin{pmatrix}0&I\\C&0\end{pmatrix}}=V_{0},\quad V_{k}={\begin{pmatrix}0&R_{k}\\P_{k}&0\end{pmatrix}},\,k\geq 0,}

mit N = 2 n {\displaystyle N=2n} . Die Iteration für die V k {\displaystyle V_{k}} lautet dann für die Einzelblöcke

P k + 1 = 1 2 ( P k + R k − 1 ) , R k + 1 := 1 2 ( R k + P k − 1 ) ,   k = 0 , 1 , … . {\displaystyle P_{k+1}={\frac {1}{2}}(P_{k}+R_{k}^{-1}),\quad R_{k+1}:={\frac {1}{2}}(R_{k}+P_{k}^{-1}),\ k=0,1,\ldots .}

Falls C {\displaystyle C} keine reellen und nicht-positiven Eigenwerte besitzt, konvergiert die Iteration gegen die eindeutige Wurzel lim k → ∞ P k = C 1 / 2 {\displaystyle \lim _{k\to \infty }P_{k}=C^{1/2}} , deren Eigenwert-Realteile positiv sind.

Bei der allgemeinen Gleichung ist die Signum-Funktion mit N = m + n {\displaystyle N=m+n} einsetzbar, wenn die Riccati-Gleichung eine Lösung X {\displaystyle X} besitzt, für die A + B X {\displaystyle A+BX} und X B − D {\displaystyle XB-D} asymptotisch stabil sind, also beide nur Eigenwerte mit negativem Realteil besitzen. Unter dieser Voraussetzung ist S ( A + B X ) = − I n {\displaystyle S(A+BX)=-I_{n}} und S ( D − X B ) = + I m {\displaystyle S(D-XB)=+I_{m}} und für die Blockmatrix M {\displaystyle M} folgt, dass

S ( M ) ( I n 0 X I m ) = ( I n 0 X I m ) ( − I n G 0 I m ) {\displaystyle S(M){\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}={\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}{\begin{pmatrix}-I_{n}&G\\0&I_{m}\end{pmatrix}}}

ist mit einer geeigneten Matrix G {\displaystyle G} . Die ersten n {\displaystyle n} Spalten dieser Gleichung zeigen mit

( S ( M ) + I N ) ( I n X ) = 0 , {\displaystyle {\big (}S(M)+I_{N}{\big )}{\begin{pmatrix}I_{n}\\X\end{pmatrix}}=0,}

dass die Matrix ( I n X ) {\displaystyle {\begin{pmatrix}I_{n}\\X\end{pmatrix}}} eine spezielle Basis des Kerns von S ( M ) + I N {\displaystyle S(M)+I_{N}} ist. Zur Lösung der Riccati-Gleichung sind also mit der Startmatrix V 0 := M {\displaystyle V_{0}:=M} , bzw. V 0 := L {\displaystyle V_{0}:=L} bei der optimalen Steuerung, die Matrizen V k {\displaystyle V_{k}} und ihr Grenzwert S ( M ) {\displaystyle S(M)} zu berechnen. Danach bekommt man X {\displaystyle X} bei Aufteilung von S ( M ) + I N {\displaystyle S(M)+I_{N}} in Blöcke aus dem folgenden Gleichungssystem

( G 12 G 22 ) X = − ( G 11 G 21 )  mit  S ( M ) + I N = ( G 11 G 12 G 21 G 22 ) . {\displaystyle {\begin{pmatrix}G_{12}\\G_{22}\end{pmatrix}}X=-{\begin{pmatrix}G_{11}\\G_{21}\end{pmatrix}}{\text{ mit }}S(M)+I_{N}={\begin{pmatrix}G_{11}&G_{12}\\G_{21}&G_{22}\end{pmatrix}}.}

Hier sind G 11 ∈ R n × n {\displaystyle G_{11}\in \mathbb {R} ^{n\times n}} , G 12 , G 21 T ∈ R n × m {\displaystyle G_{12},G_{21}^{T}\in \mathbb {R} ^{n\times m}} , G 22 ∈ R m × m {\displaystyle G_{22}\in \mathbb {R} ^{m\times m}} .

Bei der Anwendung zur optimalen Steuerung sei mit n = 2 {\displaystyle n=2} und p = 1 {\displaystyle p=1} ,

A = ( − 2 3 − 2 − 1 − 8 3 ) ,   S = ( 1 1 2 ) ,   Q = ( 1 5 3 5 3 20 3 ) , {\displaystyle A={\begin{pmatrix}-{\frac {2}{3}}&-2\\-1&-{\frac {8}{3}}\end{pmatrix}},\ S={\begin{pmatrix}1\\{\frac {1}{2}}\end{pmatrix}},\ Q={\begin{pmatrix}1&{\frac {5}{3}}\\{\frac {5}{3}}&{\frac {20}{3}}\end{pmatrix}},}

sowie R = 1 ∈ R 1 × 1 {\displaystyle R=1\in \mathbb {R} ^{1\times 1}} . Von den Eigenwerten − 5 3 ± 3 {\displaystyle -{\frac {5}{3}}\pm {\sqrt {3}}} der Matrix A {\displaystyle A} ist einer positiv, das ungeregelte System mit u ≡ 0 {\displaystyle u\equiv 0} ist also instabil. Als Blockmatrix M {\displaystyle M} tritt hier die speziellere Form

L = ( − 2 3 − 2 − 1 − 1 2 − 1 − 8 3 − 1 2 − 1 4 − 1 − 5 3 2 3 1 − 5 3 − 20 3 2 8 3 ) {\displaystyle L={\begin{pmatrix}-{\frac {2}{3}}&-2&-1&-{\frac {1}{2}}\\-1&-{\frac {8}{3}}&-{\frac {1}{2}}&-{\frac {1}{4}}\\-1&-{\frac {5}{3}}&{\frac {2}{3}}&1\\-{\frac {5}{3}}&-{\frac {20}{3}}&2&{\frac {8}{3}}\end{pmatrix}}}

auf, sie besitzt die 4 Eigenwerte ± 13 6 ± 1 2 13 {\displaystyle \pm {\frac {13}{6}}\pm {\frac {1}{2}}{\sqrt {13}}} , von denen, wie erwähnt, tatsächlich 2 positiv und 2 negativ sind. In diesem Beispiel lässt sich die Signum-Funktion von L {\displaystyle L} noch über deren Jordan-Normalform berechnen, das Ergebnis ist

S ( L ) = ( 25 338 − 135 169 − 114 169 21 338 − 75 338 − 115 169 21 338 − 87 676 − 789 676 163 338 − 25 338 75 338 163 338 − 433 169 135 169 115 169 ) . {\displaystyle S(L)={\begin{pmatrix}{\frac {25}{338}}&-{\frac {135}{169}}&-{\frac {114}{169}}&{\frac {21}{338}}\\-{\frac {75}{338}}&-{\frac {115}{169}}&{\frac {21}{338}}&-{\frac {87}{676}}\\-{\frac {789}{676}}&{\frac {163}{338}}&-{\frac {25}{338}}&{\frac {75}{338}}\\{\frac {163}{338}}&-{\frac {433}{169}}&{\frac {135}{169}}&{\frac {115}{169}}\end{pmatrix}}.}

Tatsächlich kann man direkt verifizieren, dass S ( L ) {\displaystyle S(L)} involutorisch ist, S ( L ) 2 = I {\displaystyle S(L)^{2}=I} , und mit L {\displaystyle L} kommutiert, L S ( L ) − S ( L ) L = 0 {\displaystyle L\,S(L)-S(L)\,L=0} . Eine Basismatrix Y {\displaystyle Y} des Kerns von S ( L ) + I 4 {\displaystyle S(L)+I_{4}} , also mit ( S ( L ) + I 4 ) Y = 0 {\displaystyle (S(L)+I_{4})Y=0} ist gegeben durch

Y = ( 1 0 3 2 1 0 − 1 2 2 ) = ( 1 0 0 1 3 2 − 1 − 1 2 ) ( 1 0 3 2 1 ) = ( I 2 X ) ( 1 0 3 2 1 ) . {\displaystyle Y={\begin{pmatrix}1&0\\{\frac {3}{2}}&1\\0&-1\\2&2\end{pmatrix}}=\left({\begin{array}{cc}1&0\\0&1\\\hline {\frac {3}{2}}&-1\\-1&2\end{array}}\right){\begin{pmatrix}1&0\\{\frac {3}{2}}&1\end{pmatrix}}={\begin{pmatrix}I_{2}\\X\end{pmatrix}}{\begin{pmatrix}1&0\\{\frac {3}{2}}&1\end{pmatrix}}.}

Durch spaltenweise Elimination in den ersten beiden Zeilen von Y {\displaystyle Y} wurde dort eine Einheitsmatrix erzeugt und man kann daher im unteren Block die Lösung X {\displaystyle X} der Riccati-Gleichung ablesen mit

X = ( 3 2 − 1 − 1 2 ) ,  und es gilt  A + B X = A − S R − 1 S T X = ( − 5 3 − 2 − 3 2 − 8 3 ) . {\displaystyle X={\begin{pmatrix}{\frac {3}{2}}&-1\\-1&2\end{pmatrix}},{\mbox{ und es gilt }}A+BX=A-SR^{-1}S^{T}X={\begin{pmatrix}-{\frac {5}{3}}&-2\\-{\frac {3}{2}}&-{\frac {8}{3}}\end{pmatrix}}.}

Die gesteuerte Systemmatrix A + B X {\displaystyle A+BX} hat jetzt also 2 negative Eigenwerte und das System ist daher asymptotisch stabil.

Die Berechnung der Jordan-Normalform umgeht man mit der in Abschnitt 2.2 beschriebenen Signum-Iteration. Die Konvergenz V k → S ( L ) , k ≥ 0 , {\displaystyle V_{k}\to S(L),\,k\geq 0,} ist quadratisch, man kann dies direkt an den Eigenwerten der Matrizen V k {\displaystyle V_{k}} ablesen. Diese lauten:

k = Eigenwerte  V k 0 ± 0.363891029 ± 3.969442305 1 ± 1.555983235 ± 2.110683431 2 ± 1.099331841 ± 1.292231811 3 ± 1.004487641 ± 1.033043387 4 ± 1.000010024 ± 1.000528470 5 ± 1.000000000 ± 1.000000140 6 ± 1.000000000 ± 1.000000000 {\displaystyle {\begin{array}{c|cc}k=&{\text{Eigenwerte }}V_{k}\\\hline 0&\pm 0.363891029&\pm 3.969442305\\1&\pm 1.555983235&\pm 2.110683431\\2&\pm 1.099331841&\pm 1.292231811\\3&\pm 1.004487641&\pm 1.033043387\\4&\pm 1.000010024&\pm 1.000528470\\5&\pm 1.000000000&\pm 1.000000140\\6&\pm 1.000000000&\pm 1.000000000\end{array}}}

Tatsächlich ist ‖ V 6 2 − I 4 ‖ ≈ 10 − 14 {\displaystyle \|V_{6}^{2}-I_{4}\|\approx 10^{-14}} . Setzt man zur Berechnung der Lösung X {\displaystyle X} an Stelle von S ( L ) {\displaystyle S(L)} die Näherung V 6 {\displaystyle V_{6}} ein und teilt V 6 + I 4 {\displaystyle V_{6}+I_{4}} auf wie beschrieben,

V 6 + I N = ( G 11 G 12 G 21 G 22 ) {\displaystyle V_{6}+I_{N}={\begin{pmatrix}G_{11}&G_{12}\\G_{21}&G_{22}\end{pmatrix}}}

bekommt man mit Hilfe der reduzierten QR-Zerlegung

( G 12 G 22 ) = Q ^ ⋅ R ^ ≈ ( − 0.48245 0.43832 0.04444 − 0.13345 0.66238 − 0.36922 0.57138 0.80854 ) ( 1.39805 1.07147 0 1.32121 ) {\displaystyle {\begin{pmatrix}G_{12}\\G_{22}\end{pmatrix}}={\hat {Q}}\cdot {\hat {R}}\approx {\begin{pmatrix}-0.48245&0.43832\\0.04444&-0.13345\\0.66238&-0.36922\\0.57138&0.80854\end{pmatrix}}{\begin{pmatrix}1.39805&1.07147\\0&1.32121\end{pmatrix}}}

(Angabe aus Platzgründen mit geringer Genauigkeit) die Näherungslösung

X ~ = − R ^ − 1 Q ^ T ( G 11 G 21 ) = ( 1.499999999 − 0.9999999997 − 1.000000000 2.000000001 ) . {\displaystyle {\tilde {X}}=-{\hat {R}}^{-1}{\hat {Q}}^{T}{\begin{pmatrix}G_{11}\\G_{21}\end{pmatrix}}={\begin{pmatrix}1.499999999&-0.9999999997\\-1.000000000&2.000000001\end{pmatrix}}.}

Diese Näherung ist offensichtlich auf ca. 9 Stellen genau.

• G.W. Stewart, Error and perturbation bounds for subspaces associated with certain eigenvalue problems, SIAM Review 15, 727–764
• N.J. Higham, Functions of matrices: Theory and computation, SIAM, Philadelphia, 2008.
• J.D. Roberts, Linear model reduction and solution of the algebraic Riccati equation by use of the sign function, Intern. J. Control 32, 677–687