5 F¨ unfte ¨ Ubungseinheit
Inhalt der fünften Übungseinheit:
• Nichtlineare Datenmodelle
• Kurze Wiederholung: lineare und nichtlineare Datenmodelle
• Bonus-Material: MATLAB-Werkzeuge zum Anpassen von Funktionen an Daten
• Alternativen zur Minimierung der Fehlerquadrate: Robuste Regression,
• Anwendungen der Singulärwert-Zerlegung: Total Least Squares, Da- tenkomprimierung
5.1 Uberbestimmte nichtlineare Systeme, Gauß-Newton-Verfahren ¨
Das Gauß-Newton-Verfahren findet Näherungslösungen (im Sinn der kleinsten Fehlerquadrate) für überbestimmte nichtlineare Systeme. Die Grundidee ist, ähnlich wie bei der Lösung nichtli- nearer Gleichungssysteme mit dem Newton-Verfahren, mit derJacobimatrix des nichtlinearen Systems iterativ Korrekturterme zu berechnen. Die Aufgaben 41 , 42 und die Vorlesungsfolien der 5. Vorlesung verwenden dieses Verfahren und zeigen den Rechenweg.
Aufgabe 41: Standortbestimmung durch Trilateration
Die Abstände von drei festen PunktenA, B, C in derxy-Ebene zu einem unbekannten Punkt X sind (etwas ungenau) bekannt. Gesucht ist eine möglichst gute Positionsbestimmung.
Punkt x y Entfernung
1 1 6
8 4 3,6
5 8 4,2
Die entsprechenden Gleichungen lauten:
p(x1−1)2+ (x2−1)2 = 6 p(x1−8)2+ (x2−4)2 = 3.6 p(x1−5)2+ (x2−8)2 = 4.2
Den drei Gleichungen entsprechen drei Kreise imR2. Sie haben keinen gemeinsamen Schnitt- punkt. Dementsprechend gibt es keine exak- te Lösung des überbestimmten Gleichungssys-
tems. 00 1 2 3 4 5 6 7 8 9
1 2 3 4 5 6 7 8 9
A
B C
da=6
dc=4.2
db=3.6
X
Die Folien der 5. Vorlesung diskutieren dieses Modell und rechnen die ersten Schritte vor.
Schreiben Sie ein MATLAB-Programm, das die kleinste-Quadrate-Anpassung für den Standort findet.
Bemerkung: Das GPS-System misst zur Positionsbestimmung die Abstände (genauer: Signal- laufzeiten) zu mehreren Satelliten und bestimmt auch durch Ausgleichsrechnung die Positi- on.
Aufgabe 42: Gauß-Newton-Verfahren in Wikipedia
Die Angabe zum folgenden Beispiel stammt aus der englischen bzw. französische Wikipedia (Stichworte Gauss-Newton algorithm bzw.Algorithme de Gauss-Newton). Lösen Sie die Auf- gabe mit folgender Anleitung in MATLAB.17
Biologie-Experimente zur Beziehung zwischen Substanz-Konzentrationxund Reaktionsratey in einer durch Enzyme vermittelten Reaktion haben die Daten in folgender Tabelle ergeben:
x | 0.038 0.194 0.425 0.626 1.253 2.500 3.740 y | 0.05 0.127 0.094 0.2122 0.2729 0.2665 0.3317
Gesucht ist eine Kurve (Modellfunktion) der Form y=a x
b+x
die im Sinn der kleinsten Quadrate die Daten am besten approximiert. Die Parameteraundbsind zu bestimmen.
Einsetzen der Daten ergibt sieben nichtlineare Glei-
chungen in den beiden Unbekannten a, b. 00 1 2 3 4
0.05 0.1 0.15 0.2 0.25 0.3 0.35
[S]
reaction rate
a 0.038
b+ 0.038 −0.05 = 0 a 0.194
b+ 0.194−0.127 = 0 ...
a 3.740
b+ 3.740 −0.3317 = 0
Das Gleichungssystem in Vektor-Schreibweise:
f(x) = 0 mitf :R2→R7
Die Jacobi-Matrix Df dieses Systems ist eine 7×2-Matrix. Zeile i enthält die partiellen Ab- leitungen deri-ten Gleichung nach den Unbekanntenaundb.
(Df)i1= xi
b+xi, (Df)i2=− axi
(b+xi)2
Der weitere Rechenweg verläuft völlig analog zu Aufgabe 21 und der Musterlösung auf Seite 22 im Skriptum: Wählen Sie als Startwert [a;b]=[0.9;0.2], werten Sie f und Df aus. Die kleinste-Quadrate-Näherung aus dem überbestimmten Gleichungssystem Df∆x=−f liefert den Korrekturvektor ∆x.
17Wer dieses Beispiel in der englischen oder französischen Wikipedia-Seite durchliest und bei der Abbildung der Modellkurve auf “more details” klickt, findet fertigen MATLAB-Code zur Lösung des Beispiels und Zeichnen der Kurve. Aber das Verstehen von fremdem Code ist auch nicht einfach.
Beachten Sie: bei einer quadratischen Jacobi-Matrix liefert der MATLAB-BefehlDf(x0)\f(x0) die Lösung des Gleichungssystems; hat die Jacobi-Matrix mehr Zeilen als Spalten (überbe- stimmtes System), liefert derselbe Befehl nicht die Lösung (es gibt keine exakte Lösung!), sondern eine Anpassung mit geringstmöglichem Restfehler (die „am wenigsten falsche Ant- wort“)!
Bonus-Frage: Unser Newton-Musterprogramm findet die verbesserte Näherung mit dem MATLAB- Befehl
x = x0 - Df(x0)\f(x0); % Newton-Schritt
Das Wikipedia-Beispiel nennt die Jacobi-Matrix Jf und verwendet für denselben Schritt (in MATLAB-Code geschrieben) den Befehl
x = x0 - (J’*J)\(J’*f(x0));
Warum einerseitsDf(x0)\f(x0) und andererseits(J’*J)\(J’*f(x0))? Was wird da jeweils gerechnet? Vor- und Nachteile der beiden Varianten im Vergleich?
Aufgabe 43: Kalorimeterversuch
Der DatensatzKalorimeter.datauf der Übungsseite enthält Wertepaare (Zeittin Minuten, Temperatur T in Celsius) zu einem Kalorimeterversuch, gemessen ab der Einbringung eines Versuchsobjektes. Die Temperatur sinkt zuerst rasch, steigt dann aber wieder und gleicht sich langsam der Umgebungstemperatur an.
Der Verlauf der TemperaturT als Funktion der Zeittsoll durch ein Modell der Form T(t) =T0+C1exp(−λ1t) +C2exp(−λ2t)
beschrieben werden. Bestimmen Sie die Parameter dieses nichtlinearen Modells ausgehend von den Startwerten
T0= 6 C1=−1 λ1= 1/10 C2= 20 λ2= 2
(drei Iterationen des Newton-Verfahrens reichen aus.) Zeichnen Sie Messpunkte und Modell- funktion.
5.2 Kurze Wiederholung: lineare und nichtlineare Datenmodelle
Die Aufgaben zur Anpassung von Funktionen an Daten beschränkten sich in der vorigen Ein- heit auf lineare Modelle. Dabei bedeutet „linear“: die gesuchte Anpassungf ist eine Linear- kombination vonnBasisfunktionenφ1, φ2, . . . , φn in der Form
f =a1φ1+a2φ2+· · ·+anφn .
Wichtig: Die Basisfunktionen φ1, . . . , φn können beliebige, auch nichtlineare Funktionen sein – es sind die gesuchten Koeffizientena1, . . . , an, die nur linear im Ansatz auftreten!
Verwirrender Weise wird für die Anpassung einer Ausgleichs-Geraden auch der Begriff „lineare Regression“ verwendet. Unter „polynomiale Regression“ versteht man die Anpassung eines Polynoms an gegebene Daten, wobei es sich aber um ein lineares Datenmodell handelt.
Unterscheiden Sie:
• Typ der Anpassungsfunktionen: Gerade ("lineare Regression"), Polynom ("polynomiale Regression"), Winkel- oder Exponentialfunktionen. . .
• Verknüpfung der Ansatzfunktionen: linear, nichtlinear. Dazu gehören die Begriffe lineares/nichtlineares Datenmodell.
Vergleichen Sie: Lineare Modelle
Aufgabe 38: Ansatzz(x, y) =a1+a2x+a3x2+a4y.
Hier sind die Basisfunktionenφ1(x, y) = 1, φ2(x, y) =x, φ3(x, y) =x2, φ4(x, y) =y. Aufgabe 39: Ansatzt(d, T0) =a1+a2d+a3T0+a4dT0.
Basisfunktionen sindφ1(d, T0) = 1, φ2(d, T0) =d, φ3(d, T0) =T0, φ4(d, T0) =d·T0. Aufgabe 40: Ansatzy=a1+a2cos
xπ 6
+a3sin xπ
6 Basisfunktionen sindφ1(x) = 1, φ2(x) = cos xπ6
, φ3(x) = sin xπ6. Nichtlineare Modelle
Aufgabe 42: Ansatzy=ab+xx . Hier sind die gesuchten Koeffizientenaundb nichtlinear verknüpft!
Aufgabe 43: AnsatzT(t) =T0+C1exp(−λ1t) +C2exp(−λ2t).
Wärenλ1undλ2 bekannt, dann wäre es ein linearen Modell mit ParameternT0, C1und C2. Weil auch die beidenλi gesucht sind, wird die Aufgabe nichtlinear.
Kapitel 5.1 in diesen Unterlagen zeigt einen Standard-Lösungsweg zur nichtlinearen Anpas- sung. Der folgende Abschnitt 5.3 zeigt weitere Möglichkeiten und stellt MATLAB-Werkzeuge vor. Vielleicht finden Sie diese Werkzeuge einfacher und intuitiver zu verwenden als den Lösungsweg aus Kapitel 5.1.
5.3 MATLAB-Werkzeuge zum Anpassen von Funktionen an Daten
(Bonus-Material für Interessierte: Im Vergleich zu älteren Versionen bietet MATLAB inzwi- schen tolle neue Werkzeuge. Schaut euch ’s an, es zahlt sich aus!)
Zur MATLAB-Grundausstattung gehört dasBasic Fitting Tool. Damit sollten Sie unbedingt umgehen können. Machen Sie sich gleich einmal damit vertraut: Finden Sie die entsprechende Anleitung in der MATLAB-Hilfe unter MATLAB Mathematics Elementary Math Descriptive Statistics
Examples and How To Interactive Fitting oder suchen Sie nach dem Stichwort Interactive Fit- ting.
Sie finden dort ein durchgearbeitetes Bei- spiel, das die Verwendung des Basic Fitting Tools für interaktives Anpassen von Kurven an Datenpunkte erklärt.
Arbeiten das Beispiel in der MATLAB-Hilfe durch und stellen Sie die Approximation gra- phisch dar, etwa so wie in der nebenstehen- den Abbildung.
Darüber hinaus bietet MATLAB in der Sta- tistics and Machine Learning Toolbox und in derCurve Fitting Toolbox viele Werkzeu- ge zur Datenanpassung, weit mehr als wir in diesen Übungen behandeln können. Die Aufgaben 44 und 46 zeigen beispielhaft, wel- che Möglichkeiten die Befehle cftool und polytool bieten. Aufgabe 47 erwähnt den Befehlrobustfit.
Aufgabe 44: Basic Fitting Tool im Vergleich zu Curve Fitting Toolbox
Hier lösen Sie das nichtlineare Ausgleichsproblem aus Aufgabe 42 mit MATLAB-Toolboxen.
Legen Sie mit MATLABsbasic fitting tool eine Ausgleichsgerade, ein quadratisches und dann ein kubisches Ausgleichspolynom durch folgende Datenpunkte und plotten Sie die Ergebnis- se.
X = [0.038 0.194 0.425 0.626 1.253 2.500 3.740];
Y = [0.05 0.127 0.094 0.2122 0.2729 0.2665 0.3317];
Es handelt sich hier um die Daten aus Aufgabe 42; dort finden Sie auch eine Abbildung mit gut angepasster Kurve. Welche anderen Optionen des basic fitting tool liefern plausible Kurven?
Was bedeuten die Fehlermeldungen, die Sie ab Grad 7 bekommen?
Es stellt sich heraus: Kein Ausgleichspolynom kann die Datenpunkte befriedigend approximie- ren. Das liegt teils daran, dass die Datenpunkte aufgrund der Messunsicherheit weit gestreut sind, aber auch daran, dass sich im gemessenen Experiment die y-Werte mit zunehmendem x asymptotisch einem konstanten Wert nähern. Kein Polynom kann so ein asymptotisches Verhalten beschreiben.
Der Befehl cftool öffnet die Curve Fitting App. Hier können Sie neben Polynomen auch weitere Modellfunktionen anpassen. Versu- chen Sie es: Wählen Sie links oben unter X DataundY Datadie entsprechenden Daten- vektoren aus den im Workspace vorhandenen Variablen. Oben Mitte können Sie den Funk- tionstyp wählen. Probieren Sie zuerst Poly- nomial mit verschiedenen Polynom-Graden.
Wählen Sie anschließend Custom Equation und geben Sie den Funktionsterna*x/(b+x) ein. (Das ist die Modellfunktion aus Aufga- be 42!) Ihre Anpassung sollte so aussehen wie hier gezeigt.
Zusätzlich zu den berechneten Werten der Koeffizienten aund bliefert MATLAB auch ein Konfidenzintervall. Interpretation: Falls den Daten tatsächlich ein Modell der Form y = ax/(b+x) zugrunde liegt, aber die y-
Daten zufallsbedingt verrauscht sind, dann berechnet MATLAB Schätzungen ˆaund ˆb für die Parameterabeziehungsweiseb. Zum Beispiel:
Coefficients (with 95% confidence bounds):
a = 0.3618 (0.2363, 0.4874) b = 0.5563 (-0.05629, 1.169)
Die Interpretation „Die tatsächlichen Werte vonaundbliegen mit 95%-iger Sicherheit im be- rechneten Intervall“ ist so nicht korrekt; die tatsächlichen Werte von aundbliegen entweder drinnen oder nicht – da ist kein Spielraum für Wahrscheinlichkeiten. Es ist umgekehrt: die be- rechneten Grenzen des Konfidenzintervalls sind unsicher. Jede unabhängige Wiederholung der Messung liefert (für die immer gleichen Modellparameter aundb) einen neuen Satz von Da- tenpunkten mit anderen zufallsbedingten Fehlern. MATLABs Schätz-Methode berechnet dazu Konfidenzintervalle mit jeweils etwas anderen Grenzen. Meistens (in 95% der Fälle) schätzt Matlab die Fehlergrenzen richtig, aber in 5% der Fälle überdeckt das berechnete Konfidenzin- tervall die wahren Werte nicht.
Aufgabe 45: Der Befehl regress aus der Statistics Toolbox
Die Matlab-Hilfe Zum Stichwort „regress“ zeigt ein Beispiel mit Daten zu Gewicht, Motorleis- tung und Benzinverbrauch von Autos und fin- det dazu eine Anpassungs-Funktion. Verglei- chen Sie mit Aufgabe 39: dort wird ein Daten- modell derselben Art berechnet, die Form der Matrix ist gleich.
Sie können in der MATLAB-Hilfe aufOpen Li- ve Script klicken und gelangen so in den Live Editor, der mehr Möglichkeiten bietet als der Standard-Editor. Die Dateiendung*.mlskenn- zeichnet Live-Skript-Dateien.
Wenn Sie neugierig sind, probieren Sie den Live-Editor aus. Wenn Sie in der gewohnten Arbeitsumgebung bleiben wollen, speichern Sie dieses Musterbeispiel mittelssave asals Datei- typMATLAB Code Files, Dateiendung*.m Orientieren Sie sich an diesem Muster und erstellen Sie für die Daten aus Aufgabe 39 (Durch- messer, Ausgangstemperatur, Kochzeit) mit dem regress-Befehl eine Anpassung und eine ähnliche Graphik wie oben (mit Kochzeit als z-Achse).
Aufgabe 46: Der Befehl polytool aus der Statistics toolbox
Die folgenden MATLAB-Befehlszeilen erzeugen entlang der Funktiony= 2−3x+ 2x2Daten- punkte, die durch zufällige, normalverteilte Störungen verrauscht sind.
n = 2 0 ;
x =l i n s p a c e(0 ,1 , n );
y =2 - 3* x + 2* x .^2 + 0 . 1 *r a n d n(1 , n );
Der Befehlpolytool(x,y)öffnet ein Fenster ähnlich dem basic fitting tool. Erzeugen Sie damit ein Bild wie die nebenstehende Abbildung und erklären Sie:
• Was bedeuten die verschiedenen Kurven?
• Wozu dient das verschiebbare Achsenkreuz, was sind die angezeigten X-Values und Y-Values?
• Wenn Sie auf Export. . . klicken, was bedeuten
„Parameters“ und Parameters CI“?
5.4 Lineare Regression, Ausgleichsgerade, -ebene: weitere Methoden
Aufgabe 47: Robuste Regression
Lineare Regression nach der Methode der kleinsten Quadrate reagiert empfindlich auf Aus- reißer in den Datenpunkten. Die folgenden MATLAB-Befehle erzeugen einen Datensatz von Punkten, die–bis auf einen Ausreißer–annähernd linearen Trend zeigen.
x = r a n d(4 0 , 1) ;
y = x + 0 . 2 *r a n d n( 40 ,1 );
x (1) = -1; y ( 1 ) = 1 ;
Die Funktion randn erzeugt normalverteilte Abweichungen!
Schreiben Sie eine Script-Datei, die einen solchen Datensatz erzeugt und führen Sie folgende Schritte durch:
1. Stellen sie die Datenpunkte mit dem Symbol’*’in einem Diagramm dar.
2. Kapitel 6.3 des Vorlesungsskriptums beschreibt, wie man (nach der Methode der kleins- ten Quadrate) die Parameter a und b der Ausgleichsgeraden y = a+bx bestimmt.
Implementieren Sie diese Formeln in Ihrem Skriptfile und berechnen Sie damitaundb.
3. Zeichnen Sie mit den Werkzeugen des Menues „Tools–Basic Fitting“ eine Ausgleichsge- rade ein. Stellen Sie die Formel im Diagramm dar und vergleichen Sie mit den vorher berechneten Wertenaundb.
4. Minimieren der Fehlerquadrateist das Standardverfahren zur Regression. Warum macht es keinen Sinn, dieSumme der Fehler zu minimieren? Kapitel 6.6 des Vorlesungsskrip- tums beschreibt ein Regressionsverfahren, das die Summe der Fehlerbeträge minimiert.
Dazu gibt es ein Musterprogramm als Dateilinregrob.m. Zeichnen Sie die so bestimmte Gerade im Diagramm ein.
Ihr Diagramm sollte etwa so aussehen:
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−1.5
−1
−0.5 0 0.5 1 1.5
y = 0.52*x + 0.19
Datenpunkte Ausgleich nach kl. Quadr.
Ausgl. mit minim. Fehlerbetragssumme
Übrigens: Die MATLAB-Hilfe zum Befehlrobustfitaus der Statistics Toolbox zeigt als Bei- spiel einen ganz ähnlichen Fall wie die Aufgabe 47. Der Befehlrobustfitbietet im Vergleich zur Datei linregrob.mmehr Optionen. Wenn es Sie interessiert, können Sie die Aufgabe 47 auch mit robustfitlösen und die Ergebnisse mit den linregrob.m-Resultaten vergleichen.
Auch cftool aus der Curve Fitting Toolbox, siehe Aufgabe 44 beherrscht robustes Fitting:
wählen Sie dortPolynomial, Degree 1, Robust Bisquare oderLAR. Aufgabe 48: Total Least Squares
Das Standardverfahren zur linearen Regression minimiert Summe der Abstandsquadratein y- Richtung. Wenn in den Daten nur diey-Werte mit zufälligen, normalverteilten Fehlern behaftet sind, liefert dieses Verfahren statistisch gesehen eine optimale Schätzung für die Ausgleichsge- rade.
Wenn aber auch diex-Werte durch Fehler gestört sind (Einfachster Fall: zufällig, unabhängig, normalverteilt, Mittelwert 0, gleiche Standardabweichung wie diey-Werte), ist es sinnvoller, die Quadratsumme der tatsächlichen geometrischen Abstände zwischen Gerade und Datenpunkten zu minimieren.
DieTotal Least Squares-(TLS-) Methode minimiert Quadratsumme der Normalabstände zwi- schen Gerade und Datenpunkten. Sie verwendet dazu die Singulärwertzerlegung.
Vorgangsweise für Datenx= [x1x2. . . xn], y= [y1y2, . . . yn]:
• BestimmeSchwerpunkt [¯x,y¯] der Daten.
¯ x= 1
n X
i=1,n
xi, y¯= 1 n
X
i=1,n
yi
• Verschiebe die Daten
∆xi=xi−x,¯ ∆yi=yi−y¯
• Bilde Singulärwertzerlegung
U·S·VT =
∆x1 ∆y1 ... ...
∆xn ∆yn
Der entsprechende MATLAB-Befehl lautet
> > [ U S V ]=svd([ d e l x d e l y ]) Eingabe delx und dely als Spaltenvektoren!
• TLS-Gerade geht durch den Schwerpunkt [¯x,y¯] in Richtung des ersten Spaltenvektors vonV und hat als Normalvektor den zweiten Spaltenvektor vonV. Die Normalform der Geradengleichung lautetax+by=c. Dabei ist
a b
der zweite Spaltenvektor vonV und c=ax¯+by¯ergibt sich durch Einsetzen des Schwerpunktes in die Geradengleichung.
Die Matlab-Befehle
n = 2 0 ;
x = l i n s p a c e( -1 ,1 , n ) ’;
y = 3* x +2;
y = y + 0 . 2 *r a n d n( n , 1 ) ; x = x + 0 . 2 *r a n d n( n , 1 ) ;
erzeugen Punkte auf der Geradeny= 3x+2 mit Fehlern in x- und y-Richtung. Finden Sie die Ausgleichsgerade mit dem klassischen Ansatz und mit TLS.
Wiederholen Sie die Berechnung einige Male (mindestens 10 mal). Sie testen also für immer denselben Zusammen- hangy= 3x+ 2 unterschiedliche verrauschte Datensätze und schätzen aus den fehlerbehafteten Daten die Werte k undd iny =kx+d. Welche Methode liefert im Mit- tel die besseren Schätzungen für Steigungk(hier 3) und konstanten Termd(hier 2) der Geraden?
Ihre Berechnungen sollten zeigen: Schätzungen mit klas- sischer kleinste-Quadrate-Regression liefern im Mittel zu kleine Werte fürk, währed die TLS-Schätzungen im Mit- tel nahe am korrekten Wert liegen. Man sagt: der TLS- Schätzer isterwartungstreu, der andere Schätzer istnicht erwartungstreu. Die Wahl der korrekten Methode zum Anpassen einer Ausgleichsgeraden hängt also auch von den Annahmen über die Art der zufallsverteilten Fehler ab.
Die Abbildung rechts illustriert die Idee der Total-Least- Squares-Approximation: Zu den Datenpunkten sind die Normalabstände zur Ausgleichsgeraden (durchgezogene Linie) eingezeichnet. Die klassische kleinste-Quadrate- Ausgleichsgerade ist strichliert eingezeichnet. Sie unter- schätzt die Steigung; die TLS-Ausgleichsgerade schätzt besser.
-1 -0.5 0 0.5 1
-2 -1 0 1 2 3 4 5
Aufgabe 49: TLS-Ausgleichsebene
Wenn an Datenpunkte ∈ R3 eine Ausgleichsebene angepasst werden soll, und die Daten in x-,y- undz-Richtung durch zufällige Fehler verrauscht sind (Annahme: unabhängige normal-
verteilte Fehler, Mittelwert 0, gleiche Standardabweichung), dann liefert die TLS-Anpassung (optimiert Normalabstände zur Ausgleichsebene) eine bessere Schätzung als die klassische Kleinste-Quadrate-Anpassung (optimiert Fehler inz-Richtung).
Die Vorgangsweise ist analog zur vorigen Aufgabe. Vorgangsweise für Datenx= [x1x2. . . xn], y= [y1y2, . . . yn], z= [z1z2, . . . zn]:
• BestimmeSchwerpunkt [¯x,y,¯ z¯] der Daten.
• Verschiebe Koordinatenursprung in den Schwerpunkt
∆xi =xi−x,¯ ∆yi=yi−y,¯ ∆zi=zi−z¯
• Bilde Singulärwertzerlegung
U·S·VT =
∆x1 ∆y1 ∆z1 ... ... ...
∆xn ∆yn ∆zn
• TLS-Ausgleichsebene geht durch den Schwerpunkt [¯x,y,¯ z¯] und wird von den beiden ersten Spaltenvektoren von V aufgespannt. Sie hat als Normalvektor den letzten Spal- tenvektor vonV. Die Normalform der Ebenengleichung lautet
ax+by+cz=d . Dabei ist
a b c
der letzte Spaltenvektor von V undd =ax¯+by¯+cz¯ergibt sich durch Einsetzen des Schwerpunktes in die Ebenengleichung.
Berechnen Sie für die Daten aus Abschnitt 4.4.3
%% Datenpunkte: Variable x1, x2 und Messwert y x1 = [.2 .5 .6 .8 1.0 1.1]’;
x2 = [.1 .3 .4 .9 1.1 1.4]’;
y = [.17 .26 .28 .23 .27 .24]’;
die TLS-Ausgleichsebene und vergleichen Sie mit der klassischen Kleinsten-Quadrate-Ausgleichsebene (Musterlösung aus Abschnitt 4.4.3).
5.5 Singul¨ arwertzerlegung zur Datenkompression
Die 5. Vorlesung enthält Material zu Singulärwertzerlegung und Bildkompression. Hier folgt eine Arbeitsanleitung. Der Arbeitsauftrag lautet:
• Suchen Sie sich ein nettes Foto
• komprimieren Sie es
• Stellen Sie Original und komprimierte Version gegenüber und vergleichen Sie die Daten- größe
Die hier erklärte Idee, Approximation durch Matrizen mit niedrigem Rang, wird in der Signal- verarbeitung, Bildbearbeitung und im big data-Bereich eingesetzt. (Wenn es Sie interessiert:
Suchen Sie in Wikipedia zum Stichwort „Singulärwertzerlegung - Anwendungen“, vor allem auch auf der englischen Seite.)
Aufgabe 50: Wie man eine Katze komprimiert
Nein, wir wollen der Katze nicht wehtun. Suchen Sie sich deshalb selber ein nettes Foto aus.
Zuerst lesen Sie die Bilddatei in MATLAB ein. Wir arbeiten der Einfachheit halber in Schwarz- Weiß.
Die folgenden Code-Zeilen führen eine Standard-Umrechnung von Farbe zu Schwarz- Weiß durch.
Am besten übernehmen Sie den Code mit copy- paste; Nur die einfachen Anführungszeichen müssen Sie ersetzen; und sie sollen hier ihre ei- gene Bilddatei verwenden.
photoarray = imread(’katze2.jpg’);
redpix = single(photoarray(:,:,1))/255;
greenpix = single(photoarray(:,:,2))/255;
bluepix = single(photoarray(:,:,3))/255;
graypix = 0.299*redpix + 0.587*greenpix + 0.114*bluepix;
imagesc(graypix) colormap(gray) axis off image
Nun berechnen Sie die Singulärwertzerlegung der Pixel-Matrix
[U,S,V] = svd(graypix);
Die Diagonale derS-Matrix enthält die Singulärwerte. In diesem Zusammenhang gilt (sehr ver- einfacht): Große Singulärwerte gehören zu wichtigen Datenkomponenten. Kleine Singulärwerte gehören zu eher unwichtigen Datenmaterial. Deswegen plotten wir hier die Singulärwerte.
semilogy(diag(S))
In der halblogarithmischen Darstellung ist deutlich zu sehen: Abfall um zwei Größenord- nungen innerhalb der ersten 100 S.-Werte, da- nach langsameres Abklingen.
Je nachdem, welches Bild Sie bearbeiten, sieht Ihr Diagramm anders aus. Wenn dieser steile Abfall bei Ihren Daten auch auftritt, stehen die Chancen gut für eine Datenkompression.
Das Rezept zur SVD-Kompression: behalten Sie nur wenige Spaltenvektoren vonU undV, die zu den größten Singulärwerten gehören. Daraus können Sie das Originalbild näherungsweise rekonstruieren nach dem Rezept
Spaltenvektoren vonU malS mal transponierte Spalten von V
FürnSpalten liefert dieses Rezept eine Matrix mit Rangn. (Man spricht von einer „Rang-n- Approximation“. Die Singulärwert-Zerlegung liefert die bestmögliche Rang-n-Approximation, gemessen in der 2-Norm.) Wir probieren das für unterschiedliche n und definieren dazu eine Funktion, die auch gleich das rekonstruierte Bild anzeigt.
rekonstr = @(n) imagesc(U(:,1:n)*S(1:n,1:n)*V(:,1:n)’);
colormap(gray)
rekonstr(1); axis off image rekonstr(2); axis off image rekonstr(4); axis off image rekonstr(8); axis off image ...
rekonstr(64); axis off image
Und so wird das Kätzchen immer deutlicher sichtbar. . . Suchen Sie sich für Ihr Bild ein n aus, das zufriedenstellende Bildqualität liefert. Schmeißen Sie die nicht benötigten U- und V-Spalten weg. Wie viele Byte Speicherplatz brauchen die Daten nun?
Wie viel Speicherplatz braucht die Originaldatei? Speichern Sie dazu die Original-Schwarz- Weiß-Datei als JPEG-Datei. Die folgenden Befehle speichern und schreiben Dateiinformation aus.
imwrite(graypix, ’graysamp.jpg’, ’jpg’) d = dir(’graysamp.jpg’); d.bytes
Natürlich verwendet der JPEG-Algorithmus noch viel schlauere Tricks, aber im Prinzip sucht er auch nach niedrig-dimensionalen Approximationen von Datenmatrizen. Immerhin erhalten wir auch ein vergleichbare Kompressions-Verhältnis bei akzeptabler Bildqualität.
