1. Phase:Methode der kleinsten Quadrate

Einführung

Im Vortrag über das CT-Verfahren hat Herr Köckler schon auf die Methode der kleinstenQuadrate hingewiesen. Diese Lösungsmethode, welche bei überbestimmten Gleichungssystemen(d.h. mehr Gleichungen als Variablen) zur Anwendung kommt, werden wir als erstes näherkennenlernen.

�MatrizenA, b

�

�

�

�

Lösen mit der

KQ-Methode:

x = lsqr(A,b)

�Lösungsvektorx

Bevor wir allerdings die Rekonstruktion eines Objektes mittels der Methode der kleinsten Qua-drate anwenden, machen wir uns einige Gedanken zur Vorverarbeitung und Vereinfachung.Hierbei werden wir auch gleich auf Matlab zurückgreifen.

1 Vorbereitung

Bestimmen Sie eine Gerade, welche durch die Punkte P (2, 1) und Q(6, 9) verläuft.

Hinweise:

Wie lautet die allgemeine Geradengleichung?Wo werden die Koordinaten der Punkte eingetragen?Welches sind die Unbekannten?

Lösung:Die allgemeine Geradengleichung lautet y = mx + n. In diese werden jeweils die x− und die y-Koordinate eingesetzt. Die Unbekannten sind m und n. Für diese erhält man ein lineares Gleichungs-system:

I. 1 = 2m + n

II. 9 = 6m + n II.−I.=⇒I. 1 = 2m + n

II. 8 = 4m II./4=⇒

I. 1 = 2m + n

II. 2 = m

I.−2II.=⇒ I. −3 = n

II. 2 = m

In Matrix-Vektor-Notation sieht das wie folgt aus:(19

)=

(2 16 1

)(mn

)⇒

(18

)=

(2 14 0

)(mn

)⇒

(12

)=

(2 11 0

)(mn

)⇒

(−32

)=

(0 11 0

)(mn

)

Daher lautet die Geradengleichung y = 2x − 3.

1

2 KQ-Methode mit einer Unbekannten

Das Gewicht G einer Probe wird viermal mit verschiedenen Instrumenten und Methoden gemessen,das ergibt die Messungen G1, G2, G3 und G4.

Eine möglichst gute Nährung für das Gewicht soll jetzt nach der Methode der kleinsten Quadratebestimmt werden, d.h. es soll der Ausdruck (G−G1)2 + (G−G2)2 + (G−G3)2 + (G−G4)2 bezüglichG minimiert werden.

Beschreiben Sie, wie man auf den genannten Ausdruck kommt.Bestimmen Sie die Näherung für G.

Bemerkung: Man kann das Verfahren der Ableitung - wie im Vortrag erlebt - auf mehrere Variablenausdehnen. In diesem Fall erhält man ein Gleichungssystem, welches zur Bestimmung der einzelnenVariablen gelöst wird.

Hinweise:

Aufstellen der Gleichungen und anschließendes Umstellen.Bilden der Ableitung!Nullstellen der Ableitung als mögliche Extremstellen.

Lösung:Das überbestimmte Gleichungssystem lauten zunächst

G = G1

G = G2

G = G3

G = G4

bzw. in Matrix-Vektor-Notation

⎛⎜⎜⎝

1111

⎞⎟⎟⎠ G =

⎛⎜⎜⎝

G1

G2

G3

G4

⎞⎟⎟⎠

Es ist dann f(G) = (G−G1)2 + (G−G2)2 + (G−G3)2 + (G−G4)2. Dann erhält man durch Ableitenf ′(G) = 2(G − G1) + 2(G − G2) + 2(G − G3) + 2(G − G4).Nullsetzen und nach G auflösen ergibt somit:

8G − 2G1 − 2G2 − 2G3 − 2G4 = 0 =⇒ 8G = 2(G1 + G2 + G3 + G4) =⇒ G = 14(G1 + G2 + G3 + G4)

Hätte man das Ergebnis einfacher errechnen können?Welcher bekannte Begriff liefert diese einfachere Rechnung?

Lösung:Ja.Der Mittelwert der Messwerte ergibt gerade das Minimum der quadratischen Abweichungen.

2

3 KQ-Methode in der CT -ein einführendes Beispiel

Wir schauen uns nun ein erstes - einfaches - Beispiel an, welches schon dem Bereich der Computertomo-graphie zugeordnet werden kann. Hierzu sei die folgende Anordnung gegeben:

9

9

9

9

9

9

4 2

2

5

7

8

Das Ziel dieser Aufgabe ist, ein Gleichungssystem aufzustellen, mit welchem die beste Näherung für dieroten und blauen Kästchen berechnet werden kann.

Erinnern wir uns an den Vortrag von Herrn Köckler. Dort wurde die Formel zum Aufstellen des Glei-chungssystemes eingeführt. Für einen Strahl lautet die allgemeine Gleichung:

L1µ1 + L2µ2 + . . . + Lnµn = − ln(

I1

I0

)

mit I0 als gegebener Eingangsintensität und I1 als gemessener Ausgangsintensität. Weiter ist n dieAnzahl der durchstrahlten Zellen, welche voneinander verschieden sind. In unserem Beispiel gibt esnur rote und blaue Zellen, somit ist n = 2. Da wir sechs parallele Strahlen im Beispiel sehen undpro Strahl eine Gleichung aufgestellt wird, erhalten wir am Ende ein Gleichungssystem bestehendaus sechs Gleichungen der obigen Form. D.h. eine Gleichung des Gleichungssystems hat also dieForm

L1µ1 + L2µ2 = − ln(

I1

I0

).

Die Werte L1 und L2 beinhalten die berechneten Längen des Strahles in der jeweiligen Zelle – in un-serem Beispiel in den roten bzw. blauen Kästchen. Wird eine Zelle nicht getroffen, erhält der Lj-Wertfür diese Zelle den Wert Null. Die µ-Werte sind die unbekannten Abschwächungskoeffizienten, die wirbestimmen wollen.

3

Aufgabenteil 1:In diesem Aufgabenteil machen wir uns die Rechnung etwas einfacher und setzen Lj entweder gleich 1,falls eine Zelle getroffen wurde, oder gleich 0, falls eine Zelle nicht getroffen wurde. Als Zelle bezeichnenwir dabei ein Kästchen, wie im Bild zu sehen.

Stellen Sie das Gleichungssystem für diesen vereinfachten Fall auf.Lösung:

µrot = − ln(

49

)

3µrot = − ln(

29

)

2µrot + µblau = − ln(

29

)

µrot + 2µblau = − ln(

59

)

3µblau = − ln(

79

)

µblau = − ln(

89

)

In Matrizenschreibweise Aµ = b mit µ = (µrot, µblau) erhält man⎛⎜⎜⎜⎜⎜⎜⎝

1 03 02 11 20 30 1

⎞⎟⎟⎟⎟⎟⎟⎠

·(

µrot

µblau

)= −

⎛⎜⎜⎜⎜⎜⎜⎝

ln 49

ln 29

ln 29

ln 59

ln 79

ln 89

⎞⎟⎟⎟⎟⎟⎟⎠

Ist eine eindeutige Lösung des Systems möglich?Lösung:Nein.

Nun wollen wir das System analog zur Aufgabe 2 für µblau = 0.11 lösen. Welches Ergebnis für µrot

erhalten Sie?

Lösung:Sei µblau = 0.11.Zu lösen ist

min(

(µrot + ln49)2 + (3µrot + ln

29)2 + (2µrot + µblau + ln

29)2+

+(µrot + 2µblau + ln59)2 + (3µblau + ln

79)2 + (µblau + ln

89)2

)

4

Ableiten liefert

0 = 2(µrot + ln49) + 6(3µrot + ln

29) + 4(2µrot + µblau + ln

29) + 2(µrot + 2µblau + ln

59).

Umstellen ergibt 30µrot = 16.9582 und somit ist (gerundet)

µrot = 0.5653.

Als nächstes wollen wir nun das Gleichungssystem mit beiden Unbekannten bestmöglichst lösen. Hierzunutzen wir die kleinste Quadrate Lösungsmethode lsqr, welche in Matlab implementiert ist.

Einige Vorbereitungen:Das oben aufgestellte Gleichungssystem muss in Matrix-Vektor-Schreibweise geschrieben werden, d.h.es sollte die Form Aµ = b mit µ = (µrot, µblau) besitzen. In der Matrix A stehen dabei die Koeffizientenvon µrot und µblau und in b die Einträge der rechten Seite des Gleichungssystems. Sowohl A als auchb müssen dann in Matlab über die Befehle MatrixEingabe bzw. VektorEingabe eingegeben werden(siehe Übersicht der Matlab-Befehle).

Ist dies geschehen, können Sie mit Hilfe des Befehls mu=lsqr(A,b) das Gleichungssystem mit derMethode der kleinsten Quadrate lösen. Wie ist das Ergebnis?

Eingabe der Werte:

A = MatrixEingabe(2, 6, [1 0 3 0 2 1 1 2 0 3 0 1]);b = VektorEingabe(6, [-log(4/9) -log(2/9) -log(2/9) -log(5/9) -log(7/9) -log(8/9)]);

Lösen:

mu = lsqr(A,b); liefert µrot = 0.5722 und µblau = 0.0842.

Um die Ergebnisse auf ihre Genauigkeit hin überprüfen zu können, schaut man sich die sogenannte Feh-lerquadratsumme oder auch Residuum genannt etwas näher an. Hierzu berechnet man

‖b − Aµ‖2 =

√√√√ n∑i=1

(bi − (Ax)i)2.

Dieser Ausdruck kann per Hand berechnet werden, aber um uns Zeit und Rechnerei zu ersparen, nutzenwir direkt wieder Matlab mit dem Befehl

fehler = norm(b - A*mu)

Berechnen Sie mit diesem Befehl zunächst den Fehler für den Fall, in welchem µblau vorher auf 0.11gesetzt wurde. Anschließend wollen wir diesen Fehler mit dem vergleichen, der beim Lösen des Glei-chungssystems mit Lj = {0, 1} entstanden ist. D.h. wir setzen für mu die Werte ein, welche vonMatlab mittels lsqr berechnet wurden.

Residuum einfacher Fall:µblau = 0.11, µrot = 0.5653

5

fehler1 = norm(b - A*mu) = 0.4602

Residuum für Fall Lj = {0, 1}:µrot = 0.5722 und µblau = 0.0842

fehler2 = norm(b - A*mu) = 0.4500

Aufgabenteil 2:Das Ziel ist es nun, die eben berechneten Fehler zu verbessern, indem wir die zu Beginn des Aufgabentei-les 1 gemachte Einschränkung nun durch bessere (exaktere) Werte für die Lj ersetzen.

Stellen Sie erneut ein Gleichungssystem auf. Die Werte für die verschiedenen Lj entnehmen Sie direktdem Bild, indem Sie die Länge der verschiedenen Strahlen messen.

Lösung:

1.1µrot = − ln(

49

)

2.3µrot = − ln(

29

)

2.4µrot + 0.2µblau = − ln(

29

)

0.5µrot + 2.1µblau = − ln(

59

)

2.6µblau = − ln(

79

)

1.1µblau = − ln(

89

)

In Matrizenschreibweise Aµ = b mit µ = (µrot, µblau) erhält man⎛⎜⎜⎜⎜⎜⎜⎝

1.1 02.3 02.4 0.20.5 2.10 2.60 1.1

⎞⎟⎟⎟⎟⎟⎟⎠

·(

µrot

µblau

)= −

⎛⎜⎜⎜⎜⎜⎜⎝

ln 49

ln 29

ln 29

ln 59

ln 79

ln 89

⎞⎟⎟⎟⎟⎟⎟⎠

Um nun die Werte für µrot und µblau erneut berechnen zu können, muss die Matrix A neu in Matlabeingeben werden. Ist es notwendig den Vektor b noch einmal einzugeben?Lösung:Nein.

Nachdem Sie die Matrix eingegeben haben, berechnen Sie bitte wie zuvor mit Matlab die Lösung mitder Methode der kleinsten Quadrate.

6

Eingabe der Werte:

Aneu = MatrixEingabe(2, 6, [1.1 0 2.3 0 2.4 0.2 0.5 2.1 0 2.6 0 1.1]);

Lösen:

muneu = lsqr(Aneu,b); liefert µrot = 0.6468 und µblau = 0.1070.

Wie bereits angekündigt, wollen wir dieses Ergebnis mit den vorherigen vergleichen. Dazu berechnenSie bitte erneut die Fehlerquadratsumme. Was stellen Sie fest?

Residuum für die verbesserten Eingabewerte:µrot = 0.6468, µblau = 0.1070

fehler3 = norm(b - Aneu*muneu) = 0.1316

Es fällt auf, dass dieses zum Schluss berechnete Residuum das kleinste ist. Dies kommt daher, dass indiesem Fall die bestmöglichsten Eingabewerte benutzt worden.

7