Seminar zu aktuellen Themen der Numerikim Wintersemester 2010/2011

Sequential-Subspace-Optimization-Methodein der Tomographie

Andreas Platen

03. Dezember 2010

1

Inhaltsverzeichnis

Inhaltsverzeichnis

1 Einleitung 3

2 Elektronentomographie 62.1 Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Glattung der TV-Norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Sequential-Subspace-Optimization 73.1 SESOP-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83.2 Wahl der Unterraume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.2.1 Nemirovski-Richtung . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.2 Erweiterter Unterraum . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.3 SESOP in Kombination mit Truncated-Newton (SESOP-TN) . . . 113.2.4 Parallel-Coordinate-Descent und Seperable-Surrogate-Functional . 143.2.5 Ubersicht aller Unterraume . . . . . . . . . . . . . . . . . . . . . . 14

3.3 Konvergenzresultat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.4 Alternative Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

4 Numerische Experimente 224.1 Vergleich der Unterraume . . . . . . . . . . . . . . . . . . . . . . . . . . . 224.2 Anzahl der Newton-Schritte und Parameter mu . . . . . . . . . . . . . . . 234.3 Rekonstruktion aus weniger Messdaten . . . . . . . . . . . . . . . . . . . . 24

5 Ausblick 29

2

1 Einleitung

Die Computertomographie, wie sie beispielsweise in Krankenhausern zu finden ist, liefertunter Verwendung von Rontgenstrahlen Bilder vom Inneren des Objektes. Das Prinzipder Computertomographie unterscheidet sich dabei nicht von dem der Elektronentomo-graphie, bei der das Objekt mit Elektronenstrahlen beschossen wird. Die transmittier-ten Elektronen werden gemessen, woraus sich eine Intensitatsverteilung der detektiertenElektronen ergibt. Aus diesen Messdaten werden zweidimensionale Bilder von Quer-schnitten des Objektes rekonstruiert.

Es gibt zwei verschiedene Grundansatze fur die Rekonstruktion der Bilder. Ein An-satz beruht auf dem Fourier-Slice-Theorem. Die Messdaten werden als Projektionen deszweidimensionalen Querschnitts auf eindimensionale Intervalle gesehen. Ein Elektronen-strahl der das Bild unter einem bestimmten Winkel trifft, liefert eine Messung. Dieserwird senkrecht zur Messrichtung vorwarts geruckt, so dass das ganze Bild quer durch-laufen wird. Auf diese Weise erhalt man eine eindimensionale Intensitatsverteilung, vonder anschließend eine diskrete Fouriertransformation durchgefuhrt wird. Die Fourierko-effizienten entsprechen denen des Bildes in der Richtung in der vorgeruckt wird. Analogverfahrt man fur weitere Winkel, so dass man fur viele Richtungen die Fourierkoeffizi-enten des Bildes erhalt. Durch eine Rucktransformation erhalt man schließlich das Bild.

Der Vorteil dieses Verfahrens liegt im sehr geringen Rechenaufwand, wenn die schnelleFouriertransformation verwendet wird. Es ist wichtig, ein großes Spektrum von Winkelrund um das Objekt mit den Messungen abzudecken. Andernfalls fehlt ein breites Spek-trum von Fourierkoeffizienten. Werden diese vernachlassigt, kann es beispielsweise zurAusbildung Gibbsscher Phanomene kommen. In dem rekonstruierten Bild ist dement-sprechend mit Artefakten zurechnen. Es ist schwierig die fehlenden Fourierkoeffizientengeeignet zu wahlen, weshalb in der Regel Messungen von vielen Seiten notig sind. Diesist in der Elektronentomographie jedoch haufig nicht moglich. Aus diesem Grund wirdin dieser Ausarbeitung nicht naher darauf eingegangen.

Alternativ bietet sich die sogenannten”Algebraic Reconstruction Technique“ (ART)

an. Die Messdaten entsprechen dem Linienintegral des Objektes, das heißt der Dichte,entlang des Elektronen- oder Rontgenstrahls. Die Integrale des Bildes entlang den Mess-richtungen muss also gleich dem Messwert sein. Wenn das Bild als endliche Anzahl vonPixelwerten angesehen wird, vereinfacht sich das Integral zu einer gewichteten Summedieser Pixelwerte. Man erhalt also ein lineares Gleichungssystem.

Werden viele Messungen durchgefuhrt, erhalt man ein uberbestimmtes Gleichungs-system. Die Existenz einer Losung ist fur uberbestimmte Systeme nicht sichergestellt.Daher wird derjenige Vektor gesucht, der das Residuum minimiert und gleichzeitig diekleinste `2-Norm besitzt. Es wird also die Methode der kleinsten Quadrate verwendet.Das Verfahren der Wahl ist hierfur haufig der sogenannte Kaczmarz-Algorithmus, siehe[11].

3

1 Einleitung

Um die Leistungsfahigkeit tomographischer Inversionsmethoden fur bestimmte An-wendungen zu testen, werden haufig sogenannte

”Phantome“ benutzt, die sich an die

Struktureigenschaften der realen Objekte anlehnen, siehe Abbildung 1a. Dies sind alsoObjekte, die aus verschiedenen geometrischen Formen und wenig Details bestehen. Eszeigt sich, dass die Ergebnisse der Methode mit den kleinsten Quadraten nur dann gutsind, wenn genugend viele Messungen durchgefuhrt wurden. Wenn das Gleichungssys-tem unterbestimmt ist, bekommt man Ergebnisse, welche verrauscht und unscharf sind,siehe beispielsweise Abbildung 1b. Im letzten Bild wurde dabei ein Winkelbereich von±80◦ verwendet. Verkleinert man diesen Winkelbereich auf ±60◦ ohne die Anzahl derMessungen zu reduzieren, erhalt man nur noch die Abbildung 1c als Ergebnis. Wirdauch noch die Anzahl der Messungen verringert, ergibt sich Abbildung 1d.

(a) Originalbild (b) Rekonstruktion mit Kaczmarzmit 1940 Gleichungen bei 4096Unbekannten

(c) Rekonstruktion mit Kaczmarzmit 1962 Gleichungen bei 4096Unbekannten

(d) Rekonstruktion mit Kaczmarzmit 870 Gleichungen bei 4096Unbekannten

Abbildung 1: Phantom 3, 64× 64 Pixel.

4

Aktuell wird in der Forschung untersucht, ob sich die ART auch fur die Rekonstruktioneines Bildes aus wenigen Messdaten eignet. Wenige Messungen sind beispielsweise inder Computertomographie fur die Gesundheit von Mensch und Tier forderlich. UnterUmstanden ist es aber auch notig, da das Objekt sonst einer zu hohen Strahlungsdosisausgesetzt ist und zerstort werden kann. Das Problem ist dabei, ein Bild mit vielenInformationen aus so wenig Informationen wie moglich zu rekonstruieren. Dem aus denMessungen zur Verfugung stehendem Informationsgehalt entsprechend, erhalt man einunterbestimmtes Gleichungssystem. Wenn dieses Gleichungssystem losbar ist, existierenunendlich viele Losungen. Es ist nun die Losung zu wahlen, die dem realen Bild amnachsten kommt.

Einen moglicher Ansatz hierfur konnte die sich derzeit herausbildende mathematischeTheorie des Compressed Sensing liefern. Um die Idee dieses Ansatzes zu verstehen, wer-den kurz einige Eigenschaften eines digitalen Bildes beleuchtet werden. Jeder Pixel einesrealen digitalen Bildes hat einen Wert. Die Farbwerte von Pixeln in einer kleinen Umge-bung haben mit hoher Wahrscheinlichkeit einen ahnlichen Wert. Transformiert man einsolches Bild beispielsweise in Waveletdarstellung, so ergibt sich, dass die uberwiegendeZahl der Koeffizienten vernachlassigbar klein ist. Dementsprechend wird das Bild durcheinen dunnbesetzten Koeffizientenvektor charakterisiert.

Es sei nun ein solches Signal x ∈ RN nach einer geeigneten Basistransformation sogegeben, dass lediglich k � N Werte ungleich Null sind. Das Ziel des Compressed Sen-sing ist es, aus n � N Messungen des Signals einen Informationsvektor y ∈ Rn zuerhalten, so dass sich x wieder exakt rekonstruieren lasst. Der Informationsvektor wirddabei mit einem Kodierer Φ ∈ Rn×N und y = Φx konstruiert. Die Wiederherstellungdes Signals fuhrt auf ein unterbestimmtes lineares Gleichungssystem. Da Losungen die-ses Gleichungssystems nicht eindeutig sind, muss es folglich unendlich viele Dekodierergeben. Gesucht ist jedoch das Signals, welches die gewunschte dunnbesetzte Strukturhat. Dazu wird ein spezieller nichtlinearer Dekodierer ∆ verwendet, der das Gleichungs-system lost und zusatzlich noch einen Strafterm minimiert. Auf diese Weise hofft man,die Losungen zu finden, die der Eigenschaft des Signals am ehesten entspricht.

Um solche Minimierungsprobleme in einer sehr großen Dimension N zu losen, sindVerfahren notig, deren Speicher- und Rechenaufwand nicht zu stark mit der Dimensi-on N steigt. In dieser Ausarbeitung werden wir die sogenannte Sequential-Subspace-Optimization-Methode (SESOP-Methode) [7, 8] kennenlernen.

Diese Ausarbeitung ist wie folgt gegliedert: Im folgenden Abschnitt werden wir alsEinfuhrung das Beispiel der Elektronentomographie konkretisieren. Dabei wird die Pro-blemstellung auf eine Minimierung einer Funktion reduziert. Danach folgt in Abschnitt 3die Beschreibung des SESOP-Verfahrens, welches sich auf solche Minimierungsproblemeanwenden lasst. Es handelt sich dabei um ein iteratives Verfahren, welches eine vorge-gebene Funktion schrittweise in einem Unterraum minimiert. In Abschnitt 4 werden wirdieses Verfahren auf Beispiele aus der Elektronentomographie anwenden. Zum Schlussfolgt eine Zusammenfassung dieser Arbeit.

5

2 Elektronentomographie

2 Elektronentomographie

In der Elektronentomographie stoßt man mit den herkommlichen Verfahren auf Proble-me, wenn man nur wenige Messungen in einem eingeschrankten Winkelbereich durch-fuhrt. Aus diesem Grund werden wir in Abschnitt 2.1 einen weiteren ART-Ansatzkennenlernen. Es wird dabei versucht, ein Minimierungsproblem aufzustellen, welchesLosungen mit den Eigenschaften von Bild 1a liefern soll. Anschließend modifizieren wir inAbschnitt 2.2 diese Problemstellung, um auf diese Rekonstruktion das SESOP-Verfahrenanwenden zu konnen. Die SESOP-Methode folgt im Anschluss in Abschnitt 3.

2.1 Problemstellung

Bei der Rekonstruktion eines Bildes aus wenigen Messungen erhalten wir zunachst einunterbestimmtes Gleichungssystem der Form

Ax = b. (1)

Die Matrix A ∈ Rn×N hangt dabei von den Winkeln und Positionen ab, von denendie Messungen durchgefuhrt werden. Die Messergebnisse wiederum werden in b ∈ Rnfestgehalten. Die gesuchten Pixelwerte xi,j der i-ten Zeile und j-ten Spalte des Bildeswerden durch die Unbekannte x ∈ RN angegeben. Die Dimension N kann dabei sehrschnell sehr groß werden, wenn eine hohe Auflosung des Bildes erreicht werden soll. Furein Bild mit der Auflosung von 256× 256 Pixeln hat man bereits 65536 Unbekannten.

Da das Material unter Umstanden sehr empfindlich ist und man keine beliebigenWinkel verwenden kann, erhalt man moglicherweise deutlich weniger Daten, als fur einBild mit N Pixelwerten benotigt wird. Die Folge ist, dass das Gleichungssystem (1)unterbestimmt ist und es daher viele Losungen dieser Gleichung gibt. Man erwartetjedoch, dass das Bild stuckweise konstant und damit die Variation der Pixelwerte geringist. Die totale Variation (TV) eines Bildes mit einer Breite von B und Hohe von HPixeln kann beispielsweise mit der TV-Norm

‖x‖TV :=

H−1∑i=1

B−1∑j=1

(|xi,j − xi+1,j |+ |xi,j − xi,j+1|) (2)

gemessen werden. Dies stellt ein Maß fur den Unterschied zweier direkt benachbarterPixelwerte dar. Wenn das Bild eine gleichfarbige Flache zeigt, ist die TV-Norm gleichNull. Je mehr Objekte hinzukommen und je detailreicher das Bild wird, desto großer wirddiese Norm. Da wir jedoch ein Bild mit stuckweise konstanten Flachen erwarten ist unserZiel, die TV-Norm des Bildes zu minimieren und gleichzeitig das Gleichungssystem (1)zu losen. Ein Ansatz dieses Problems mit Nebenbedingung ist die Strafterm-Methode,in dem man fur ein geeignetes µ ∈ R das Problem

arg minx∈RN

1

2‖Ax− b‖22 + µ‖x‖TV (3)

6

2.2 Glattung der TV-Norm

lost. Das SESOP-Verfahren, welches in Abschnitt 3 vorgestellt wird, ist dabei ein itera-tives Verfahren, um solche Minimierungsprobleme mit sehr hoher Dimension zu losen.Dabei muss die zu minimierende Funktion mindestens einmal stetig differenzierbar sein.Der TV-Norm (2) fehlt diese Eigenschaft auf Grund der Betragsfunktionen. Daher folgtim nachsten Unterabschnitt eine Moglichkeit, mit diesem Problem umzugehen.

2.2 Glattung der TV-Norm

Eine Variante, um die Differenzierbarkeit zu gewahrleisten, ist eine Glattung der Be-tragsfunktion. Narkiss und Zibulevsky, siehe Gleichung (56) in [7], schlagen daher vor,die Betragsfunktion |s| durch eine der folgenden glatten Funktionen mit einem Glattungs-parameter ε > 0 zu ersetzen:

ψ1(s) :=√s2 + ε2, (4)

ψ2(s) :=ε|s| − log(ε|s|+ 1)

ε, (5)

ψ3(s) := ε

(∣∣∣sε

∣∣∣+1∣∣ s

ε

∣∣+ 1− 1

). (6)

Alle drei Funktionen konvergieren fur ε → 0 gegen die Betragsfunktion und sind glattgenug, um sie fur das SESOP-Verfahren zu verwenden.

3 Sequential-Subspace-Optimization

Unser allgemeines Ziel besteht darin, eine konvexe und genugend glatte Funktion

f : RN → R

ohne Nebenbedingung zu minimieren, also das Problem

arg minx∈RN

f(x) (7)

zu losen. Es kann dabei vorkommen, dass die Anzahl der Unbekannten N die Großenord-nung 104 uberschreitet, siehe zum Beispiel Abschnitt 2. Daher sind Algorithmen wich-tig, bei denen der Speicher- und Rechenaufwand pro Iteration nicht mehr als linear inN steigt. Narkiss und Zibulevsky stellten daher 2005 erstmals die Sequential-Subspace-Optimization-Methode (SESOP-Methode) [7, 8] vor.

Die SESOP-Methode basiert auf der Idee der Krylow-Teilraumverfahren, siehe zumBeispiel [6], Abschnitt 4.3 und [12], Kapitel 6–7. Dies sind iterative Verfahren zum Losenvon großen dunnbesetzten linearen Gleichungssystemen. Eine der ersten Methoden die-ser Klasse ist das Verfahren der konjugierten Gradienten (CG-Verfahren), welches 1952

7

3 Sequential-Subspace-Optimization

von Hestenes und Stiefel [5] veroffentlicht wurde. Fur eine symmetrisch positiv definiteMatrix A wird dabei anstelle des linearen Gleichungssystems Ax = b das aquivalenteProblem

arg minx∈RN

(1

2xTAx− bTx

)(8)

gelost. In jedem Iterationsschritt k wird der quadratische Ausdruck lediglich in einemspeziellen Teilraum Vk ⊂ RN minimiert, das heißt

xk+1 = arg minx∈Vk

(1

2xTAx− bTx

).

Dabei werden V0 ( V1 ( ... ( Vk jeweils durch A-orthogonale Basen aufgespannt. DieLosung xk+1 kann explizit aus der vorherigen Iteration xk berechnen werden. DiesesVerfahren terminiert, wenn schließlich Vk = RN gilt.

Dieses Prinzip macht sich die SESOP-Methode zu nutze. Anstatt das Minimum von faus (7) im gesamten Raum RN zu suchen, wird die Suche auf einen Unterraum des RNeingeschrankt. Dadurch erhalt man ein Teilproblem mit einer sehr kleinen Dimension.Daruber wird solange iteriert, bis eine Abbruchbedingung erfullt ist.

Um diesen Vorgang zu konkretisieren, folgt zunachst in Abschnitt 3.1 der Algorithmus,bevor wir in Abschnitt 3.2 auf die Wahl des Unterraums eingehen. In Abschnitt 3.3werden wir schließlich ein Konvergenzresultat vom SESOP-Algorithmus kennenlernen.Zuletzt werden in Abschnitt 3.4 alternative Losungsverfahren angesprochen.

3.1 SESOP-Algorithmus

Um gleich zu Beginn einen Uberblick uber den Algorithmus zu haben, werden wir unszunachst mit dem Aufbau des SESOP-Verfahrens beschaftigen.

Das CG-Verfahren bestimmt in jeder Iteration die nachste Naherung xk+1 des Pro-blems (8) aus der aktuellen Iterierten xk. Dies kann als Minimierung des quadratischenAusdrucks aus (8) in eine vorgegebene Richtung mit dem Ausgangspunkt xk gesehenwerden.

Dieses Problem konnen wir uns vorstellen als einen Bergwanderer, der einen Wegins Tal sucht. Dabei stellen wir uns den Funktionsgraphen als ein Gebirge vor, indemsich ein Tal befindet, welches das gesuchte Minimum reprasentiert. Zu Beginn stehtder Bergwanderer ganz oben auf einem Berg. Das Tal kann er jedoch auf Grund vonNebel nur erahnen. Um zum Tal zu gelangen wird er sich also fur eine Abstiegsrichtungentscheiden und so lange in diese Richtung gehen, bis es nicht mehr weiter abwarts geht.Nun ist er naher am Tal und kann seine Abstiegsrichtung anpassen und wieder ein Stuckweit entlang dieser neuen Richtung gehen und so weiter. Er sucht das Tal also nichtsofort im gesamten Gebiet, sondern schrittweise von der aktuellen Position aus in einebestimmte Richtung.

8

3.2 Wahl der Unterraume

Diese Art der Suche verwendet auch das SESOP-Verfahren. Anstelle von einer Ab-stiegsrichtung konnen hier jedoch mehrere Richtung verwendet. Von der aktuellen Ite-rierten xk aus wird in diesen Richtungen nach dem Minimum von f gesucht. Die Such-richtungen werden im Folgenden durch die Spalten einer Matrix D ∈ RN×r ausgedrucktund spannen den Unterraum auf, in dem das Minimum gesucht wird. Der folgende Al-gorithmus beschreibt das SESOP-Verfahren.

Algorithmus 1 SESOP-Algorithmus zum Losen von (7), siehe [7], Abschnitt 2.2, bzw.[8], Abschnitt 3.2.

1: Wahle Startvektor x0 ∈ RN2: for k = 0 to K do3: Wahle Unterraum D ∈ RN×r . siehe Abschnitt 3.24: Normiere Spalten von D5: α∗ ← arg min

α∈Rrf(xk + Dα) . z.B. mit gedampftem Newton-Verfahren

6: xk+1 ← xk + Dα∗ . Korrektur der letzten Iteration7: Prufe Abbruchbedingung8: end for9: return xk+1

Die Dimension N des ursprunglichen Problems (7) wurde nun auf das Losen von KProblemen der Form

arg minα∈Rr

f(xk + Dα)

mit typischerweise sehr viel kleinerer Dimension r � N reduziert. Aufgrund der kleinerenDimension dieser Probleme kann auf bekannte Verfahren, wie dem Newton-Verfahrenoder auch Quasi-Newton-Verfahren, siehe zum Beispiel [2], zuruckgegriffen werden. DieKonvergenz des SESOP-Algorithmus 1 hangt im Wesentlichen von der Wahl von D ab.

3.2 Wahl der Unterraume

Sei xk ∈ RN stets die k-te Iterierte vom SESOP-Algorithmus 1. Die Wahl der Spal-ten von D ist keineswegs eindeutig. Ein erster Ansatzpunkt ist das CG-Verfahren, wel-ches schnell fur quadratische Funktionen konvergiert. Aufgrund dieser Tatsache werdenwir in Abschnitt 3.2.1 und 3.2.2 versuchen, Ideen des CG-Verfahrens in die SESOP-Methode einzubinden. Im Abschnitt 3.2.3 und 3.2.4 werden wir alternative Moglichkeitensehen, die jeweils aus dem sogenannte Truncated-Newton-Verfahren bzw. der Parallel-Coordinate-Descent- und Seperable-Surrogate-Functional-Methode stammen. Um unseinen Uberblick zu verschaffen, werden im letzten Unterabschnitt 3.2.5 die hier vorge-stellten Moglichkeiten zusammengefasst.

9

3 Sequential-Subspace-Optimization

3.2.1 Nemirovski-Richtung

Das CG-Verfahrens ist eines der effizientesten Verfahren fur die Minimierung quadrati-scher Funktionen f . Aus diesem Grund bildet dieses Verfahren eine sehr wichtige Grund-lage fur die Wahl der Unterraume und den weiteren Verlauf dieser Ausarbeitung.

Zunachst klaren wir die Frage nach den Grunden fur die Effizienz des Verfahrens.Folgende drei Eigenschaften sind entscheidend, fur nahere Informationen sei auf [13]verwiesen:

a) Der aktuelle Gradient ist orthogonal zu den Richtungen aller bisherigen Schritte.

b) Der aktuelle Gradient ist orthogonal zu allen vorherigen Gradienten.

c) Der Abstieg der Funktion f in der k-ten Iteration ist mindestens von der OrdnungO(‖∇f(xk)‖22), wie sich spater noch zeigen wird.

Das CG-Verfahren kann auch fur nicht-quadratische Funktionen umformuliert werden.Die Konvergenzrate, wie sie fur quadratische Funktionen gultig ist, kann hier jedoch nichtmehr erwartet werden und unter Umstanden ist die Konvergenz nicht mehr sichergestellt,siehe [13], Kapitel 14.

Um eine geeignete Methode fur konvexe und glatte Optimierungsprobleme mit nicht-quadratischer Funktion zu konstruieren, schlug Nemirovski 1982 vor, die drei Bedingun-gen abzuschwachen. (Da die originale Quelle von Nemirovski [10] in russisch geschriebenwurde, konnte ich dies nicht nachlesen. Daher wird im Folgenden Narkiss und Zibulevsky[7], Abschnitt 1.1, zitiert.) Folgende Eigenschaften wurden von Nemirovski vorgeschla-gen:

a) Der aktuelle Gradient sollte orthogonal zu der Summe aller vorherigen Schrittesein, das heißt orthogonal zu d1

k = xk − x0.

b) Der aktuelle Gradient sollte orthogonal zu einer gewichteten Summe aller vorheri-gen Gradienten d2

k =∑k−1

i=0 wi∇f(xi) mit zuvor definierten Gewichten wi sein.

c) Die Minimierung sollte von der aktuellen Iteration aus uber einen Unterraum ge-schehen, welche den aktuellen Gradienten ∇f(xk) enthalt.

Wahlt man zum Beispiel als Spalten von D in Algorithmus 1 die drei Vektoren d1k,

d2k und ∇f(xk), dann erfullt diese Methode die eben genannten Bedingungen.

3.2.2 Erweiterter Unterraum

Mit den drei Richtungen aus dem vorherigen Abschnitt benotigt der SESOP-Algorithmus1 nur eine Auswertung des Gradienten pro Iteration. In vielen Anwendungen, wie derTomographie, siehe Abschnitt 2, benotigt diese Auswertung wenig Rechenaufwand. Da-durch kann weitere Rechenkapazitat genutzt werden, um die Suche auf einen großeren

10

3.2 Wahl der Unterraume

Unterraum des RN auszudehnen, indem man weitere Suchrichtungen vorgibt. Folglichwird der Aufwand fur das Losen des Teilproblems in Algorithmus 1 Zeile 5 ansteigen,aber dafur kann unter Umstanden die Anzahl der Iterationen reduziert werden.

Um uns fur weitere Suchrichtungen zu entscheiden, betrachten wir wieder das CG-Verfahren. Fur quadratische Funktionen wird dabei in einem Unterraum minimiert, dersamtliche letzten Schritte

pl = xl − xl−1

und die zugehorigen Gradienten ∇f(xl) fur l = 1, ..., k enthalt. Diese Bedingung wurdezuvor in Abschnitt 3.2.1 abgeschwacht, so dass dies nicht mehr erfullt sein muss. Da alsSuchrichtung auch die Richtung des Gradienten der aktuellen Iterierten zugelassen ist,wird gewahrleistet, dass der Abstieg mindestens so groß ist wie beim Gradientenverfah-ren. Um uns nun auch dem Verfahren der konjugierten Gradienten zu nahern, konnenwir als zusatzliche Suchrichtungen eine feste Anzahl letzter Schritte pk−i und letzterGradienten ∇f(xk−i) fur i = 0, 1, ..., imax verwenden. Dabei sollten nicht alle letztenSchritte und Gradienten als Suchrichtung vorgegeben werden, da sonst die Dimensiondes Problems in Algorithmus 1 Zeile 5 linear mit der Anzahl der Iterationen k wachst.Unser Ziel ist ein Verfahren zu erhalten, dessen Rechenaufwand pro Iteration linear inN steigt und nicht zusatzlich vom Iterationsindex k abhangt.

Wenn wir eine zweidimensionale Unterraum-Minimierung mit dem SESOP-Verfahrenund den Richtungen des letzten Gradienten ∇f(xk) und letzten Schrittes pk durch-fuhren, erhalten wir eine Methode, die fur quadratische Funktionen mit dem CG-Ver-fahren ubereinstimmt. Diese Eigenschaft kann vorteilhaft sein, wenn die Funktion f einegute quadratische Approximation besitzt. Wenn man mehr letzte Schritte und Gradien-ten in den Unterraum mit aufnimmt, kann dies dazu beitragen die Anzahl der Iterationenzu reduzieren. Der Anstieg des Aufwands fur eine Iteration fallt dabei verhaltnismaßiggering aus.

3.2.3 SESOP in Kombination mit Truncated-Newton (SESOP-TN)

Eine Variante des SESOP-Methode stammt von Zibulevsky [14] aus dem Jahr 2008.Diese ist eine Kombination aus der SESOP-Methode und dem sogenannten Truncated-Newton-Verfahren, siehe zum Beispiel [9]. Letzteres ist ebenfalls ein Verfahren, um dasProblem (7) zu losen.

Truncated-Newton: Dieses Verfahren besteht weitestgehend aus einem Newton-Ver-fahren. Sei die zu minimierende Funktion f konvex und zweimal stetig differenzierbar.Aus der Konvexitat folgt, dass ein lokales Minimum gleichzeitig das globale Minimum

11

3 Sequential-Subspace-Optimization

sein muss. Wegen der Differenzierbarkeit verschwindet der Gradient von f im Minimum.Das heißt die Minimierungsaufgabe ist aquivalent zu der Nullstellensuche

∇f(x) = 0.

Die Idee des Newton-Verfahrens ist es, in jedem Schritt die Nullstelle einer Linearisie-rung von ∇f zu bestimmen. Die Linearisierung von ∇f entsteht dabei aus der Taylor-entwicklung und der aktuellen Iterierten xk als Entwicklungspunkt. Es wird der Vektorx bestimmt, der die Gleichung

0 = ∇f(xk) +∇2f(xk)(x− xk) (9)

erfullt. Dazu wird das lineare Gleichungssystem ∇2f(xk)s = −∇f(xk) gelost. Die neueIterierte wird dann mit Hilfe einer Liniensuche in α ∈ R durch xk+1 = xk+αs bestimmt.Die Unbekannte α wird so bestimmt, dass die Funktion f in Richtung s minimiert wird,vergleiche dazu Zeile 5 und 6 in Algorithmus 1 mit D := s. Fur diese Wahl von D wirdder SESOP-Algorithmus 1 zum eben beschriebenen Newton-Verfahren.

Diese Nullstellensuche (9) entspricht der Minimierung der quadratischen Taylor-Ent-wicklung qk(x) von f mit dem Entwicklungspunkt xk. Die Taylor-Entwicklung ist durch

qk(x) = f(xk) +∇f(xk)T (x− xk) +

1

2(x− xk)

T∇2f(xk)(x− xk) (10)

gegeben. Der Gradient von qk entspricht gerade der rechten Seite von (9).

Der Unterschied vom Truncated-Newton-Verfahren zum gerade beschriebenen New-ton-Verfahren ist, dass diese Minimierung nur annahernd geschieht. Die lineare Glei-chung ∇2f(xk)s = −∇f(xk) wird dazu mit einer festen Anzahl von Schritten des CG-Verfahren gelost wird. Das CG-Verfahren wird also fruhzeitig abgebrochen, wodurchdas Gleichungssystem nicht exakt gelost wird. Dies erklart auch den Namen Truncated-Newton-Verfahren.

In diesem Fall endet die Newton-Iteration mit einer Liniensuche von der aktuellenIterierten xk aus in Richtung des Ergebnisses des CG-Verfahrens. Dies garantiert einenAbstieg der Funktion f . Die Effektivitat des Truncated-Newton-Verfahren hangt dabeivon der Wahl der Abbruchbedingung des CG-Verfahrens ab. Um diese Schwierigkeitzu umgehen, werden wir diese Liniensuche auf eine geeignete Unterraumminimierungausgedehnt.

Unteraumminimierung anstelle der Liniensuche: Angenommen xk sei die k-te Iteriertedes Truncated-Newton-Verfahren und das CG-Verfahren wurde nach l Schritten mit derLosung xlk abgebrochen. Zuruck vom quadratischen Model qk wollen wir den nachstenCG-Schritt imitieren und f anstelle von qk verwenden. Ausgehend von der aktuelleninneren Iterierten xlk wurde der nachste CG-Schritt l + 1 eine Minimierung von qk im

12

3.2 Wahl der Unterraume

Unterraum V lk durchfuhren. Dieser Unterraum wird dabei durch den letzten CG-Schritt

xlk − xl−1k und dem aktuellen Gradienten ∇qk(xlk) aufgespannt.Im Truncated-Newton-Verfahren wurde f an dieser Stelle lediglich in Richtung s =

xlk − xk minimiert werden. Anstelle dieser eindimensionalen Minimierung, soll mit Zeile5 aus dem SESOP-Algorithmus 1 die Funktion f vom Punkt xk aus in einem erweitertenUnterraum D := Vk ⊃ V l

k minimiert werden. Um einen monotonen Abstieg von f zuerhalten, nehmen wir die Richtung des Truncated-Newton-Verfahren

dkTN = xlk − xk

zum Unterraum Vk bzw. D hinzu.Auf diese Weise kann die Objektfunktion relativ zu f(xk) verkleinert werden. Falls

die TN-Richtungen keinen zufriedenstellenden Abstieg liefern, konnen wir optional dieRichtungen aus Abschnitt 3.2.1 und 3.2.2 in Vk aufnehmen.

An dieser Stelle sei hervorgehoben, dass das Truncated-Newton-Verfahren ein Spezi-alfall des SESOP-Algorithmus 1 ist. Die oben beschriebenen Newton-Iterationen konnendamit als SESOP-Iterationen angesehen werden.

SESOP-TN: In der k-ten Iteration des SESOP-Verfahrens wird zunachst das Newton-System

∇2f(xk)dkTN = −∇f(xk)

mit l Schritten des CG-Verfahrens annahernd gelost und damit das quadratische Modellqk(x) naherungsweise minimiert. Sei xlk die letzte Iterierte des CG-Verfahrens in derk-ten Iteration vom SESOP-Algorithmus 1. Eine weitere wunschenswerte Eigenschaftist, die Abfolge des CG-Verfahrens uber die Iterationen des SESOP-Verfahrens hinauseinzuhalten. Der SESOP-Schritt in Zeile 5 und 6 von Algorithmus 1 ist also ein Zwischen-schritt des CG-Verfahrens, welcher die Funktion f anstelle von qk in einem Unterraumminimiert. Daher kann man diesen Zwischenschritt auch als einen uberarbeiteten CG-Schritt ansehen. Fur das CG-Verfahren in der nachsten Iteration vom SESOP-Verfahren,ist es daher sinnvoll, an diesen letzten Schritt anzuknupfen. Im ersten CG-Schritt kanndazu die Minimierung des quadratischen Modells qk in dem zweidimensionalen Unter-raum durchgefuhrt, der durch xk − xlk−1 und ∇f(xk) aufgespannt wird. Dadurch wird

der Wert xl+1k−1 durch den SESOP-Schritt xk ausgetauscht.

Wir wahlen nun als Unterraum fur den SESOP-Schritt, das heißt als Spalten von D,die TN-Richtung dkTN = xlk − xk, den letzten Gradienten des quadratischen Modells qk,das heißt ∇qk(xlk), die letzte Richtung des CG-Verfahrens xlk − xl−1k und optional eineKombination der Richtungen aus Abschnitt 3.2.1 und 3.2.2.

Diese Methode lost das Problem des Truncated-Newton-Verfahren, falls das CG-Verfahren zu fruh abbricht. Fur quadratische Funktionen f stimmt der Ablauf desSESOP-TN mit dem vom CG-Verfahren, welches direkt auf die Funktion angewandt

13

3 Sequential-Subspace-Optimization

wird, uberein. Dies ist dabei unabhangig von der Abbruchbedingung des CG-Verfahrens.Dem normalen Truncated-Newton-Verfahren fehlt diese Eigenschaft, so dass hier die Ab-bruchbedingung weiterhin eine große Rolle spielt.

3.2.4 Parallel-Coordinate-Descent und Seperable-Surrogate-Functional

Der Vollstandigkeit halber werden an dieser Stelle noch weitere Varianten der SESOP-Methode erwahnt, auf die wir jedoch nicht genauer eingehen werden.

Elad, Matalon und Zibulevsky [3] schlagen noch weitere Richtungen fur den UnterraumD vor. Das Problem (7) ist dabei auf Funktionen der Form

f(x) = ‖Ax− b‖22 + ρ(x) (11)

beschrankt, wobei A ∈ Rn×N , b ∈ Rn und ρ(x) =∑N

i=1 ρi(xi) seien. Die Funktion ρkann dann beispielsweise die `1-Norm oder auch die TV-Norm (2) sein.

3.2.5 Ubersicht aller Unterraume

An dieser Stelle werden alle Richtungen zusammengefasst, die wir fur den Unterraum Daus dem SESOP-Algorithmus 1 kennengelernt haben.

Die ursprungliche Wahl von Narkiss und Zibulevsky [7], Abschnitt 2.1 bzw. [8], Ab-schnitt 3.1, ist eine Teilmenge der folgenden Richtungen als Spalten fur D, siehe Ab-schnitt 3.2.1 und 3.2.2:

• Aktueller Gradient:

∇f(xk). (12)

• Nemirovski-Richtungen:

d(1)k := xk − x0,

d(2)k :=

k∑i=0

wi∇f(xi),(13)

wobei

wi :=

{1 fur i = 0,

12 +

√14 + w2

i−1 fur i > 0.(14)

• Vorherige Richtungen:

pk−1 = xk−i − xk−i−1, i = 0, 1, ..., s1. (15)

14

3.3 Konvergenzresultat

• Vorherige Gradienten:

∇f(xk−i), i = 1, 2, ..., s2. (16)

Zusatzlich zu diesen Richtung konnen die Richtungen von SESOP-TN [14] gemaß Ab-schnitt 3.2.3 gewahlt werden. Dazu muss zunachst die Gleichung

∇2f(xk)dkTN = −∇f(xk)

mit l Schritten des CG-Verfahrens gelost werden. Dabei wird im ersten Schritt mit demzweidimensionalen Unterraum begonnen, der durch xk − xlk−1 und ∇f(xk) aufgespanntwird. Als Unterraum ergeben sich zusatzlich zu den bereits aufgelisteten Richtungen:

• Truncated-Newton-Richtung:

dkTN = xlk − xk.

• Letzter Gradient des quadratischen Modells qk aus (10):

∇qk(xlk).

• Letzte Richtung des CG-Verfahrens:

xlk − xl−1k .

Wenn f die Form (11) hat, kann zusatzlich eine der Richtungen von Elad, Matalonund Zibulevsky [3] verwendet werden:

• Richtung aus dem Parallel-Coordinate-Descent- oder Seperable-Surro-gate-Functional-Verfahren.

3.3 Konvergenzresultat

In diesem Abschnitt werden wir sehen, dass fur eine genugend glatte Funktion f derSESOP-Algorithmus 1 gegen die exakte Losung konvergiert. Als Suchrichtung mussenwir dazu den aktuellen Gradienten (12) und die Nemirovski-Richtungen (13) verwenden.Der folgende Satz liefert bereits diese Aussage.

Satz 1. (siehe [7] Theorem 1)Sei f ∈ C2(RN ) konvex und ∇f Lipschitz-stetig mit Lipschitz-Konstante L, das heißt

‖∇f(x)−∇f(y)‖2 ≤ L‖x− y‖2 fur alle x,y ∈ RN . (17)

15

3 Sequential-Subspace-Optimization

Weiter sei der Unterraum {Dα : α ∈ Rr} im SESOP-Algorithmus 1 durch den aktuel-len Gradienten (12), den Nemirovski-Richtungen (13) und optional durch verschiedeneweitere Richtungen aufgespannt. Die exakte Losung x∗ von Problem (7) sei maximal mitdem Abstand R vom Startvektor x0 entfernt, das heißt

‖x∗ − x0‖2 ≤ R. (18)

Dann gilt fur den SESOP-Algorithmus 1 angewendet auf die Funktion f die Ungleichung

εK+1 := f(xK+1)− f(x∗) ≤ LR2

K2. (19)

Satz 1 wird im Folgenden in vier Schritten bewiesen. Auf Grund einer besseren Uber-sicht, werden die ersten drei Schritte in folgenden drei Hilfssatzen formuliert.

Lemma 1. (siehe [7] Proposition 1)Unter den Voraussetzungen von Satz 1 gilt, dass

f(xk+1) ≤ f(xk)−‖∇f(xk)‖22

2L.

Beweis. Zuerst zeigen wir, dass

vT∇2f(x0)v ≤ L‖v‖22 fur alle x0,v ∈ RN (20)

gilt. Dazu betrachten wir die Taylor-Entwicklung von ∇f(x) mit Entwicklungspunkt x0,welche fur ein t 6= 0 und v ∈ RN durch

∇f(x0 + t · v) = ∇f(x0) +∇2f(x0)T tv +O(‖tv‖22)

gegeben ist. Durch Umformung der Gleichung erhalten wir

∇f(x0 + t · v)−∇f(x0)

t=∇2f(x0)

T tv +O(‖tv‖22)t

. (21)

Als Grenzwert der rechten Seite von (21) fur t→ 0 ergibt sich

limt→0

∇2f(x0)T tv +O(‖tv‖22)

t= ∇2f(x0)v, (22)

denn es gilt

limt→0

∣∣∣∣O(‖tv‖22)t

∣∣∣∣ = ‖v‖2 limt→0

O(‖tv‖22)|t|‖v‖2

= ‖v‖2 limt→0

O(‖tv‖22)‖tv‖2

= ‖v‖2 · 0 = 0.

16

3.3 Konvergenzresultat

Fur die linke Seite von (21) gilt wegen (17), dass∥∥∥∥∇f(x0 + t · v)−∇f(x0)

t

∥∥∥∥2

≤ L‖v‖2.

Da die obere Schranke unabhangig von t ist, folgt zusammen mit (21) und (22) fur denGrenzwert t→ 0, dass

‖∇2f(x0)v‖2 ≤ L‖v‖2 fur alle x0 ∈ RN .

Die Cauchy-Schwarz-Ungleichung liefert

|vT∇2f(x0)v| ≤ ‖v‖2‖∇2f(x0)v‖2 ≤ L‖v‖22

und damit Ungleichung (20). Die ersten beiden Richtungsableitungen von f(x) in Rich-tung des Gradienten ∇f(x) =: g(x) sind durch

f ′g(x) = ∇f(x)T∇f(x) = ‖∇f(x)‖22 (23)

und

f ′′gg(x) = ∇f(x)T∇2f(x)∇f(x) ≤ L‖∇f(x)‖22 (24)

gegeben, wobei die letzte Ungleichung aus (20) folgt.

Wir betrachten nun die Minimierung von f(x) entlang der Abstiegsrichtung −∇f(x)an einem beliebigen Punkt x ∈ RN . Sei ϕ(α) := f(x−α∇f(x)) eine Funktion und q(α)ein quadratisches Polynom mit den Eigenschaften

q(0) = ϕ(0),

q′(0) = ϕ′(0)

und

q′′(α) = L‖∇f(x)‖22

fur beliebiges α, da die zweite Ableitung konstant ist. Wegen (20) gilt die Ungleichung

q′′(α) = L‖∇f(x)‖22 ≥ ∇f(x)T∇2f(x− α∇f(x))∇f(x) = ϕ′′(α)

und damit

q(α) ≥ ϕ(α)

17

3 Sequential-Subspace-Optimization

fur jedes α ∈ R. Dadurch ist sichergestellt, dass die Funktion f bei exakter Liniensucheum einen Wert ∆f(x) fallt. Dieser Wert ist großer als der Abstieg von q bis zu seinemMinimum, das heißt

∆f(x) ≥ q′(0)2

2q′′. (25)

Dies wird in Abbildung 2 verdeutlicht. Mit (23) und (24) erhalten wir

∆f(x) ≥ ‖∇f(x)‖222L

.

Da der Gradient in den Richtungen des SESOP-Verfahrens enthalten ist, folgt die Be-hauptung von Lemma 1.

q(α)

∆f

q′(0)2

2q′′

ϕ(α)

Abbildung 2: Funktion ϕ und seine Majorante q liefern (25), siehe [7] Abbildung 1.

Lemma 2. (siehe [7] Proposition 2)Unter den Voraussetzungen von Satz 1 gilt, dass

εk ≤ (∇f(xk))T (x0 − x∗). (26)

Beweis. Sei εk := f(xk)− f(x∗), wie in (19). Aus Lemma 1 erhalten wir

εk − εk+1 ≥‖∇f(xk)‖22

2L. (27)

18

3.3 Konvergenzresultat

Wir betrachten nun die Funktion f entlang der Richtung xk − x∗. Sei dazu

ϕ(α) := f(xk − α(xk − x∗)).

Aus der Konvexitat von f folgt, dass auch ϕ konvex und damit speziell die Ungleichung

ϕ(0)− ϕ(1) ≤ ϕ′(0) · 1

erfullt ist. Durch die Substitution f(xk) = ϕ(0) und f(x∗) = ϕ(1) erhalten wir

εk = f(xk)− f(x∗) ≤ ∇f(xk)T (xk − x∗). (28)

Bei der Bestimmung von xk wird die Funktion f in die Nemirovski-Richtung xk−1 − x0

und die tatsachlich verwendete Richtung xk − xk−1 minimiert. Daher folgt, dass

∇f(xk) ⊥ xk − x0 (29)

ist. Mit (28) und (29) folgt die Behauptung von Lemma 2, denn es gilt

εk ≤ ∇f(xk)T (xk − x0 + x0 − x∗)

= ∇f(xk)T (xk − x0) +∇f(xk)

T (x0 − x∗)

= ∇f(xk)T (x0 − x∗).

Lemma 3. (siehe [7] Proposition 3)Unter den Voraussetzungen von Satz 1 mit gegebenen Gewichten wk gilt, dass

K∑k=0

wkεk ≤ R

√√√√2LK∑k=0

w2k(εk − εk+1). (30)

Beweis. Die Summe von (26) uber alle Schritte mit Gewichten wk ergibt

K∑k=0

wkεn ≤

(K∑k=0

wk∇f(xk)

)T(x0 − x∗).

Mit der Cauchy-Schwarz-Ungleichung und (18) erhalten wir

K∑k=0

wkεn ≤

∥∥∥∥∥k∑k=0

wk∇f(xk)

∥∥∥∥∥2

·R. (31)

19

3 Sequential-Subspace-Optimization

Da die Nemirovski-Richtung∑K

k=0wk∇f(xk) eine Richtung des Unterraums des Algo-rithmus ist, gilt

∇f(xk) ⊥k−1∑l=0

wl∇f(xl).

Der Satz von Pythagoras, siehe Abbildung 3, zusammen mit (27) liefert∥∥∥∥∥K∑k=0

wk∇f(xk)

∥∥∥∥∥2

2

=K∑k=0

w2k‖∇f(xk)‖22 ≤ 2L

K∑k=0

w2k(εk − εk+1). (32)

Einsetzen von (32) in (31) liefert die Behauptung von Lemma 3.

x1

x2

x3

y1

y2

Abbildung 3: Illustration zu Gleichung (32): Wenn y2 = x1 +x2 +x3 und jeder der dreiVektoren xi rechtwinklig zur Summe der vorherigen Vektoren ist, dannfolgt ‖y2‖22 = ‖y1‖22+‖x3‖22 = (‖x1‖22+‖x2‖22)+‖x3‖22, siehe [7] Abbildung2.

Es folgt der letzte Schritt fur den Beweis von Satz 1.

Beweis von Satz 1. Aus der Monotonie der εk von Lemma 3 folgt die Ungleichung

εK ≤ R√

2L

√∑Kk=0w

2k(εk − εk+1)∑K−1

k=0 wk.

Mit dem Ziel, Satz 1 zu beweisen, zeigen wir, dass man mit optimalen Gewichten einbesseres Ergebnis erzielen kann. Durch eine Umordnung von (30) erhalten wir

K∑k=0

wkεk ≤ R√

2L√w20ε0 + (w2

1 − w20)ε1 + ...+ (w2

K − w2K−1)εK − w2

KεK+1.

20

3.4 Alternative Verfahren

Sei s :=∑K

k=0wkεk. Zur Abschatzung von εK , wird gezeigt, dass die Ungleichung

s ≤ R√

2L√s− w2

KεK+1 (33)

erfullt ist. Fur diesen Zweck konnen wir die Gewichte so wahlen, dass

wk =

{w20 fur k = 0,

w2k − w2

k−1 fur k > 0(34)

gilt. Eine Losung sind dabei die Gewichte aus (14), dass heißt

wk =

{1 fur k = 0,

12 +

√14 + w2

k−1 fur k > 0,

wobei der letzte Term die großere Wurzel von (34) ist. Fur große k ist wk ≈ k2 . Wegen

(33) gilt

w2KεK+1 ≤ s−

s2

2LR2.

Unser Ziel ist, eine obere Schranke fur εK zu finden. Wir sind also an dem”schlechtesten“

Wert von s, das heißt an

s = arg maxs

{s− s2

2LR2

},

interessiert. Das Maximum wird an der Stelle s = LR2 angenommen, so dass wir schließ-lich

εK+1 ≤LR2

4w2K

≈ LR2

N2

erhalten. Dies vervollstandigt den Beweis von Satz 1.

3.4 Alternative Verfahren

Elad und Zibulevsky geben in [4] eine Ubersicht uber Methoden zur Losung von Proble-men der Form (11). Unter anderem wird dort auch die SESOP-Methode aufgegriffen.

Ein weiteres Verfahren ist NESTA von Becker, Bobin und Candes [1]. Die AbkurzungNESTA steht fur

”Nesterov’s algorithm“. Es handelt sich um ein iteratives Verfahren

welches sich beispielsweise auf Probleme der Form

arg minx∈RN

1

2‖Ax− b‖22 + µ‖x‖1,

21

4 Numerische Experimente

anwenden lasst, wobei µ ein vorgegebener Parameter ist. Um die Differenzierbarkeit derFunktion zu gewahrleisten, glattet NESTA die Funktion und versucht das Problem

arg minx∈RN

1

2‖Ax− b‖22 + µ‖x‖1+λ

zu losen. Der Glattungsparameter λ > 0 wird dabei im Laufe der Iterationen so weit wiemoglich verringert, so dass die originale Funktion angenahert wird.

4 Numerische Experimente

An dieser Stelle werden wir einige Resultate fur das Problem (3) aus der Elektronenmi-kroskopie sehen. Zur Glattung der TV-Norm wurde die Funktion ψ3 aus (6) mit ε = 10−3

verwendet. Der exakte Querschnitt, den wir aus Messungen rekonstruieren wollen, ist inAbbildung 1a zu sehen.

Dieses Bild wird nun virtuell in einen Tomographen gelegt. Dabei wird davon aus-gegangen, dass keine Messfehler gemacht werden. Wenn nichts anderes angegeben ist,werden die Messungen aus 30 verschiedenen Winkeln zwischen ±80◦ vorgenommen. DerAbstand zwischen den Messungen aus gleichen Winkeln ist 0.02 fur ein Bild aus demEinheitsquadrat [0, 1]2. Damit ergeben sich 1940 Gleichungen, das heißt A ∈ R1940×4096

und b ∈ R1940. Als Startvektor fur den SESOP-Algorithmus wurde stets der Nullvektor,das heißt ein schwarzes Bild, verwendet. In samtlichen Beispielen wurden 1000 SESOP-Iterationen berechnet. Die Iterationen sind in keinem dieser Durchlaufe fruhzeitig abge-brochen.

Zunachst werden wir in Abschnitt 4.1 den Unterraum D aus dem SESOP-Algorithmus1 mit verschiedenen Suchrichtungen untersuchen. Fur die Zeile 5 in Algorithmus 1 wirdein gedampftes Newton-Verfahren fur Systeme verwendet. Zudem muss der Parameter µin Problem (3) geeignet gewahlt werden. In Abschnitt 4.2 werden wir daher Anzahl derNewton-Schritte und den Parameter µ variieren. Als letzte numerische Experimente wer-den wir in Abschnitt 4.3 versuchen, das Bild aus weniger Messdaten zu rekonstruieren,siehe Abschnitt 2.1.

Notation: Mit SESOP-s1 wird das SESOP-Verfahren mit s1 letzten Richtungen (15)bezeichnet. Falls die Nemirovski-Richtungen (13) verwendet werden, schreiben wir imFolgenden SESOP-s+1 , ansonsten SESOP-s−1 . Das Verfahren SESOP-TN aus Abschnitt3.2.3 mit l CG-Schritten wird mit SESOP-s1-TN-l bezeichnet.

4.1 Vergleich der Unterraume

Im Folgenden werden wir verschiedene Kombinationen der Richtungen aus Abschnitt3.2 fur das Problem (3) mit µ = 0.01 testen. Der aktuelle Gradient (12) wird stets in die

22

4.2 Anzahl der Newton-Schritte und Parameter mu

Unterraumminimierung mit aufgenommen. In Tabelle 1 wurden zusatzlich unterschied-lich viele letzte Richtungen (15) und die Nemirovski-Richtungen (13) als Suchrichtungenverwendet. Fur das Newton-Verfahren wurde lediglich eine Iteration verwendet. Es zeigtsich, dass die Nemirovski-Richtungen die Minimierung der Funktion stark verlangsa-men. Dies verbessert sich, wenn die Newton-Iterationen erhoht werden. Dennoch sinddie Ergebnisse bei diesem Beispiel ein wenig schlechter als ohne diese Richtungen.

Die Hinzunahme der zuletzt verwendeten Richtungen (15) beschleunigt den Abstiegund sorgt auch fur subjektiv bessere Ergebnisse. Das Bild von SESOP-8− und SESOP-64− sind kaum von Abbildung 4c zu unterscheiden. Das heißt mehr als vier Richtun-gen machen fur dieses Beispiel scheinbar kaum noch einen Unterschied, so dass in denweiteren Experimenten uberwiegend SESOP-4− verwendet wird. Diese Losung durchSESOP-4− ist dabei praktisch nicht mehr von dem originalen Bild 1a zu Unterscheiden.Die Losung wurde dabei in nur 7 Sekunden bestimmt. Fur die Losung des Kaczmarz-Algorithmus aus Abbildung 1b wurden hingegen 105 Iterationen berechnet, welche deut-lich langer als 7 Sekunden auf der gleichen Maschine benotigen.

In Tabelle 2 wurde SESOP-TN verwendet. Diese neuen Richtungen sorgen fur SESOP-0− fur einen schnelleren Abstieg und verbessern sichtbar das Ergebnis. Die Funktion fdes Problems (3) ist nach 1000 SESOP-Schritten deutlich kleiner, wenn man mehr CG-Schritte durchfuhrt. Bei SESOP-4− hingegen verschlechtert sich das Ergebnis ein wenig,je mehr CG-Schritte man verwendet.

4.2 Anzahl der Newton-Schritte und Parameter µ

Bisher wurde nur eine Newton-Iteration fur die Unterraum-Minimierung in Zeile 5 vonAlgorithmus 1 und der Parameter µ = 0.01 verwendet. Mit SESOP-4− wurden im letz-ten Abschnitt bereits sehr gute Ergebnisse erzielt. Daher dient dieses als

”Referenz“-

Verfahren fur diesen Abschnitt.

Der Tabelle 3 ist zu entnehmen, dass die Anzahl der Newton-Iterationen fur diesesBeispiel keinen Einfluss auf das Ergebnis hat. Verkleinert man hingegen den Wert µ, wirdauch das Residuum kleiner. Die TV-Norm ist immer von der gleichen Großenordnung,wobei hier eine ansteigende Tendenz zu erkennen ist, wenn µ verkleinert wird. DieseBeobachtung ist dadurch zu erklaren, dass die TV-Norm bei einem großen µ einen großenEinfluss auf die zu minimierende Funktion in (3) hat. Man erlaubt dem Residuum alsoeinen großeren Spielraum fur Fehler. Bei einem kleinen µ wird die TV-Norm unwichtigersorgt dafur, dass das Residuum entsprechend kleiner bleiben muss.

Zu beachten ist, dass wir von perfekten Messungen ohne Fehler ausgehen. Sind dieMessdaten hingegen verrauscht, so ist die exakte Losung keine Losung der GleichungAx = b. Dem Residuum in (3) muss daher ein großerer Spielraum gewahrleistet werden,da die Gleichung nicht exakt gelost werden soll. Dies kann durch erhohen des Parametersµ erreicht werden.

23

4 Numerische Experimente

Das Bild 6d mit SESOP-4− mit einem Newton-Schritten und dem Parameter µ =0.001 erzielt dabei das beste Ergebnis und ist vom originalen Bild 1a nicht mehr zuunterscheiden.

4.3 Rekonstruktion aus weniger Messdaten

Im vorherigen Unterabschnitt haben wir ein sehr gutes Ergebnis erzielt. Wir hatten dabei1940 Gleichungen fur 4096 Unbekannte zur Verfugung. Das Ziel ist es, die Anzahl derGleichungen so weit wie moglich zu reduzieren.

In Tabelle 4 wurde abwechselnd die Spanne der Winkel von ±80◦ auf ±60◦ und dieAnzahl der Winkel von 30 auf 20 verringert und der Abstand der Messungen von h = 0.02auf h = 0.03 erhoht. Die Werte der Tabelle 4 sehen dabei sehr ahnlich zu denen ausTabelle 3 aus, wobei dort das Originalbild 1a nahezu exakt rekonstruiert wurde. Wennman nur die Anzahl der Winkel reduziert oder den Messabstand erhoht hat, dann sinddie subjektiven Ergebnisse in Form des Bildes ebenfalls sehr gut, siehe Abbildung 7dund 7c. Hingegen ist die Spanne, in der die Winkel gewahlt werden, ein sehr wichtigerParameter. Wenn man diesen auf ±60◦ reduziert, erhalt man die Abbildungen 7b und7a. Es lassen sich deutliche Probleme bei dem

”h“-formigen weißen Gebilde unten links

erkennen. Die schmale Kante kann auf diese Weise nicht mehr rekonstruiert werden.Wenn die Objekte großer und breiter sind, werden diese auch hier noch gut erkannt.

In der Einleitung wurden Ergebnisse mit dem Kaczmarz-Algorithmus und einem Win-kelbereich von±60◦ gezeigt. Mit den gleichen Messungungen verbessert sich das Ergebnisvon Kaczmarz aus Abbildung 1c mit Hilfe von SESOP zu Abbildung 7c. Verwendet mannur 870 Gleichungen ergab sich als Ergebnis zunachst 1d. SESOP lieferte eine deutlicheVerbesserung, siehe Abbildung 7d.

24

4.3 Rekonstruktion aus weniger Messdaten

Zeit Residuum TV-Norms±1 (in s) ‖AxK − b‖2 ‖xK‖TV f(xK) Bild

0− 3.38 28.7756 57373.4 602.51 4a0+ 5.31 354.968 136792 1722.88 4b4− 7.07 1.41091 44396.6 445.377 4c4+ 9.24 63.7326 105287 1116.6 4d8− 10.64 2.54125 44497.4 447.5158+ 13.01 44.9916 87677.5 921.76764− 111.61 1.05297 44423.2 445.28564+ 116.53 39.6631 82187.2 861.535

Tabelle 1: Vergleich zwischen SESOP-s±1 mit K = 1000 SESOP-Iterationen und jeweilseinem Newton-Schritt fur Problem (3) mit µ = 0.01.

(a) SESOP-0− (b) SESOP-0+

(c) SESOP-4− (d) SESOP-4+

Abbildung 4: Vergleich zwischen SESOP-s±1 mit und ohne Nemirovski-Richtungen, sieheTabelle 1.

25

4 Numerische Experimente

Zeit Residuum TV-Norms1 l (in s) ‖AxK − b‖2 ‖xK‖TV f(xK) Bild

0− 0 3.4 28.7756 57373.4 602.51 4a0− 1 10.63 30.7264 58956.2 620.2880− 8 30.66 25.3926 55318.4 578.577 5a0− 32 97.33 8.69904 46886 477.559 5b

4− 0 7.13 1.41091 44396.6 445.377 4c4− 1 14.7 2.91001 44808.2 450.9924− 8 34.84 4.22879 44940.5 453.634 5c4− 32 102.75 4.83312 45296 457.793 5d

Tabelle 2: Vergleich zwischen SESOP-s−1 -TN-l mit K = 1000 SESOP-Iterationen undjeweils einem Newton-Schritt fur Problem (3) mit µ = 0.01.

(a) SESOP-0−-TN-8 (b) SESOP-0−-TN-32

(c) SESOP-4−-TN-8 (d) SESOP-4−-TN-32

Abbildung 5: Vergleich zwischen SESOP-s−1 -TN-l, siehe Tabelle 2.

26

4.3 Rekonstruktion aus weniger Messdaten

Newton- Zeit Residuum TV-Norms1 Schritte (in s) ‖AxK − b‖2 ‖xK‖TV µ f(xK) Bild

4− 1 6.89 154.758 42408.8 0.1 4395.64 6a4− 4 20.97 172.88 42389.5 0.1 4411.834− 16 33.96 170.829 42384.7 0.1 4409.3 6b

4− 1 6.98 1.41091 44396.6 0.01 445.377 4c4− 4 21.01 3.13235 44602.7 0.01 449.164− 16 34.86 3.11136 44578.3 0.01 448.895

4− 1 6.9 0.0657363 44615 0.001 44.6808 6c4− 4 20.95 0.0740766 44632.6 0.001 44.70674− 16 31.49 0.0777603 44640.2 0.001 44.7179 6d

Tabelle 3: SESOP-4− mit K = 1000 Schritten fur verschiedene Werte µ und unterschied-lich viele Newton-Schritte in Zeile 5 von Algorithmus 1.

(a) Ein Newton-Schritt, µ = 0.1 (b) 16 Newton-Schritte, µ = 0.1

(c) Ein Newton-Schritt, µ = 0.001 (d) 16 Newton-Schritte, µ = 0.001

Abbildung 6: SESOP-4− mit verschiedenen µ-Werten und unterschiedlicher Anzahl anNewton-Schritten, siehe Tabelle 3.

27

4 Numerische Experimente

Winkel Mess- Zeit Residuum TV-NormBereich # abstand (in s) ‖AxK − b‖2 ‖xK‖TV f(xK) Bild

±80◦ 30 0.02 6.9 0.0657363 44615 44.6808 6c±80◦ 20 0.02 6.46 0.0724943 44604.4 44.6769 7a±80◦ 30 0.03 6.44 0.0775321 44669.3 44.7468 7b±60◦ 30 0.02 6.98 0.359255 43455.5 43.8148 7c±60◦ 20 0.03 6.03 0.185723 42221.1 42.4068 7d

Tabelle 4: SESOP-4− mit K = 1000 SESOP-Iterationen, µ = 0.001 und 1 Newton-Schritten mit weniger Gleichungen. In der Tabelle von oben nach unten istdie Anzahl der Gleichungen: 1940, 1294, 1296, 1962, 870.

(a) s = 80, r = 20, h = 0.02 (b) s = 80, r = 30, h = 0.03

(c) s = 60, r = 30, h = 0.02 (d) s = 60, r = 20, h = 0.03

Abbildung 7: SESOP-4− fur das Problem mit r Winkeln in einer Spanne von s und mitMessabstand h, siehe Tabelle 4.

28

5 Ausblick

Numerische Verfahren fur Minimierungsprobleme mit sehr hoher Dimension gewinnenimmer mehr an Bedeutung. In dieser Ausarbeitung haben wir dazu das Sequential-Subspace-Optimization-Verfahren von Narkiss und Zibulevsky [7, 8] kennengelernt. Da-bei wird die zu minimierende Funktion nicht direkt im gesamten RN , sondern schrittweisein einem kleinen Unterraum minimiert.

Untersucht haben wir diese Methode an einem Beispiel aus der Elektronentomogra-phie. Es hat sich dabei gezeigt, dass die Wahl der Unterraume eine große Rolle spielt.Selbst wenn mit einem Unterraum gute Ergebnisse erzielt werden, erhalt man damitnoch keine Garantie, dass dies auch bei Vergroßerung dieses Unterraums der Fall ist.

Speziell in der Tomographie ist die Spanne der Winkel, aus denen die Messungendurchgefuhrt werden sehr wichtig. Auf Grund der numerischen Experimente aus Ab-schnitt 4.3 sollten an dieser Stelle nach Verbesserungen gesucht werden sollte.

Zunachst konnen die Richtungen von Elad, Matalon und Zibulevsky [3] aus der Paral-lel-Coordinate-Descent- und Seperable-Surrogate-Functional-Methode untersucht wer-den. Diese lassen sich auf das Problem (3) aus der Tomographie anwenden, vergleicheGleichung (11).

Eine weitere Verbesserung konnte ein Quasi-Newton-Verfahren, siehe beispielsweise[2], anstelle des gedampften Newton-Verfahrens in Zeile 5 von Algorithmus 1 sein. Dieskonnte die Probleme beheben, die auftreten konnen, wenn man den Unterraum durchweitere Suchrichtungen erganzt. Die Hinzunahme von Suchrichtungen sollte das Ergebnishochstens verbessern, aber mindestens so gut sein, wie ohne diese neuen Richtungen. Dieshat sich in den Experimenten jedoch nur teilweise bestatigt.

Da die Wahl der Unterraume von großer Bedeutung ist, konnte man zudem versuchenin jedem SESOP-Schritt neu zu entscheiden, welche Suchrichtungen verwendet werdensollen. Wenn sich beispielsweise herausstellt, dass die zu minimierende Funktion f indem aktuellen Bereich gut mit der quadratischen Approximation ubereinstimmt, kannes sinnvoll sein das SESOP-TN-Verfahren zu verwenden, da dieses das CG-Verfahrengut imitieren kann. Andernfalls greift man auf andere Richtungen zuruck.

Literatur

[1] Becker, Stephen, Jerome Bobin und Emmanuel Candes: NESTA: A Fast andAccurate First-order Method for Sparse Recovery. Technischer Bericht, CaliforniaInstitute of Technology, Pasadena, 2009.

[2] Dennis, John E. und Robert B. Schnabel: Numerical Methods for Unconstrai-ned Optimization and Nonlinear Equations. Prentice-Hall series in computationalmathematics. Prentice-Hall, Englewood Cliffs, 1983.

29

Literatur

[3] Elad, Michael, Boaz Matalon und Michael Zibulevsky: Coordinate andSubspace Optimization Methods for Linear Least Squares with Non-Quadratic Re-gularization. Applied and Computational Harmonic Analysis, 23(3):346–367, 2007.

[4] Elad, Michael und Michael Zibulevsky: L1-L2 Optimization in Signal andImage Processing: Iterative Shrinkage and Beyond. IEEE Signal Processing Maga-zine, 27(3):78–88, 2010.

[5] Hestenes, Magnus R. und Eduard Stiefel: Methods of Conjugate Gradients forSolving Linear Systems. Journal of Research of the National Bureau of Standards,49(6):409–436, 1952.

[6] Meister, Andreas: Numerik linearer Gleichungssysteme. Vieweg, Wiesbaden,2005.

[7] Narkiss, Guy und Michael Zibulevsky: Sequential Subspace Optimization Me-thod for Large-Scale Unconstrained Problems. Technischer Bericht CCIT Num-mer 559, Electrical Engineering Department, Technion, Haifa, 2005. http://ie.

technion.ac.il/~mcib/sesop_report_version301005.pdf.

[8] Narkiss, Guy und Michael Zibulevsky: Support Vector Machine via SequentialSubspace Optimization. Technischer Bericht CCIT Nummer 557, Electrical Engi-neering Department, Technion, Haifa, 2005. http://ie.technion.ac.il/~mcib/

sesop_svm_report.pdf.

[9] Nash, Stephen G.: A Survey of Truncated-Newton Methods. Journal of Compu-tational and Applied Mathematics, 124(1–2):45–59, 2000.

[10] Nemirovski, Arkadi: Orth-method for smooth convex optimization. Izvestiia Aka-demii Nauk SSSR, Ser. Tekhnicheskaya Kibernetika, 2, 1982. Russisch.

[11] Popa, Constantin und Rafal Zdunek: Kaczmarz extended algorithm for tomo-graphic image reconstruction from limited-data. Mathematics and Computers inSimulation, 65(6):579–598, 2004.

[12] Saad, Yousef: Iterative methods for sparse linear systems. Society for Industrialand Applied Mathematics, Philadelphia, 2. Auflage, 2007.

[13] Shewchuk, Jonathan R.: An Introduction to the Conjugate Gradient Me-thod Without the Agonizing Pain. School of Computer Science Carnegie Mel-lon University, Pittsburgh, 1994. http://www.cs.cmu.edu/~quake-papers/

painless-conjugate-gradient.pdf.

30

Literatur

[14] Zibulevsky, Michael: SESOP-TN: Combining Sequential Subspace Optimizationwith Truncated Newton method. Technischer Bericht, Electrical Engineering De-partment, Technion, Haifa, 2008. http://ie.technion.ac.il/~mcib/sesoptn_

paper1.pdf.

31