Lineare Ausgleichsrechnung, QR-Zerlegung

Gegeben: Datenpunkte (tj, yj), j = 1, . . . , M Ferner seien vorgegeben “Ansatz-funktionen”ϕ₁(t), . . . , ϕ_n(t), z.B. ϕ_i(t) = tⁱ⁻¹.

Es sei M ≥n. (In der Regel ist M n.) Gesucht sind Koeffizienten (x₁, . . . , x_n), so daß

( _M X

j=1

(y_j −(

n

X

i=1

x_i ·ϕ_i(t_j)))² )

= min

x1,...,xn

d.h. die Summe der Abweichungsquadrate der Ordinaten ist zu minimieren.

Mit den Setzungen r_j =y_j −

n

X

i=1

x_iϕ_i(t_j) also~r=





 r₁

... r_M





 Residuenvektor ist also folgendes Problem zu l¨osen

M

X

j=1

rj2

=k~rk²₂ = min^!

x1,...,xn

.

4.6.1 L¨osungsansatz mittels Differentialrechnung: Gauß-sche Normal-gleichungen

Man setze

s(x₁, . . . , x_n) :=

M

X

j=1

(y_j −(

n

X

i=1

x_i·ϕ_i(t_j)))² (=k~rk²₂)

(sist also die ”Fehlerquadratsumme”). Das notwendige Extremalkriterium lautet:

∂

∂x_ks(x₁, . . . , x_n) = 0 f¨ur k = 1, . . . , n

Mit Hilfe der Kettenregel ergibt sich die partielle Ableitung

∂

∂x_ks(x₁, . . . , x_n) =

M

X

j=1

2 (y_j −(

n

X

i=1

x_i·ϕ_i(t_j)))

| {z }

=rj

(−ϕ_k(t_j)) = 0 f¨ur k = 1, . . . , n

Setzt man

~ ϕk =







ϕ_k(t₁) ... ϕ_k(t_M)







dann ergibt sich diese partielle Ableitung aus

−2·PM

j=1r_j·ϕ_k(t_j) = 0 |: (−2)

⇔

M

X

j=1

rj ·ϕk(tj) = 0

⇔ ~r^T ·ϕ~_k = 0 f¨ur k = 1, . . . , n Wir f¨uhren die Matrix aus den Ansatzfunktionen, ausgewertet auf dem Gitter der t₁, . . . , t_M, ein:

Φ = (~ϕ1, . . . , ~ϕn)∈R^M×n . Dann liest sich die notwendige Extremalbedingung als

Φ^T~r= 0∈Rⁿ Mit

~y =





 y₁

... y_M





 ~x=





 x₁

... x_n





 und Φ^T~r = 0 erhalten wir

~r = ~y−Φ~x ⇒(~y−Φ~x)^TΦ =~0 |()^T

⇔Φ^T(~y−Φ~x) =~0 |+ Φ^TΦ~x

⇔ Φ^TΦ~x= Φ^T~y Gauß’sches

Normalgleichungssystem

Die Matrix Φ^TΦ ∈ R^n×n ist symmetrisch und ~x ∈ Rⁿ, Φ^T~y ∈ Rⁿ. Falls die Spaltenϕ~₁, . . . , ~ϕ_n von Φ linear unabh¨angig sind, dann gilt:

Φ~u=~0⇔~u=~0 .

Dann ist Φ^TΦ positiv definit, d.h. ~u^T(Φ^TΦ)~u >0 f¨ur ~u6=~0. Beweis:

~

u^T(Φ^TΦ)~u = (~u^TΦ^T)(Φ~u) = (Φ~u)^T(Φ~u) =kΦ~uk²₂>0

⇔ Φ~u6=~0∀~u6=~0 .

~xergibt~r als “optimales Residuum” mit

Φ^T~r =~0 d.h.~r⊥{~ϕ1, . . . , ~ϕn} .

D.h. dasoptimale Residuum ist orthogonal zum Bildraum von Φ, (der Menge aller Linearkombinationen der Spalten von Φ.). Wir betrachten die 2. Ableitung:

∂² hinreichende Optimalkriterium erf¨ullt.

Zur Erinnerung: Ists: Rⁿ→R,∇s= (_∂x^∂s

positiv definit, so istxeine strenge lokale Minimalstelle von s.

Beispiel:

Die gesuchte Funktion ist also 3.157t+ 0.65.

Dieser eigentlich elegante Zugang leidet unter einem Problem: Die Matrix Φ^TΦ ist oft sehr schlecht konditioniert. Die Rundungsfehler bei der Aufstellung der Normalgleichungen haben einen großen Einfluß auf die berechnete L¨osung.

Wir stellen uns folgende Frage:

Kann man den Cholesky-Faktor L^T der Zerlegung Φ^TΦ = LL^T (mit L = ...

... ) direkt berechnen, ohne Φ^TΦ zu bilden? Die Antwort lautet ”ja” und der n¨achste Abschnitt zeigt den Weg dazu. Man kann dann die schlechte Kondition des Nor-malgleichungssystems weitgehend vermeiden.

4.6.2 QR-Zerlegung

Annahme: Gegeben sei eine unit¨are M ×M Matrix Q, so daß Q·Φ = ^R₀

∈ R^M×n.

M M M

M n n

0 R

...

...

... · ...

...

... = ...

0 Dann folgt

Φ =Q^T R

0

, weil Q unit¨ar, gilt Q⁻¹ =Q^T und

Φ^TΦ = (R^T|0)QQ^T

| {z }

IM

R 0

=R^TR .

Bemerkung 4.11. Mit Φ^TΦ = LL^T (Cholesky) folgt hieraus nun R = DL^T, wobei D eine Diagonalmatrix mit Elementen ±1 (im Reellen) ist. 2

Definition 4.9. Sei Q∈ R^M×M unit¨ar und Φ∈R^M×n mit M ≥n. Eine Zerle-gung der Form

QΦ =



 R

· · · 0



∈R^M×n

mit einer oberen Dreiecksmatrix R ∈ R^n×n nennen wir QR-Zerlegung. F¨ur M =n gilt QΦ =R.

2 Q wird konstruktiv in n Schritten gebildet. Falls M = n ben¨otigt man n −1 Schritte. Bei diesem Rechengang wirdQ aber nicht explizit aufgestellt, weil dies den Aufwand nur unn¨otig vergr¨ossern w¨urde. Stattdessen konstruiert man Q als ein Produkt von n bzw. n−1 einfachen unit¨aren Matrizen.

Wir wenden uns zun¨achst dem Spezialfalln = 1 zu:

Sei Φ = (ϕ~1)∈R^M×1. Dann erreichen wir Q·Φ =















∗ 0 ... 0















mit ∗ = R∈R^1×1

| ∗ | = k~ϕ₁k₂ wenn wir Q als geeignete Spiegelung w¨ahlen. Q = I_n− _~uT²~u~u~u^T beschreibt eine Spiegelung an der Hyperebene H im R^M mit Normalenvektor~u:

Denn

~

x=λ~u: Q~x = (In− 2

~

u^T~u~u~u^T)~x

= ~x− 2

~

u^T~u~u(~u^T~x)

= ~x− 2

~

u^T~u(~u^T~x)~u

= λ~u− 2

~

u^T~u(~u^Tλ~u)~u

= λ~u−2λ~u

= −λ~u=−~x

~x⊥~u: ~u^T~x = 0

⇒Q~x = ~x Man bezeichnet diese Matrizen auch als

Householdermatrix:

U =I− 2

~

u^H~u~u~u^H

(benannt nach A.S. Householder, der sie zuerst in diesem Zusammenhang benutz-te.)

Zu einer gegebenen Spalte~xwollen wir nun solch eine Spiegelung, d.h.~u konstru-ieren, die diese in ein Vielfaches des ersten Koordinateneinheitsvektors ¨uberf¨uhrt.

Ist~xgegeben, so kann man ein solches ~usofort angeben: x= (x₁, . . . , x_n)^T

~u =











(|x1|+k~xk2)σ x₂

... x_n











Dabei istσ das verallgemeinerte Vorzeichen vonx1: σ= sign ₀(z)^def=

1 z = 0 z/|z| sonst.

Beispiel 4.14. Wir nehmen

~x = (−8,3,1,5,−1)^T . Dann leistet offenbar

~

u = (−18,3,1,5,−1)^T das Gew¨unschte, denn

~

u^T~u = 360, ~u^T~x = 180, also ~x− 2

~

u^T~u(~u^T~x)~u = (10,0,0,0,0)^T .

Diese Methode wird nun systematisch auf Φ angewendet: Die erste Transformation U₁ transformiert die erste Spalte von Φ auf ein Vielfaches des 1. Einheitsvektors. U₁ wird auf alle Spalten von Φ angewendet und auf den Vektor der Messwerte~y. Danach wird die gleiche Vorgehensweise wiederholt, jetzt mit den Komponenten 2, . . . , M der zweiten Spalte vonU₁Φ . Dies definiertU₂. AuchU₂ wird auf alle ¨ubrigen Spalten von U₁Φ angewendet und aufU₁~y usw.

Allgemein lautet der Algorithmus:

Hier werden im Schritt i die Spalten in zwei Teilspalten zerlegt: Die erste Teilspalte, gekennzeichnet durch die Doppel-Tilde, ¨andert sich nicht (die ersten i−1 Zeilen des Systems bleiben unge¨andert) und die zweite Teilspalte (Tilde) wird mit der Househol-dermatrix multipliziert, wobei deren Struktur explizit ausgenutzt wird.

Beispiel 4.15.

Durch ~u₁ wird U₁ gegeben: Schliesslich hat man

U_n. . . U₁(Φ~x−~y) =

worin ~c₁ die ersten n Komponenten der transformierten rechten Seite sind. Die L¨osung der Ausgleichsaufgabe bestimmt sich dann aus

R~x = ~c₁

und die L¨ange des optimalen Residuenvektors ist||~c₂||₂. Mit allgemeinen Bezeich-nungen haben wir

Satz 4.13. QR Zerlegung und Anwendung Es sei A ∈ R^m×n mit m ≥ n. Dann existiert eine orthonormale Matrix Q ∈ R^m×m mit QA =



 R

· · · 0



, R n×n obere Dreiecksmatrix. Ist A vom Rang n, dann ist R in-vertierbar und die Aufgabe:

Bestimme~x^∗ :

kA~x^∗ −~bk²₂ ≤ kA~x−~bk²₂ f¨ur alle ~x∈Rⁿ besitzt eine eindeutig bestimmte L¨osung ~x^∗, die sich aus

R~x^∗ =~c₁ errechnet, wo Q~b =





~c₁

· · ·

~c2



 mit~c₁ ∈Rⁿ.

(R ist in diesem Falle regul¨ar) 2

Beispiel 4.16.

(A,b)= [ -4 1 , 4.5 ] [ 2 2 , -1.0 ] [ 2 2 , 2.0 ] [ 1 1 , -1.5 ] ; u_1= [ -9 , 2 , 2 , 1 ]’ ;

U_1 = I - (2/(u_1’u_1))*u_1*u_1’ ;

U_1(A,b) = [ 5.0000 1.0000 , -3.5000 ] [ 0.0000 2.0000 , 0.7778 ] [ 0.0000 2.0000 , 3.7778 ] [ 0.0000 1.0000 , -0.6111 ] ; u_2= [ 0 , 5 , 2 , 1 ]’ ;

U_2= I - (2/(u_2’*u_2))*u_2*u_2’ ;

U_2*U_1*(A,b) = [ 5.0000 1.0000 , -3.5000 ] [ 0.0000 -3.0000 , -2.8333 ] [ 0.0000 0.0000 , 2.3333 ] [ 0.0000 0.0000 , -1.3333 ] ; x2= (-2.8333)/(-3.0000)= 0.9444 ;

x1= (-3.5000 - 0.9444)/5.0000 = -0.8889 ;

Residuenlaenge = 2.6874

2 NUMAWWW lineare Gleichungssyteme, QR-Zerlegung

Bemerkung: Die Methode der kleinsten Quadrate ist nat¨urlich nicht auf Messda-ten mit einem “freien” Parameter (im Beispielt) beschr¨ankt. Man kann sie w¨ ort-lich auch auf ganz allgemeine Ans¨atze

yi = a1φ1(ξi, ηi, . . .) +. . .+anφn(ξi, ηi, . . .) , i= 1, . . . N anwenden, wobeiξ, η, . . . die “Messstellen” repr¨asentieren.

4.7 Zusammenfassung

Das Standardverfahren zur L¨osung linearer Gleichunggsysteme ist der Gauss’sche Algorithmus. Um den Einfluss von Rundungsfehlern auf die berechnete L¨osung unter Kontrolle zu halten, ist die Anwendung von Pivotisierungsregeln unerl¨asslich, mit Ausnahme spezieller Matrizen, insbesondere der hermitisch positiv definiten.

Dieser Algorithmus erzeugt eine Faktorisierung P A = LR bzw.

P AQ = LR

mit Permutationsmatrizen P und Q (die durch Permutationsvektoren repr¨ asen-tiert werden) und einer unteren Dreiecksmatrix L mit Diagonale (1, . . . ,1), engl

”unit lower triangular”, und einer oberen Dreiecksmatrix R. Auf der Diagonalen von R stehen dann die Pivotelemente. Es ist daher

det(A) = ±detR = ±

n

Y

i=1

ρ_i,i

Diese Zerlegung ersetzt die Information ¨uber A gleichwertig und erlaubt z.B.

die sp¨atere L¨osung von Gleichungssystemen mit A bei beliebiger rechter Seite b.

Ist A hermitisch und positiv definit, dann kann man zweckm¨assig die Cholesky-Zerlegung

A = LL^H

mit einer unteren DreiecksmatrixLmit positiven reellen Diagonalelementen ver-wenden. Der Gauss’sche Algorithmus erlaubt die Ber¨ucksichtigung von Besetzt-heitsstrukturen (Bandstruktur, Hessenbergstruktur, auch ”sparsity”). Der Fehle-reinfluss bei der Anwendung dieses Algorithmus oder allgemeiner von Datenfeh-lern in Matrix und Inhomogenit¨at wird beschrieben durch die sogenannte ”Kondi-tionszahl” der Matrix. Wir haben (etwas vergr¨obert) die Aussage ”normrelativer

Fehler in der L¨osung kleinergleich Summe der normrelativen Fehler in Matrix und rechter Seite, multipliziert mit der Konditionszahl”. Die Konditionszahl ist stets gr¨ossergleich eins und oft sehr gross gegen eins. Lineare Ausgleichsaufgaben kann man ¨uber die Normalgleichungen mit Hilfe der Choleskyzerlegung l¨osen. Wegen des u.U. sehr verst¨arkten Fehlereinflusses sollte man aber besser den Weg ¨uber die QR-Zerlegung der Ansatzmatrix gehen. Die QR-Zerlegung vermittelt zugleich eine Berechnung von Orthogonalbasen von Bildraum R(A) und Kern N(A^H).

Im Dokument 4 Lösung linearer Gleichungssysteme: Direkte Methoden (Seite 32-43)