Quadratmittelapproximation durch B-Splines mit freien Knoten

Diplomarbeit

Martin-Luther-Universität Halle/Wittenberg Fachbereich Mathematik und Informatik

eingereicht von Torsten Schütze geb. am 24.09.1969 in Nordhausen

Betreuer: Prof. Dr. habil. H. Schwetlick

Technische Universität Dresden

Institut für Wissenschaftliches Rechnen

Halle, im Juli 1993

Schütze, Torsten

Quadratmittelapproximation durch B-Splines mit freien Knoten, 1993, 66 S., 40 Lit., 16 Abb. Diplomarbeit, Martin-Luther-Universität Halle/S., Fachbereich Mathematik und Informatik

Die fehlerbehafteten Meßwerte einer glatten univariaten Funktion sollen durch einen Spline approximiert werden. Der Spline wird als Linearkombination von B-Splines dargestellt. Eine Teilmenge der Splineknoten, die freien Knoten, wird in den Optimierungsprozeß einbezogen. Minimiert wird die Summe der Fehlerquadrate und ein Glättungsterm. Es wird ein nichtlineares Quadratmittelproblem in den Splinekoeffizienten und den freien Knoten entwickelt. Das resultierende separable, restringierte Quadratmittelproblem wird auf ein Problem, in dem nur noch die freien Knoten auftreten, reduziert. Das reduzierte Problem wird mit einem Verallgemeinerten Gauß-Newton-Verfahren gelöst.

An dieser Stelle möchte ich Herrn Prof. Dr. habil. H. Schwetlick für die kontinuierliche Betreuung sowie für die wertvollen Hinweise und Anregungen meinen herzlichen Dank aussprechen.

k  Ordnung der B-Splines
n  Anzahl der B-Splines
m  Anzahl der Meßpunkte
l  Anzahl der freien Splineknoten
r  Ordnung der Ableitung im Glättungsterm
(x, y)  Meßdaten (x_{i}, y_{i}), i = 1, …, m
p  p = (p(1), …, p(l)) Zeiger auf die freien Knoten
t  (finite) Knotenfolge t = (t_{1}, …, t_{n}, …, t_{n+k})^{T}
t̃  freie Splineknoten t̃ = (t_{p(1)}, …, t_{p(l)})^{T}

B_{j,k}(x, t), B_{j,k}(x) j-ter B-Spline der Ordnung k bezüglich der Knotenfolge t

an der Stelle x

α  Splinekoeffizienten α = (α_{1}, …, α_{n})^{T}
S_{k,t}  Splineraum der Ordnung k bezüglich der Knotenfolge t
P_{k}  Raum der Polynome k-ter Ordnung
B  (B)_{i,j} = B_{j,k}(x_{i}) Beobachtungsmatrix
supp(.)  Support einer Funktion, suppf(x) := \overline{\{x , f(x) ≠ 0\}}
macheps  relative Maschinengenauigkeit

N (0, σ^{2}I) Normalverteilung mit Erwartungswert 0 und Kovarianzma-

trix σ^{2}I

#t_{i}  Vielfachheit des Knotens t_{i}
∥.∥  Euklidische Norm
K (x, ρ)  K (x, ρ) :={∥y − x∥ < ρ}
K̄ (x, ρ)  K̄ (x, ρ) :={∥y − x∥ ≤ ρ}

Tabelle 1: Verwendete Bezeichnungen

Inhaltsverzeichnis

1 Einleitung 1

2 B-Splines 3 2.1 Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Spezialisierung auf ein endliches Intervall . . . . . . . . . . . . . . . 4 2.3 Berechnung der Ableitung eines Splines . . . . . . . . . . . . . . . . 6

3 Lineare Quadratmittelprobleme 8

3.1 Nichtrestringierte lineare Quadratmittelprobleme . . . . . . . . . . . 8 3.1.1 Lösung von Problem LS mittels QR-Faktorisierung . . . . . . 9 3.1.2 Lineare Quadratmittelprobleme mit Bandstruktur . . . . . . 10

3.2 Lineare Quadratmittelprobleme mit linearen Ungleichheitsnebenbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2.1 Transformation von Problem LSI in Problem LDP . . . . . . 13

4 Nichtlineare Quadratmittelprobleme 16

4.1 Freie Quadratmittelprobleme . . . . . . . . . . . . . . . . . . . . . . 16 4.1.1 Das Gauß-Newton-Prinzip . . . . . . . . . . . . . . . . . . . . 16 4.1.2 Globalisierung . . . . . . . . . . . . . . . . . . . . . . . . . . 18

4.2 Restringierte Quadratmittelprobleme . . . . . . . . . . . . . . . . . . 20 4.2.1 Optimalitätsbedingungen . . . . . . . . . . . . . . . . . . . . 20 4.2.2 Verallgemeinertes Gauß-Newton-Verfahren . . . . . . . . . . . 21 4.2.3 Globalisierung . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2.4 Nichtlineare Quadratmittelprobleme mit linearen Ungleich-

heitsnebenbedingungen . . . . . . . . . . . . . . . . . . . . . 25 4.2.5 Regularisierung . . . . . . . . . . . . . . . . . . . . . . . . . . 27

5 Separable Quadratmittelprobleme 29 5.1 Problemformulierung . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.2 Gauß-Newton-ähnliche Verfahren für separable Quadratmittelprobleme 30 5.3 Der Übergang zum reduzierten Funktional . . . . . . . . . . . . . . . 32 5.4 Die Bestimmung der Jacobi-Matrix . . . . . . . . . . . . . . . . . . . 33 5.5 Konvergenzraten und Aufwand . . . . . . . . . . . . . . . . . . . . . 36

6 Splineapproximation mit freien Knoten 38 6.1 Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 6.2 Anordnungsbedingung für die Splineknoten . . . . . . . . . . . . . . 39 6.3 Untersuchung der Struktur des Problems . . . . . . . . . . . . . . . . 41

6.3.1 Der Übergang zum reduzierten Funktional beim freien Problem 41 6.3.2 Der Übergang zum reduzierten Funktional beim restringierten

Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

6.4 Existenz einer Lösung des reduzierten Problems . . . . . . . . . . . . 42 6.5 Die Berechnung der Jacobi-Matrix . . . . . . . . . . . . . . . . . . . 42

6.5.1 Die Berechnung der Jacobi-Matrix im Golub/Pereyra-Modell 43

6.5.2 Die Berechnung der Jacobi-Matrix im Kaufman-Modell … 43 6.6 Ein Verfahren zur Lösung des Problems . . . . . . . . . . . . . . . . 44 6.7 Testbeispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

6.7.1 Titan-Daten . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.7.2 Titan-Daten II . . . . . . . . . . . . . . . . . . . . . . . . . . 49 6.7.3 Spektraldaten . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

7 Splineglättung mit freien Knoten 56 7.1 Darstellung des Glättungsterms . . . . . . . . . . . . . . . . . . . . . 56 7.2 Approximation des Glättungsterms . . . . . . . . . . . . . . . . . . . 57 7.3 Formulierung des kompletten Funktionals . . . . . . . . . . . . . . . 58 7.4 Testbeispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

7.4.1 Titan-Daten III . . . . . . . . . . . . . . . . . . . . . . . . . . 60 7.4.2 RPN-14-Daten . . . . . . . . . . . . . . . . . . . . . . . . . . 62

Literatur 64

1 Einleitung

Gegeben seien die Meßpunkte

a ≤ x_{1} < … < x_{i} < x_{i+1} < … < x_{m} ≤ b,  (1.1)

und die mit Fehlern behafteten Meßwerte y_{i}, i = 1, …, m, einer unbekannten glatten Funktion f. Die Messungen seien stochastisch unabhängig und normalverteilt

y_{i} = f(x_{i}) + ϵ_{i},  i = 1, …, m,  (1.2)
ϵ = (ϵ_{1}, …, ϵ_{m})^{T} = N (0, σ^{2}I).  (1.3)

Die Funktion f soll durch eine Funktion s approximiert werden, welche

1. möglichst einfach zu berechnen ist,

2. durch wenige Parameter beschrieben werden kann,

3. ”möglichst wenig“ von den Meßwerten abweicht.

Da die Meßwerte fehlerbehaftet sind, ist es prinzipiell nicht sinnvoll, die Meßdaten zu interpolieren. Wir wollen die fehlerbehafteten Daten durch eine Summe von B-Splines approximieren. Wir machen dazu den Ansatz

ỹ_{i} = ∑_{j=1}^{n} B_{j,k}(x_{i})α_{j}.  (1.4)

Die Forderung nach einer möglichst geringen Abweichung lautet dann

y_{i} − ỹ_{i} ”klein“ i = 1, …, m.

Unter der Voraussetzung (1.3) ist ein, im stochastischen Sinne ”optimales“ Maß der Abweichung die Summe der Fehlerquadrate

∑_{i=1}^{m} [y_{i} − ỹ_{i}]^{2} ⟶ min,  (1.5)

∑_{i=1}^{m} [y_{i} − ∑_{j=1}^{n} B_{j,k}(x_{i})α_{j}]^{2} ⟶ {min \atop α} .  (1.6)

(1.6) liefert dann eine Maximum-Likelihood-Schätzung. Zur Definition von B-Splines benötigt man eine sogenannte Knotenfolge t. Die Festlegung der Knoten erfolgt meist interaktiv oder auf Grund einer heuristischen Strategie (siehe Prozedur NEWNOT in [dB78]). Ziel dieser Arbeit soll die Einbeziehung der Knotenfolge t (oder einer Teilmenge t̃ davon) in den Minimierungsprozeß sein

∑_{i=1}^{m} [y_{i} − ∑_{j=1}^{n} B_{j,k}(x_{i}, t̃)α_{j}]^{2} ⟶ min_{α,t̃} .  (1.7)

(1.7) ist ein hochgradig nichtlineares, restringiertes (da die Knoten t nicht beliebig gewählt werden können) Quadratmittelproblem spezieller Struktur.

Neben der Splineapproximation (1.6) betrachten wir in dieser Arbeit noch die Splineglättung

∑_{i=1}^{m} [y_{i} − ∑_{j=1}^{n} B_{j,k}(x_{i})α_{j}]^{2} + μ ∫_{a}^{b} [s^{(r)}(x)]^{2} dx ⟶ {min \atop α} .  (1.8)

Mit dem Glättungsparameter μ > 0 kann die Güte der Approximation s im Vergleich zur Glattheit gesteuert werden (siehe [Wah90] für den statistischen Hintergrund und einige Verallgemeinerungen). Die Einbeziehung der Knotenfolge in den Minimierungsprozeß (1.8) liefert ein ähnlich strukturiertes nichtlineares, restringiertes Quadratmittelproblem wie (1.7). Die Arbeit ist wie folgt aufgebaut: In Kapitel 2 werden B-Splines über eine Rekursionsformel eingeführt und einige später benötigte Eigenschaften dargestellt. Die Kapitel 3 - 5 enthalten die Grundlagen und verwendeten Algorithmen zur Lösung der auftretenden speziellen linearen und nichtlinearen Quadratmittelprobleme. Kapitel 6 beschäftigt sich mit der Splineapproximation mit freien Knoten, während in Kapitel 7 die Methoden auf die Splineglättung mit freien Knoten ausgedehnt werden. Das Listing des entwickelten Programmpaketes befindet sich in einem gesonderten Anhang.

2 B-Splines

Wir werden im folgenden B-Splines über eine Rekursionsformel definieren und einige später benötigte Eigenschaften angeben.

2.1 Definitionen

Definition 1 (Knotenfolge) Eine Folge t = (… t_{−1}, t_{0}, t_{1}, …) mit t_{j+1} ≥ t_{j} ∀j und lim_{j→±∞} t_{j} = ±∞ heißt (bi-

infinite) Knotenfolge.

Definition 2 (Normalisierte polynomiale B-Splines) Der j-te B-Spline der Ordnung k für die Knotenfolge t wird mit B_{j,k} bezeichnet und definiert durch:

B_{j,1}(t) := \{ \begin{matrix}1 & t_{j} ≤ t < t_{j+1} \\ 0 & sonst \end{matrix}
B_{j,k}(t) := ω_{j,k}(t) · B_{j,k−1}(t) + (1 − ω_{j+1,k}(t)) · B_{j+1,k−1}(t) für k > 1
mit ω_{j,k}(t) := \{ \begin{matrix}\frac{t−t_{j}}{t_{j+k−1}−t_{j}} & falls t_{j} < t_{j+k−1} \\ 0 & sonst. \end{matrix}

Für t_{j} = t_{j+k} erkennt man durch wiederholte Anwendung der Rekursionsformel sofort

B_{j,k}(t) = 0.

Es ist deshalb vernünftig,

t_{j+k} > t_{j} ∀j  (2.1)

zu fordern, um die trivialen B-Splines B_{j,k} ≡ 0 von vornherein von der Betrachtung auszuschließen. Wir geben jetzt nach [dB78] die wichtigsten Eigenschaften der B-Splines an.

Satz 1 (Eigenschaften der B-Splines) 1. suppB_{j,k} = [t_{j}, t_{j+k}] (falls t_{j+k} > t_{j})

2. B_{j,k}(t) > 0 ∀t ∈ (t_{j}, t_{j+k})

3. ∑_{j} B_{j,k}(t) = 1

4. B_{j,k_{|[t_{j} ,t_{j+1}]}} ∈ P_{k_{|[t_{j} ,t_{j+1}]}}

Die Glattheitseigenschaften von B-Splines hängen von der Anzahl der zusammenfallenden Knoten ab. Wir definieren daher:

Definition 3 (Knotenvielfachheit) Der Knoten t_{j} hat die Vielfachheit ν_{j}, wenn

∃i : t_{i} < t_{j} = t_{i+1} = … = t_{i+ν_{j}} < t_{i+ν_{j}+1}

Bezeichnung: ν_{j} = #t_{j}

Definition 4 (Spline,Splineraum) Die Menge S_{k,t} := {∑_{j} B_{j,k}α_{j} , α_{j} ∈ ℝ} heißt Splineraum der Ordnung k bezüglich der Knotenfolge t. Die Elemente des Splineraumes S_{k,t} heißen Splines.

Den angekündigten Zusammenhang zwischen Knotenvielfachheit und Glattheit sowie die lineare Unabhängigkeit der B-Splines zeigt das folgende Theorem:

Theorem 1 (Curry,Schoenberg, zitiert nach [dBH85]) Wenn t_{j} < t_{j+k} für alle j gilt, dann sind die B-Splines B_{j,k} linear unabhängig und bilden eine Basis für den Raum S̃_{k,t} aller stückweisen Polynome vom Grad k − 1 mit Bruchstellen (breakpoints) t_{j}, die in t_{j} (k − 1 − #t_{j})-mal stetig differenzierbar sind.

Bemerkung 1 Die Beziehung zwischen Glattheit und Knotenvielfachheit kann auch durch die folgende Formel ausgedrückt werden:

Anzahl der Glattheitsbedingungen in t_{j} + Knotenvielfachheit #t_{j} = k

Beispiel:

Ordnung k Knotenvielfachheit #t_{j} Glattheit in t_{j}
3  1  C^{1}
3  2  C
3  3  unstetig
3  4 (ausgeschlossen)
4  1  C^{2}
4  2  C^{1}
4  3  C
4  4  unstetig
4  5 (ausgeschlossen)

2.2 Spezialisierung auf ein endliches Intervall

Da jeder B-Spline nur einen endlichen Träger besitzt, genügt es, sich auf ein endliches Intervall zu beschränken. Wir betrachten daher im folgenden ausschließlich die (finite) Knotenfolge

t = (t_{1}, t_{2}, …, t_{n}, t_{n+1}, …, t_{n+k}).  (2.2)

Dabei bezeichnet n die Anzahl der B-Splines. Später werden wir B-Splines benutzen, um gegebene Daten (x_{i}, y_{i}), i = 1, …, m, zu approximieren. Die Meßstellen x_{i} befinden sich in einem endlichen Intervall x_{i} ∈ [a, b] ⊂ ℝ. Wir beschränken uns daher auf die Knotenfolge

t = (t_{1}, t_{2}, …, t_{n}, t_{n+1}, …, t_{n+k})

mit

t_{1} ≤ … ≤ t_{k} ≤ a b ≤ t_{n+1} ≤ … ≤ t_{n+k}

und t_{k+1} ≤ … ≤ t_{n} ∈ [a, b]. Die Knoten t_{k+1}, …, t_{n} werden als innere Knoten bezeichnet. Es gelte stets #t_{j} ≤ k ∀j. Stärkere Glattheitsforderungen (etwa einfache innere Knoten, #t_{j} = 1, j = k + 1, …, n) werden an späterer Stelle diskutiert. Für die äußeren Knoten t_{1}, …, t_{k} und t_{n+1}, …, t_{n+k} gelte einer Empfehlung von [Cox82] zufolge:

t_{1} = t_{2} = … = t_{k} = a  und  b = t_{n+1} = t_{n+2} = … = t_{n+k}  (2.3)

Bemerkung 2 Die Bedingung (2.3) ändert nur die Basis des Splineraumes, ist aber nach [Cox82] günstig für die Kondition der auftretenden Matrizen. Diese Bedingung ist für die weiteren Ausführungen nicht wesentlich.

Wir betrachten jetzt die Einschränkung von S_{k,t} auf [a, b].

S_{k,t_{|[a,b]}} = {∑_{j=1}^{n} B_{j,k}(t)α_{j} , α_{j} ∈ ℝ}

Statt S_{k,t_{|[a,b]}} werden wir im folgenden kurz S_{k,t} schreiben. Die Elemente dieses Raumes bezeichnen wir mit s. Wegen Satz 1 gibt es genau eine Darstellung

s(t) = ∑_{j=1}^{n} B_{j,k}(t)α_{j}  (2.4)

oder in Vektorschreibweise

s(t) = β(t)^{T} α
mit β(t) := (B_{1,k}(t), …, B_{n,k}(t))^{T}
und α := (α_{1}, …, α_{n})^{T} .

Die Zahlen α_{1}, …, α_{n} heißen Splinekoeffizienten.

Bemerkung 3 • Unter der Voraussetzung t_{j+k} > t_{j} ∀j ist die Dimension von S_{k,t} gleich n.

• Nach Definition 2 sind die B-Splines rechtsstetig. Am rechten Intervallende definieren wir daher B_{j,k}(b) := B_{j,k}(b−).

Bei der Quadratmittelapproximation mit B-Splines spielt die Struktur folgender Matrix, der sogenannten Beobachtungsmatrix, eine Rolle:

B = [ \begin{matrix}B_{1,k}(x_{1}) & … & B_{n,k}(x_{1}) \\ ⋮ & & ⋮ \\ B_{1,k}(x_{m}) & … & B_{n,k}(x_{m}) \end{matrix} ] ∈ ℝ^{m,n}

Satz 2 (Schoenberg-Whitney 1953) Die Matrix B ∈ ℝ^{m,n} hat Vollrang genau dann, wenn es eine Teilfolge der Datenabszissen (x_{i_{j}}) gibt, so daß

B_{j,k} (x_{i_{j}}) ≠ 0  j = 1, …, n.  (2.5)

Bemerkung 4 • Eine hinreichende Bedingung (und bei t_{j} < t_{j+k} auch notwendige) Bedingung für (2.5) ist die Schoenberg-Whitney-Bedingung

t_{j} < x_{i_{j}} < t_{j+k}.  (2.6)

• In der Originalarbeit von Schoenberg/Whitney[1953] wird der Satz für den Fall der Spline-Interpolation mit abgebrochenen Potenzfunktionen bewiesen. Die zitierte Version stammt aus [dB78, Theorem XIII.1, Lemma XIV.2].

• Die Matrix B ist sogar total positiv. Aus dieser Tatsache ergeben sich neben einigen wichtigen theoretischen Eigenschaften (variation diminishing property) Vereinfachungen für die numerische Implementierung (siehe [dBP77]).

• Wegen Eigenschaft 1, Satz 1 hat die Matrix B eine Bandbreite kleiner oder gleich k.

Beispiel: n = 8, m = 9, k = 4

x_{1} x_{2} x_{3}  x_{4} x_{5} x_{6}  x_{7} x_{8}  x_{9} \begin{matrix}t_{1} \\ t_{2} \\ t_{3} \\ t_{4} \end{matrix}  t_{5}  t_{6}  t_{7}  t_{8}  \begin{matrix}t_{9} \\ t_{10} \\ t_{11} \\ t_{12} \end{matrix}

B = [ \begin{matrix}x & x & x & x \\ x & x & x & x \\ x & x & x & x \\ & x & x & x & x \\ & & x & x & x & 0 \\ & & x & x & x & x \\ & & & x & x & x & x \\ & & & x & x & x & x \\ & & & & x & x & x & x \end{matrix} ] ∈ ℝ^{9,8}

2.3 Berechnung der Ableitung eines Splines

Wir wollen nun die Ableitung eines Splines berechnen.

Lemma 1 (Ableitung eines Splines) Wenn t_{j} < t_{j+k−1}, j = 2, …, n gilt, dann ist die Ableitung eines Splines aus S_{k,t} ein Spline der Ordnung k − 1 bezüglich derselben Knotenfolge, d.h.

DS_{k,t} = S_{k−1,t},

wobei D den Differentiationsoperator bezeichnet. Gilt s(t) = ∑_{j=1}^{n} B_{j,k}(t)α_{j}, so besitzt s′ die Darstellung

s′(t) = ∑_{j=2}^{n} B_{j,k−1}(t)α′_{j}  (2.7)

mit  α′_{j} := \frac{α_{j} − α_{j−1}}{\frac{t_{j+k−1}−t_{j}}{k−1}} .  (2.8)

Bemerkung 5 Gilt t_{j} = t_{j+k−1}, so ist B_{j,k}(t_{j}−) ≠ B_{j,k}(t_{j}+). Damit ist B_{j,k} in t_{j} nicht differenzierbar und somit auch s ∈ S_{k,t} als Linearkombination von B-Splines i.a. nicht differenzierbar.

Durch wiederholte Anwendung von Lemma 1 erhält man die r-te Ableitung eines Splines. Wir setzen t_{j} < t_{j+k−r}, j = r + 1, …, n mit r < k voraus. Nach Lemma 1 existiert dann s^{(r)}(t) und hat die Darstellung

s^{(r)}(t) = ∑_{j=r+1}^{n} B_{j,k−r}(t)α_{j}^{(r)}  (2.9)

mit rekursiv bestimmbaren Splinekoeffizienten α_{j}^{(r)}. Wir stellen diese Ableitung nach [Kun89] in Vektorform dar.

2 B-SPLINES 7 Bezeichnungen:

α^{(0)} := α ∈ ℝ^{n}  (2.10)
α^{(r)} := (α_{r+1}^{(r)}, …, α_{n}^{(r)})^{T} ∈ ℝ^{n−r}  (2.11)

Splinekoeffizienten von s^{(r)}

β_{r}(t) := (B_{r+1,k−r}(t), …, B_{n,k−r}(t))^{T} ∈ ℝ^{n−r}  (2.12)
H_{r} := (k − r) [ \begin{matrix}\frac{1}{t_{k+1}−t_{r+1}} & & & & 0 \\ & . & . & . \\ 0 & & & & \frac{1}{t_{n+k−r}−t_{n}} \end{matrix} ] ∈ ℝ^{n−r,n−r} (2.13)
L_{r} := [ \begin{matrix}−1 & 1 \\ & −1 & 1 \\ & & ... & ... \\ & & & ... & ... \\ & & & & −1 & 1 \\ & & & & & −1 & 1 \end{matrix} ] ∈ ℝ^{n−r,n−r+1}  (2.14)

Der Spline s^{(r)} besitzt dann die Darstelllung

s^{(r)}(t) = β_{r}(t)^{T} α^{(r)}.

Mit (2.13), (2.14) und (2.8) erhält man

α^{(r)} = H_{r}L_{r}α^{(r−1)}
= H_{r}L_{r} … H_{1}L_{1}α
= D_{r}α
mit  D_{r} := H_{r}L_{r} … H_{1}L_{1} ∈ ℝ^{n−r,n}.

In Vektorform erhalten wir also für die r-te Ableitung eines Splines

s^{(r)}(t) = β_{r}(t)^{T} D_{r}α.  (2.15)

3 Lineare Quadratmittelprobleme

3.1 Nichtrestringierte lineare Quadratmittelprobleme

Gegeben sei das lineare Quadratmittelproblem

Problem LS:

min {∥y − Ax∥ : x ∈ ℝ^{n}}  (3.1)

mit gegebenem A ∈ ℝ^{m,n}, b ∈ ℝ^{m} und m ≥ n. Zur Untersuchung des Problems LS betrachten wir die Singulärwertzerlegung von A.

Satz 3 (Singulärwertzerlegung) Sei A ∈ ℝ^{m,n}. Dann folgt: Es existieren orthogonale Matrizen U ∈ ℝ^{m,m}, V ∈ ℝ^{n,n} und eine Diagonalmatrix

Σ = ( \begin{matrix}σ_{1} \\ & & & 0 \\ & & σ_{r} \\ & 0 & & 0 \end{matrix} ) ∈ ℝ^{m,n} r = Rang(A)

mit σ_{1} ≥ … ≥ σ_{r} > 0, so daß A = UΣV^{T}. Die Zahlen σ_{1}, …, σ_{r} sind durch A eindeutig festgelegt und heißen Singulärwerte von A.

Die Anwendung der Singulärwertzerlegung auf das Problem LS ergibt:

r = b − Ax
= UU^{T} b − UΣV^{T} x

Unter Beachtung der Orthogonalinvarianz der Euklidischen Norm folgt:

∥r∥^{2} = ∥U^{T} b − ΣV^{T} x∥^{2}
= ∑_{i=1}^{r} (u_{i}^{T} b − σ_{i}v_{i}^{T} x)^{2} + ∑_{i=r+1}^{m} (u_{i}^{T} b)^{2}
U = [u_{1}, …, u_{m}] V = [v_{1}, …, v_{n}]

Der erste Summand ist stets nichtnegativ und verschwindet genau dann, wenn

x = ∑_{i=1}^{r} (\frac{u_{i}^{T} b}{σ_{i}}) v_{i}.

Definition 5 (Moore-Penrose-Inverse) Sei A = UΣV^{T} die Singulärwertzerlegung von A ∈ ℝ^{m,n} und Σ^{+} ∈ ℝ^{n,m} definiert durch

Σ^{+} := ( \frac{1}{σ_{1}} ; 0 ; \frac{1}{σ_{r}} ; 0  0 )

Dann heißt A^{+} := VΣ^{+}U^{T} ∈ ℝ^{n,m} Moore-Penrose-Inverse von A.

Damit läßt sich die eindeutige Minimum-Norm-Lösung x^{N} durch x^{N} = A^{+}b angeben. Wir geben jetzt noch eine wichtige Eigenschaft der Moore-Penrose-Inverse an.

Satz 4 (Moore-Penrose-Bedingungen) Die Matrix A^{+} ist durch die folgenden Bedingungen eindeutig festgelegt:

(P 1) AA^{+}A = A  (P 2) A^{+}AA^{+} = A^{+}
(P 3) (AA^{+})^{T} = AA^{+}  (P 4) (A^{+}A)^{T} = A^{+}A

Bemerkung 6 Matrizen, welche nur einen Teil der Moore-Penrose-Bedingungen erfüllen, werden als verallgemeinerte Inverse bezeichnet und durch die erfüllten Moore-Penrose-Bedingungen gekennzeichnet. Eine Matrix A^{−}, welche die Bedingungen (P1) und (P3) erfüllt, heißt auch symmetrische verallgemeinerte Inverse. Obwohl Moore-Penrose-Inversen sehr oft bei der theoretischen Untersuchung von Algorithmen benutzt werden, ist es oftmals billiger, bei der numerischen Implementierung mit verallgemeinerten Inversen zu arbeiten, welche nur die für den Algorithmus notwendigen Moore-Penrose-Bedingungen erfüllen (siehe auch [GP73] und [Kau75]).

Wir werden das Problem LS nicht mittels einer Singulärwertzerlegung, sondern mit einer (billigeren) QR-Faktorisierung lösen.

3.1.1 Lösung von Problem LS mittels QR-Faktorisierung

Alle in unseren Algorithmen auftretenden Matrizen haben vollen Rang r = n, was wir nötigenfalls durch Regularisierung in den entsprechenden nichtlinearen Quadratmittelproblemen erreichen können. Wir spezialisieren daher die QR-Faktorisierung auf diesen Fall. Sei A ∈ ℝ^{m,n}, m ≥ n eine Matrix mit Rang r = n. Dann existieren

Q ∈ ℝ^{m,m} orthogonal
Q = (Q_{1}, Q_{2})  Q_{1} ∈ ℝ^{m,n} Q_{2} ∈ ℝ^{m,m−n}
R̂ = ( \begin{matrix}R \\ 0 \end{matrix} )  R̂ ∈ ℝ^{m,n}, R ∈ ℝ^{n,n} reguläre obere Dreiecksmatrix
P ∈ ℝ^{n,n} Permutationsmatrix,

so daß

Q^{T} AP = R̂
[ \begin{matrix}Q_{1}^{T} \\ Q_{2}^{T} \end{matrix} ] AP = [ \begin{matrix}R \\ 0 \end{matrix} ]

Bemerkung 7 • Im Vollrangfall ist die Permutationsmatrix P eigentlich nicht nötig. Sie wird jedoch zur Verbesserung der numerischen Stabilität der QR-Faktorisierung verwendet.

• Für vollbesetzte Matrizen A ∈ ℝ^{m,n} erzeugen wir Q durch eine Folge von Householder-Transformationen. Wir verwenden einen Algorithmus mit Spaltenpivotisierung (siehe [GvL89], [KS88]). Die Implementierung lehnt sich an die Prozedur QRDC aus LINPACK [DBMS79] an.

Unter Benutzung der QR-Faktorisierung von A können wir nun einen Algorithmus zur Lösung von Problem LS angeben. Es gilt:

∥r∥^{2} = ∥b − Ax∥^{2} = ∥Q^{T} (b − Ax)∥^{2}
= ∥Q^{T} b − (Q^{T} AP) (P^{T} x)∥^{2}
= ∥( \begin{matrix}Q_{1}^{T} b \\ Q_{2}^{T} b \end{matrix} ) − [ \begin{matrix}R \\ 0 \end{matrix} ] (P^{T} x)∥^{2}

Wir definieren y := P^{T} x ∈ ℝ^{n}, c_{1} := Q_{1}^{T} b ∈ ℝ^{n}, c_{2} := Q_{2}^{T} b ∈ ℝ^{m−n} und erhalten

∥r∥^{2} = ∥c_{1} − Ry∥^{2} + ∥c_{2}∥^{2}

y können wir nun aus dem oberen Dreieckssystem Ry = c_{1} bestimmen. Besitzt die zu transformierende Matrix eine spezielle Gestalt, so lassen sich Householder-Spiegelungen nur unter besonderen Vorkehrungen effektiv einsetzen (fill in). Wir verwenden aus diesem Grund (explizite) Givens-Drehungen. Zudem besitzen die Givens-Drehungen bei der Faktorisierung einer Matrix die deutlich besseren Fehlerkumulationskonstanten (siehe [KS88]).

3.1.2 Lineare Quadratmittelprobleme mit Bandstruktur

Bei der Approximation bzw. beim Glätten mit B-Splines treten Matrizen der folgenden Form auf:

(a) ( \begin{matrix}x & x & x & x \\ x & x & x & x \\ x & x & x & x \\ & x & x & x & x \\ & x & x & x & x \\ & x & x & x & x \\ & x & x & x & x \\ & & x & x & x & x \\ & & x & x & x & x \\ & & x & x & x & x \end{matrix} ) (b) ( \begin{matrix}x & x & x & x \\ & x & x & x & x \\ & & x & x & x & x \\ & & & x & x & x \\ & & & & x & x \\ & & & & & x \\ x & x & x \\ & x & x & x \\ & & x & x & x \\ & & & x & x & x \end{matrix} )

Für eine Matrix A ∈ ℝ^{m,n} definieren wir

f_{i} := min{j : a_{i,j} ≠ 0}
l_{i} := max{j : a_{i,j} ≠ 0}.

Die Zahlen f_{i} und l_{i} sind die Indizes des ersten bzw. letzten Nichtnullelementes in der i-ten Zeile.

Definition 6 Die Matrix A ∈ ℝ^{m,n} hat die (Zeilen) Bandbreite nbands, wenn

nbands = max_{i=1,...,m} w_{i} w_{i} := (l_{i} − f_{i} + 1)

Die Matrizen in (a) und (b) sind beide von der Bandbreite 4.

Definition 7 Eine Matrix ist in Standardform, wenn der Index des am weitesten rechts stehenden Nichtnullelementes l_{i} jeder Zeile eine nichtfallende Funktion der Zeilennummer ist.

Bemerkung 8 • Die Matrix in (a) ist offenbar in Standardform, während (b) nicht in Standardform ist.

• Durch Zeilenvertauschung kann jede Matrix in Standardform gebracht werden. Wir erhalten z.B. aus (b)

( \begin{matrix}x & x & x \\ & x & x & x \\ x & x & x & x \\ & & x & x & x \\ & x & x & x & x \\ & & & x & x & x \\ & & x & x & x & x \\ & & & x & x & x \\ & & & & x & x \\ & & & & & x \end{matrix} )

eine Matrix in Standardform, indem wir, bis auf die ersten und letzten Zeilen, abwechselnd je eine Zeile der Teilmatrizen nehmen.

Betrachten wir nun das Problem LS

min {∥y − Ax∥ : x ∈ ℝ^{n}}

mit einer Bandmatrix A ∈ ℝ^{m,n} der Bandbreite nbands. A habe Vollrang n.

Grundalgorithmus

S1: Transformation des m × n Problems in ein äquivalentes n × n Problem

S1.1: Berechnung der Givens-Faktorisierung der Matrix A

[ \begin{matrix}Q_{1}^{T} \\ Q_{2}^{T} \end{matrix} ] A = [ \begin{matrix}R \\ 0 \end{matrix} ]  R ∈ ℝ^{n,n}  reguläre obere Dreiecksmatrix (Q_{1}, Q_{2}) ∈ ℝ^{m,m}  orthogonal

S1.2: Transformation der rechten Seite b

c_{1} = Q_{1}^{T} b c_{1} ∈ ℝ^{n}

S2: Lösung des Gleichungssystems Rx = c_{1} mit regulärer oberer Dreiecksmatrix R

x = R^{−1}c_{1}

Bemerkung 9 • Die Schritte S1.1 und S1.2 können bei Bedarf zusammen realisiert werden. Die Speicherung und spätere Rekonstruktion der Drehungsparameter kann dann entfallen.

• Wenn A eine Blockbandmatrix in Standardform der Bandbreite nbands ist, hat die obere Dreiecksmatrix R diesselbe Bandbreite ([Cox81]).

• Bei vorliegender Rangdefizienz muß Schritt S2 geeignet modifiziert werden (siehe [Cox89]). Im Falle total positiver Matrizen, wie sie bei der Approximation mit B-Splines auftreten ([Cox82]), ist keine Pivotisierung in Schritt S1.1 nötig (siehe [dBP77]).

Zur Berechnung der Givens-Faktorisierung der Matrix A verwenden wir einen Algorithmus zur zeilenweisen Givens-Faktorisierung nach [Cox81].

Algorithmus zur zeilenweise Givens-Faktorisierung Arbeitsvektor (u, v) zur Speicherung der i-ten Zeile von (A, b)

S1: Initialisierung

Setze R ∈ ℝ^{n,n}, c ∈ ℝ^{n}, ∥r∥^{2} Null

S2: Berechnung

for i:=1 to m do

S2.1: Stelle (a_{i}^{T} , b_{i}) auf (u, v) bereit

S2.2: { Annuliere u durch Givens-Drehungen mit den Zeilen von R }

for j:=f_{i} to l_{i} do

if u_{j} ≠ 0 then

S2.2.1: setup GIVENS rotation( r_{jj},u_{j})

S2.2.2: apply GIVENS rotation(nbands,r_{jj},u_{j})

S2.2.3: apply GIVENS rotation(1,c_{j},v)

S2.3: Aktualisiere das Residium ∥r∥^{2} := ∥r∥^{2} + v^{2}

Bemerkung 10 • Wenn A ∈ ℝ^{m,n} eine Bandmatrix der Bandbreite nbands ist und das Problem für mehrere rechte Seiten b zu lösen ist, so müssen in Schritt S2.2.1 zusätzlich Informationen über die Drehungsparameter gespeichert werden. Diese können in einer Matrix Q der Dimension m×nbands aufgehoben werden. Die Schritte S2.2.3 und S2.3 können dann entfallen.

• Der Aufwand beträgt nach [Cox81] (siehe auch [Bjö 90]) für eine Matrix A ∈ ℝ^{m,n}

Matrix-Struktur Standardform Nicht-Standardform
volle Matrix ∼ 2m n^{2} ∼ 2m n^{2}
Bandmatrix  ∼ 2m nbands^{2} ∼ 2m n nbands

flops.

• Die Wirkungsweise des Algorithmus findet man auch in [Cox81] und [Kun89] beschrieben.

3.2 Lineare Quadratmittelprobleme mit linearen Ungleichheitsnebenbedingungen

Bei der Lösung der später betrachteten nichtlinearen Quadratmittelprobleme muß eine Folge von linearen Quadratmittelproblemen mit linearen Ungleichheitsnebenbedingungen gelöst werden. Wir betrachten daher jetzt

Problem LSI:

min {∥Ax − b∥ : Cx ≥ h, x ∈ ℝ^{n}}  (3.2)

Problem LDP (Least Distance Programming):

min {∥x∥ : Cx ≥ h, x ∈ ℝ^{n}}  (3.3)

Problem NNLS (Nonnegative Least Squares):

min {∥Ax − b∥ : x ≥ 0, x ∈ ℝ^{n}}  (3.4)

mit Matrizen A ∈ ℝ^{m,n}, C ∈ ℝ^{ncstr,n} und Vektoren h ∈ ℝ^{ncstr}, b ∈ ℝ^{m}. Das allgemeine Problem LSI mit einer Vollrangmatrix A wird durch eine orthogonale Transformation in das Problem LDP überführt. Das Problem LDP schließlich wird in das Problem NNLS transformiert, welches mit einer aktiven Mengenstrategie gelöst wird (siehe [LH74, Chapter 23]). Die Prozeduren zur Lösung von LDP und NNLS entnehmen wir ebenfalls den gleichnamigen FORTRAN-Programmen in [LH74]. Wir geben daher nur die Transformation von Problem LSI in Problem LDP an.

3.2.1 Transformation von Problem LSI in Problem LDP

Sei A ∈ ℝ^{m,n} eine Matrix mit Vollrang n, C ∈ ℝ^{ncstr,n}, b ∈ ℝ^{m} und h ∈ ℝ^{ncstr}.

Problem LSI:

min {∥Ax − b∥ : Cx ≥ h, x ∈ ℝ^{n}}  (3.5)

Sei AP = QR̂ die QR-Faktorisierung von A, wie in Abschnitt 3.1.1 beschrieben.

Q ∈ ℝ^{m,m} orthogonal
Q = (Q_{1}, Q_{2})  Q_{1} ∈ ℝ^{m,n} Q_{2} ∈ ℝ^{m,m−n}
R̂ = ( \begin{matrix}R \\ 0 \end{matrix} )  R̂ ∈ ℝ^{m,n}, R ∈ ℝ^{n,n} reguläre obere Dreiecksmatrix
P ∈ ℝ^{n,n} Permutationsmatrix.

Dann gilt

Q^{T} AP = R̂
[ \begin{matrix}Q_{1}^{T} \\ Q_{2}^{T} \end{matrix} ] AP = [ \begin{matrix}R \\ 0 \end{matrix} ]

Wir erhalten

φ(x) := \frac{1}{2} ∥Ax − b∥^{2} = \frac{1}{2} ∥(Q^{T} AP) (P^{T} x) − Q^{T} b∥^{2}
= \frac{1}{2} ∥[ \begin{matrix}R \\ 0 \end{matrix} ] (P^{T} x) − ( \begin{matrix}Q_{1}^{T} b \\ Q_{2}^{T} b \end{matrix} )∥^{2}

Wir definieren wieder y := P^{T} x ∈ ℝ^{n}, c_{1} := Q_{1}^{T} b ∈ ℝ^{n}, c_{2} := Q_{2}^{T} b ∈ ℝ^{m−n} und erhalten

φ(x) = \frac{1}{2} ∥Ry − c_{1}∥^{2} + \frac{1}{2} ∥c_{2}∥^{2} .

Wir führen eine erneute Variablentransformation aus

z := Ry − c_{1} ∈ ℝ^{n}.

Damit erhalten wir aus

φ(x) := \frac{1}{2} ∥Ax − b∥^{2} ⟶ Min!
φ̃(z) := \frac{1}{2} ∥z∥^{2} + \frac{1}{2} ∥c_{2}∥^{2} ⟶ Min!

Transformation der Nebenbedingungen

Cx ≥ h
y = P^{T} x ⇝ CPy ≥ h
z = Ry − c_{1} ⇝ CPR^{−1}(z + c_{1}) ≥ h
CPR^{−1}z ≥ h − CPR^{−1}c_{1}

Das Problem LSI ist damit äquivalent zu folgendem Least Distance Problem

∥z∥ ⟶ Min!
mit CPR^{−1}z ≥ h − CPR^{−1}c_{1}

Zusammenfassend geben wir noch einmal den Algorithmus zur Lösung von Problem LSI an.

Algorithmus LSI { Lösung des Problems ∥Ax − b∥ ⟶ Min! mit Cx ≥ h mit Vollrangmatrix A }

S1: Berechne QR-Faktorisierung von A

[ \begin{matrix}Q_{1}^{T} \\ Q_{2}^{T} \end{matrix} ] A = [ \begin{matrix}R \\ 0 \end{matrix} ]

S2: Berechne c_{1} := Q_{1}^{T} b ∈ ℝ^{n}

S3: Berechne

C̃ := CPR^{−1} ∈ ℝ^{ncstr,n}
h̃ := h − CPR^{−1}c_{1} ∈ ℝ^{ncstr}

S4: Löse das Least Distance Problem

∥z∥ ⟶ Min!
mit CPR^{−1}z ≥ h − CPR^{−1}c_{1}

S5: Bilde z + Q_{1}^{T} b = z + c_{1}

S6: Löse das obere Dreieckssystem

Ry = z + c_{1}

S7: Bilde x := Py

Bemerkung 11 • Die QR-Faktorisierung in Schritt S1 erfolgt mit Householder-Spiegelungen mit Spaltenpivotisierung.

• Vor der Lösung der oberen Dreickssysteme in S3 und S6 erfolgt eine Schätzung der Konditionszahl von R. Der verwendete Algorithmus folgt Algorithmus A3.3.1 CONDEST in [DS83] (siehe auch [CMSW79]). Wenn cond(R) > \frac{1}{\sqrt{macheps}} bricht der Algorithmus ab. Die Matrix A ist dann (fast) rangde-

fizient. Der aufrufende Algorithmus zur Lösung nichtlinearer Quadratmittelprobleme kann dann Maßnahmen zur Regularisierung (Levenberg-Marquardt o.ä.) treffen.

• In [Bjö 90] wird die Transformation von Problem LSI mit rangdefizienter Matrix A in ein verallgemeinertes Least Distance Problem betrachtet. Diese Transformation beruht auf einer kompletten orthogonalen Faktorisierung

Q^{T} APV = [ \begin{matrix}T & 0 \\ 0 & 0 \end{matrix} ]

T ∈ ℝ^{r,r} obere reguläre Dreiecksmatrix, V ∈ ℝ^{n,n} orthogonal (siehe [GvL89]) und führt auf ein verallgemeinertes Least Distance Problem

∥x_{1}∥ ⟶ Min!
mit C_{1}x_{1} + C_{2}x_{2} ≥ h.

4 Nichtlineare Quadratmittelprobleme

Die in dieser Arbeit auftretenden Optimierungsprobleme sind entweder freie Quadratmittelprobleme

min \{\frac{1}{2} ∥F(x)∥^{2} : x ∈ ℝ^{n}\}  (4.1)

mit F : ℝ^{n} ⟶ ℝ^{m}, m ≥ n oder restringierte Quadratmittelprobleme

min \{\frac{1}{2} ∥F_{1}(x)∥^{2} : F_{2}(x) = 0, F_{3}(x) ≥ 0, x ∈ ℝ^{n}\}  (4.2)

mit F_{i} : ℝ^{n} ⟶ ℝ^{m_{i}}. Wir können uns dabei hier sogar auf nichtlineare Quadratmittelprobleme mit linearen Ungleichheitsnebenbedingungen beschränken.

min \{\frac{1}{2} ∥F(x)∥^{2} : Cx ≥ h, x ∈ ℝ^{n}\}  (4.3)

mit F : ℝ^{n} ⟶ ℝ^{m}, C ∈ ℝ^{ncstr,n}, h ∈ ℝ^{ncstr}. Die zur Lösung solcher nichtlinearer Quadratmittelprobleme verwendeteten Algorithmen sind grundsätzlich iterativ und von der Form

x^{k+1} := x^{k} + s^{k},

wobei x^{k} die aktuelle Iterierte bezeichnet und s^{k} der Schritt (step) zur nächsten Iterierten x^{k+1} ist.

4.1 Freie Quadratmittelprobleme

4.1.1 Das Gauß-Newton-Prinzip

Gegeben sei eine Funktion F : ℝ^{n} ⟶ ℝ^{m}, F(x) = (F_{1}(x), …, F_{m}(x))^{T} mit m ≥ n. Wir betrachten das nichtlineare Quadratmittelproblem

min {∑_{i=1}^{m} [F_{i}(x)]^{2} : x ∈ ℝ^{n}}  (4.4)

Der einfacheren Beschreibung wegen definieren wir durch die Festlegung

g(x) := \frac{1}{2} ∑_{i=1}^{m} [F_{i}(x)]^{2} = \frac{1}{2} ∥F(x)∥^{2} ∀ x ∈ ℝ^{n}  (4.5)

eine Funktion g : ℝ^{n} ⟶ ℝ, mit der (4.4) in der Form

min {g(x) : x ∈ ℝ^{n}}  (4.6)

geschrieben werden kann. Falls die Funktionen F_{i}, i = 1, …, m, genügend glatt sind, erhalten wir

∇g(x) = F′(x)^{T} F(x)
∇^{2}g(x) = F′(x)^{T} F′(x) + S(x)
mit S(x) := F(x) ∘ F′′(x) = ∑_{i=1}^{m} F_{i}(x)F′′_{i}(x).

Wir bilden im Punkt x = x^{k} eine quadratische Modellfunktion

g(x + s) ≈ μ_{N}(x + s) := g(x) + ∇g(x)^{T} s + \frac{1}{2}s^{T} ∇^{2}g(x)s.

Eine naheliegende Möglichkeit zur Lösung von (4.6) besteht darin, in jedem Schritt s^{k} als Lösung des quadratischen Ersatzproblems

min \{μ_{N}(x + s) = g(x) + ∇g(x)^{T} s + \frac{1}{2}s^{T} ∇^{2}g(x)s : s ∈ ℝ^{n}\}  (4.7)

zu bestimmen. Dies führt auf das Newton-Verfahren zur Lösung von (4.6). Aus Aufwandsgründen muß jedoch die Berechnung aller partiellen Ableitungen zweiter Ordnung \frac{∂^{2}F_{i}}{∂x_{j}∂x_{k}} i = 1, …, m, j = 1, …, n, k = 1, … n entfallen.

Eine weitere Möglichkeit besteht darin, den Term S(x) in der Hesse-Matrix ∇^{2}g(x) zu vernachlässigen, da dieser Term bei Problemen mit kleinen Residuen bzw. bei schwach nichtlinearen Problemen klein ist. Dieses Vorgehen führt auf das Gauß-Newton-Modell

min \{μ_{GN}(x + s) = g(x) + ∇g(x)^{T} s + \frac{1}{2}s^{T} H(x)s : s ∈ ℝ^{n}\}  (4.8)

mit H := F′(x)^{T} F′(x). Es ist äquivalent mit einer Ersetzung der Funktion F im Punkt x^{k} durch eine linearisierte Funktion

F̃_{k} := F(x^{k}) + F′(x^{k}) (x − x^{k}).

Äquivalent zu (4.8) ist

min \{μ_{GN}(x + s) = \frac{1}{2} ∥F + F′s∥^{2} : s ∈ ℝ^{n}\} .  (4.9)

Wenn F′(x^{k}) Vollrang hat, so besitzt das Problem (4.9) die eindeutige Lösung

s^{k} := − [F′(x^{k})^{T} F′(x^{k})]^{−1} F′(x^{k})^{T} F(x^{k})  (4.10)

und wir erhalten die Iterationsvorschrift

x^{k+1} := x^{k} + s^{k}.  (4.11)

Diese Iterationsvorschrift bezeichnet man als Gauß-Newton-Verfahren. Aus den zahlreichen Konvergenzsätzen für das Gauß-Newton-Verfahren zitieren wir eine affin-invariante Version aus [DH91].

Satz 5 (Lokale Konvergenz des Gauß-Newton-Verfahrens) Sei D ⊂ ℝ^{n} offen, konvex und F : D ⊂ ℝ^{n} ⟶ ℝ^{m}, m ≥ n, eine stetig differenzierbare Abbildung, deren Jacobi-Matrix F′(x) für alle x ∈ D vollen Rang n habe. Es existiere eine Lösung x^{∗} ∈ D des zugehörigen nichtlinearen Quadratmittelproblems ∥F(x)∥ ⟶ min. Ferner gebe es Konstanten ω > 0 und 0 ≤ κ^{∗} < 1, so daß

∥F′(x)^{+} [F′(x + sv) − F′(x)] v∥ ≤ sω ∥v∥^{2}

für alle s ∈ [0, 1], x ∈ D und v ∈ ℝ^{n} mit x + v ∈ D, und

∥[F′(x)^{+} − F′(x^{∗})^{+}] F′(x^{∗})∥ ≤ κ^{∗} ∥x − x^{∗}∥

für alle x ∈ D. Falls dann für einen gegebenen Startwert x^{0} ∈ D

ρ := ∥x^{0} − x^{∗}∥ < 2 \frac{1 − κ^{∗}}{ω} =: σ

gilt, so bleibt die durch das Gauß-Newton-Verfahren (4.10), (4.11) definierte Folge {x^{k}} in der offenen Kugel K (x^{∗}, ρ) und konvergiert gegen x^{∗}, d.h.

∥x^{k} − x^{∗}∥ < ρ für k > 0 und lim_{k→∞} x^{k} = x^{∗}.

Die Konvergenzgeschwindigkeit läßt sich abschätzen durch

∥x^{k+1} − x^{∗}∥ ≤ \frac{ω}{2} ∥x^{k} − x^{∗}∥^{2} + κ^{∗} ∥x^{k} − x^{∗}∥ .

Insbesondere ergibt sich quadratische Konvergenz für kompatible nichtlineare Quadratmittelprobleme ( ∥F(x^{∗}∥ = 0). Darüber hinaus ist die Lösung x^{∗} eindeutig in K (x^{∗}, ρ).

4.1.2 Globalisierung

Wir verwenden in unseren Algorithmen eine global konvergente Modifikation des Gauß-Newton-Verfahrens. Wir führen ein gedämpftes Gauß-Newton-Verfahren mit Armijo-Goldstein-Schrittweitenstrategie aus (siehe [Sch79]).

Gedämpftes Gauß-Newton-Verfahren mit akzeptiertem Schätzwert

S1: Wähle Startwert x^{0} ∈ ℝ^{n}, α ∈ (0, 1), δ ∈ (0, \frac{1}{2})

Setze k := 0; γ_{−1} := 1;

S2: Iteration

Repeat

S2.1: F := F(x^{k}); Berechnung des Residiums

J := F′(x^{k}); Berechnung der Jacobi-Matrix S2.2: Bestimme Suchrichtung p^{k}(= −s^{k}) aus

Jp^{k} ≅ F

S2.3: γ̃_{k} := γ_{k}; x^{+} := x^{k} − γ̃_{k}p^{k}; F^{+} := F(x^{+}); S2.4: Abstiegstest

if \frac{1}{2} ∥F∥^{2} − \frac{1}{2} ∥F^{+}∥^{2} ≥ δγ̃_{k}F^{T} Jp^{k}

then

x^{k+1} := x^{+};
γ_{k+1} := min(γ_{k}/α, 1);
k := k + 1;

else

γ_{k} := γ_{k} ∗ α; goto S2.3;

until (k > kmax) or Konvergenz

Bemerkung 12 • Die Berechnung der Suchrichtung p^{k} in Schritt S2.2 kann über die Normalgleichungen (J^{T} J)p^{k} = J^{T} F mittels Cholesky-Zerlegung erfolgen. Wir wählen jedoch die direkte Lösung mittels QR-Faktorisierung der Matrix J. Dazu verwenden wir Householder-Spiegelungen mit Spaltenpivotisierung (siehe Abschnitt 3.1.1).

• Die analytische Berechnung der Jacobi-Matrix J gestaltet sich wegen der komplizierten Gestalt der später untersuchten Funktionen F äußerst schwierig. Wir gehen deshalb einen anderen Weg und verwenden eine Approximation A_{k} der Jacobi-Matrix. Die einfachste und nächstliegende Möglichkeit besteht darin, die j-te Spalte von A_{k} als Differenzenquotient darzustellen.

A_{k}(x^{k})e^{j} = \frac{F(x^{k} + h_{j}^{k}e^{j}) − F(x^{k})}{h_{j}^{k}}  (4.12)

Die Berechnung einer Jacobi-Matrix kostet n + 1 Funktionsauswertungen von F. Wir benutzen diese Approximation mit h_{j}^{k} := ε_{1}^{d} ∗ {|x_{j}^{k}| + ε_{2}^{d}}. Wir wählen dabei ε_{1}^{d} in der Größenordnung von \sqrt{macheps} und ε_{2}^{d} ≈ 10^{−10}.

In einigen späteren Anwendungen (siehe Abschnitt 5) wird diese Strategie innerhalb eines bekannten analytischen Ausdrucks für J angewendet. Konvergenzaussagen für sogenannte konsistente Approximationen, d.h.

lim_{k→∞} ∥A_{k} − F′(x^{k})∥ = 0,

findet man in [Sch79].

• Das Abbruchkriterium konkretisieren wir wie folgt:

Repeat … until

(k > kmax) or (∥F(x^{k})∥ ≤ ε_{0}^{a}) or (∥x^{k+1} − x^{k}∥ ≤ ε_{1}^{a} (∥x^{k}∥ + ε_{2}^{a})) or

(∥F′(x^{k})^{T} F(x^{k})∥ ≤ ε_{3}^{a}) or (|F(x^{k})^{T} F′(x^{k})s| ≤ ε_{4}^{a})

Das Verfahren bricht ab, wenn die maximale Schrittanzahl kmax erreicht ist. Die Defektschranke ε_{0}^{a} wird vom Nutzer vorgegeben.

\frac{∥x^{k+1}−x^{k}∥}{∥x^{k}∥} gibt für Q-superlinear konvergente Verfahren den relativen Fehler ε_{1}^{a} der Iterierten x^{k} asymptotisch richtig wieder. ε_{2}^{a} macht den Test auch im Fall ∥x^{k}∥ = 0 sinnvoll und wird etwa in der Größenordnung 10^{−3} gewählt.

Der Test (∥F′(x^{k})^{T} F(x^{k})∥ ≤ ε_{3}^{a}) stellt fest, ob ein stationärer Punkt erreicht

ist. (|F(x^{k})^{T} F′(x^{k})s| ≤ ε_{4}^{a}) stellt schließlich fest, ob überhaupt noch eine Abstiegsrichtung s gefunden werden kann.

Die Defektschranke ε_{0}^{a}, die Schranke ε_{3}^{a} für die Norm des Gradienten sowie die Schranke ε_{4}^{a} sind abhängig von der Größenordnung der Funktionswerte (und von der Skalierung von x) und können vom Nutzer vorgegeben werden.

• Der Abstiegstest sichert einen hinreichenden Abstieg der Funktion g(x) = \frac{1}{2} ∥F(x)∥^{2} vom Punkt x^{k} zum Punkt x^{+} = x^{k} − γ̃_{k}p^{k}. Der Schrittweitenparameter γ̃_{k} wird so gewählt, daß die Armijo-Goldstein-Bedingung

g(x^{k}) − g(x^{k} − γ̃_{k}p^{k}) ≥ δγ̃_{k}∇g(x^{k})^{T} p^{k}

erfüllt ist. Die Wahl des Verfahrensparameter δ ∈ (0, \frac{1}{2}) sichert, daß die Minimumstelle einer quadratischen Funktion nicht ausgeschlossen wird. Die Bestimmung eines zulässigen γ̃_{k} kann durch eine Kombination einer Reduktionsund einer Interpolationsformel (etwa durch eine quadratische oder kubische Parabel) effizienter erfolgen. Wir haben uns in unseren Algorithmen auf die einfache Reduktionstechnik beschränkt.

4.2 Restringierte Quadratmittelprobleme

Wir betrachten jetzt das restringierte nichtlineare Quadratmittelproblem (NP)

min \{\frac{1}{2} ∥F_{1}(x)∥^{2} : F_{2}(x) = 0, F_{3}(x) ≥ 0, x ∈ ℝ^{n}\}  (4.13)

mit F_{i} : D ⊂ ℝ^{n} ⟶ ℝ^{m_{i}}. J_{i} bezeichne die entsprechenden Jacobi-Matrizen F′_{i}. Zur Charakterisierung einer Optimallösung von (4.13) zitieren wir Standardaussagen der nichtlinearen Optimierung. Unsere Darstellung folgt [Boc87] (siehe auch [Fle87], [GMW81], [Spe93]).

4.2.1 Optimalitätsbedingungen

Definition 8 (Lagrange-Funktion) Die Funktion L : ℝ^{n} × ℝ^{m_{2}} × ℝ_{+}^{m_{3}} ⟶ ℝ mit

L(x, λ, μ) := \frac{1}{2} ∥F_{1}(x)∥^{2} − λ^{T} F_{3}(x) − μ^{T} F_{3}(x)

x ∈ ℝ^{n}, λ ∈ ℝ^{m_{2}}, μ ∈ ℝ_{+}^{m_{3}} heißt Lagrange-Funktion zu (4.13).

Wir fassen die an einer Stelle z ∈ D aktiven Ungleichheitsnebenbedingungen durch die Indexmenge I(z) zusammen:

I(z) := {i | F_{3i}(z) = 0, i = 1, …, m_{3}}

Die entsprechenden Komponenten von F_{3} und ihre Ableitungen seien mit F̃_{3} bzw.

J̃_{3} bezeichnet. Ferner bezeichne F_{c} := ( F_{2} F̃_{3} ) die zusammengesetzte Abbildung

F_{c} : D ⊂ ℝ^{n} ⟶ ℝ^{m_{c}}, m_{c} := m_{2} + m̃_{3}, m̃_{3} := #I und J_{c} := [ J_{2} J̃_{3} ] die entspre-

chende Jacobi-Matrix.

Definition 9 (Regulärer Punkt) Der Punkt z ∈ D heißt regulär, wenn

Rang J_{c}(z) = m_{c}.  (4.14)

Bemerkung 13 Die in der Definition geforderte Regularitätsbedingung (auch Slater-Bedingung, constraint qualification genannt) kann abgeschwächt werden. Es muß gefordert werden, daß der Tangentenkegel F an einer Stelle z mit der Menge F aller zulässigen Richtungen der linearisierten Nebenbedingungen übereinstimmt (siehe [Fle87, S.202]). Hinreichende Bedingungen dafür sind

• Die Vektoren J_{2i}, i = 1, …, m_{2} und J̃_{3i}, i = 1, …, m̃_{3} sind linear unabhängig ( ≡ 4.14).

• Die aktiven Nebenbedingungen F_{2} = 0 F̃_{3} = 0 sind linear.

Einen Überblick über gebräuchliche Restriktionsqualifikationen der nichtlinearen Optimierung gibt [Spe93].

Wir können nun die notwendigen Optimalitätsbedingungen erster Ordnung formulieren.

Satz 6 (Notwendige Optimalitätsbedingungen erster Ordnung) Sei x^{∗} ∈ D ein regulärer Punkt und Lösung des restringierten nichtlinearen Quadratmittelproblems (4.13). Dann ist x^{∗} zulässig

F_{2}(x^{∗}) = 0 F_{3}(x^{∗}) ≥ 0  (4.15)

und es gilt: Es existieren Lagrange-Parameter λ^{∗}, μ^{∗} als eindeutige Lösung der Stationaritätsbedingungen

\frac{∂}{∂x}L(x^{∗}, λ^{∗}, μ^{∗}) = F_{1}(x^{∗})^{T} J_{1}(x^{∗}) − λ^{∗}^{T} J_{2}(x^{∗}) − μ^{∗}^{T} J_{3}(x^{∗}) = 0  (4.16)
μ^{∗} ≥ 0  i ∉ I(x^{∗}) ⟹ μ_{i}^{∗} = 0  (4.17)

Ein Tupel (x^{∗}, λ^{∗}, μ^{∗}), das die Zulässigkeitsbedingung (4.15) und die Stationaritätsbedingungen (4.16), (4.17) erfüllt, heißt Kuhn-Tucker-Punkt, die Bedingungen (4.15), (4.16), (4.17) Kuhn-Tucker-Bedingungen.

Definition 10 Die Richtung p ∈ T (x^{∗}) ⊂ ℝ^{n}

T (x^{∗}) := {p ≠ 0 | J_{2}(x^{∗})p = 0, μ_{i}J_{3i}(x^{∗})p = 0 i ∈ I(x^{∗}), J̃_{3}(x^{∗})p ≥ 0}

heißt zulässige Richtung der linearisierten Nebenbedingungen.

Schließlich geben wir noch die notwendigen Optimalitätsbedingungen zweiter Ordnung an.

Satz 7 (Notwendige Optimalitätsbedingungen zweiter Ordnung) Falls x^{∗} ∈ D ein regulärer Punkt und Lösung des restringierten nichtlinearen Quadratmittelproblems (4.13) ist, so gilt für alle zulässigen Richtungen der linearisierten Nebenbedingungen p ∈ T (x^{∗})

p^{T} H(x^{∗}, λ^{∗}, μ^{∗})p ≥ 0 mit H(x^{∗}, λ^{∗}, μ^{∗}) := \frac{∂^{2}}{∂x^{2}} L(x^{∗}, λ^{∗}, μ^{∗})  (4.18)

Satz 8 (Hinreichende Optimalitätsbedingungen zweiter Ordnung) Sei (x^{∗}, λ^{∗}, μ^{∗}) Kuhn-Tucker-Punkt von (4.13) und für alle zulässigen Richtungen der linearisierten Nebenbedingungen p ∈ T (x^{∗}) sei die Hesse-Matrix der Lagrange-Funktion positiv definit

p^{T} H(x^{∗}, λ^{∗}, μ^{∗})p > 0.  (4.19)

Dann ist x^{∗} striktes lokales Minimum von (4.13).

4.2.2 Verallgemeinertes Gauß-Newton-Verfahren

Zur numerischen Lösung des restringierten nichtlinearen Quadratmittelproblems (4.13) betrachten wir eine Verallgemeinerung des Gauß-Newton-Verfahrens (siehe [Boc87]). Grundlage ist die wiederholte Lösung eines linearisierten restringierten Quadratmittelproblems.

Linearisiertes restringiertes Quadratmittelproblem LP (x):

min \{\frac{1}{2} ∥F_{1}(x) + J_{1}(x)s∥^{2} : F_{2}(x) + J_{2}(x)s = 0, F_{3}(x) + J_{3}(x)s ≥ 0, s ∈ ℝ^{n}\}

(4.20)

Grundalgorithmus Verallgemeinertes Gauß-Newton-Verfahren

S1: Wähle Startpunkt x^{0} ∈ ℝ^{n}

S2: For k = 0, 1, … do

S2.1: Bestimme s^{k} als Lösung des linearisierten restringierten Quadratmittelproblems LP (x^{k})

∥F_{1}(x^{k}) + J_{1}(x^{k})s∥^{2} ⟶ min

unter den Nebenbedingungen

F_{2}(x^{k}) + J_{2}(x^{k})s = 0
F_{3}(x^{k}) + J_{3}(x^{k})s ≥ 0
S2.2: Setze x^{k+1} := x^{k} + s^{k};

Betrachten wir die hinreichenden Optimalitätsbedingungen (4.19) für das linearisierte restringierte Quadratmittelproblem (4.20). Die Hesse-Matrix der Lagrange-Funktion ist gerade J_{1}(x)^{T} J_{1}(x). Die Bedingung

p^{T} J_{1}^{T} J_{1}p > 0 für alle p ∈ Z := {p ≠ 0 : J_{c}p = 0}  (4.21)

ist äquivalent mit

Rang [ \begin{matrix}J_{1} \\ J_{c} \end{matrix} ] = n.  (4.22)

Sei (s^{∗}, λ^{∗}, μ^{∗}) Kuhn-Tucker-Punkt des linearisierten restringierten Quadratmittelproblems (4.20) an der Stelle x^{∗} und die hinreichende Bedingung (4.21) bzw. (4.22) sei erfüllt. Bock [Boc87] zeigt, daß es dann eine Umgebung X von x^{∗} gibt, so daß ∀x ∈ X das linearisierte restringierte Quadratmittelproblem LP (x) genau einen Kuhn-Tucker-Punkt (s, λ, μ) besitzt und das s striktes Minimum ist. Unter der Voraussetzung der strikten Komplementarität der Kuhn-Tucker-Vorzeichenbedingung (4.17)

i ∈ I(x^{∗}) ⟹ μ_{i}^{∗} > 0  (4.23)

gilt: I_{LP (x)}(s) ist konstant auf X, d.h. die Indexmenge der aktiven Ungleichheitsnebenbedingungen ändert sich nicht mehr.

Bemerkung 14 Insbesondere kann bei Vorliegen strikter Komplementarität kein Zigzagging, d.h. das Oszillieren zwischen verschiedenen Teilmengen der Nebenbedingungen, auftreten, welches die Konvergenz eines numerischen Algorithmus drastisch verschlechtern kann (siehe [Fle87, S.268]).

Damit ist gezeigt, daß das lokale Konvergenzverhalten des ungleichungrestringierten Problems dem des gleichungsrestringierten Problems entspricht, in dem die Gleichungs- und aktiven Ungleichungsnebenbedingungen zusammengefaßt werden. Wir beschränken uns deshalb im weiteren auf solche Probleme mit dem Grundgebiet D = X in der Umgebung der Lösung.

Voraussetzungen: Für alle x ∈ D gelte

Rang J_{c}(x) = m_{c}  (4.24)
Rang J(x) = n  J(x) := [ \begin{matrix}J_{1}(x \\ J_{c}(x) \end{matrix} ]  (4.25)

Bock [Boc87] zeigt nun, daß sich die Lösung dieses gleichungsrestringierten linearen Quadratmittelproblems durch eine verallgemeinerte Inverse ausdrücken läßt. Mit dem folgenden Satz aus [Boc87] ist damit eine einheitliche Konvergenztheorie für freie und restringierte Quadratmittelprobleme möglich.

Satz 9 (Existenz einer verallgemeinerten Inversen) Betrachte das gleichungsrestringierte lineare Quadratmittelproblem

∥J_{1}s + F_{1}∥^{2} ⟶ min  (4.26)

unter der Nebenbedingung

J_{c}s + F_{c} = 0.  (4.27)

J_{1} und J_{c} erfüllen die Voraussetzungen (4.24), (4.25). Dann gilt:

(i) Für beliebiges F := ( \begin{matrix}F_{1} \\ F_{c} \end{matrix} ) ∈ ℝ^{m_{1}+m_{c}} existiert genau ein Kuhn-Tucker-

Punkt (s, λ_{c}) von (4.26), (4.27), und s ist striktes Minimum.

(ii) Es existiert eine lineare Abbildung J^{−} : ℝ^{m_{1}+m_{c}} ⟶ ℝ^{n}, so daß s = −J^{−}F Lösung von (4.26), (4.27) ist für beliebiges F ∈ ℝ^{m_{1}+m_{c}}.

(iii) Der Lösungsoperator J^{−} ist eine verallgemeinerte Inverse, d.h.

(P 2) J^{−}JJ^{−} = J^{−}.

Bemerkung 15 • Unter den Voraussetzungen (4.24) und (4.25) gelten auch die Moore-Penrose Bedingungen (P1) und (P4).

• Unter den Voraussetzungen (4.24), (4.25) ist J(x)^{−} stetig differenzierbar auf X (vgl. Abschnitt 5).

Die lokale Konvergenz des verallgemeinerten Gauß-Newton-Verfahrens folgt nun aus Satz 9 und dem folgenden Satz.

Satz 10 (Kontraktionssatz, [Boc87]) Sei J(x)^{−} verallgemeinerte Inverse zur Ableitung J(x) von F ∈ C^{1}(D), und die Lipschitzbedingungen an J bzw. J^{−}

∥J(y)^{−} [J (x + t(y − x)) − J(x)] (y − x)∥ ≤ ωt ∥y − x∥^{2} mit ω < ∞ (4.28)
∥[J(z)^{−} − J(x)^{−}] R(x)∥ ≤ κ(x) ∥z − x)∥ ≤ κ ∥z − x∥ (4.29)
R(x) := F(x) − J(x)J(x)^{−}F(x) mit κ < 1

seien erfüllt ( ∀t ∈ [0, 1], ∀x, y, z ∈ D, x − y = J(x)^{−}F(x)). Für alle Startwerte x^{0} ∈ D mit

δ_{0} := \frac{α_{0}ω}{2} + κ < 1  α_{k} := ∥J(x^{k})^{−}F(x^{k})∥
D_{0} := K̄ (x^{0}, \frac{α_{0}}{1 − δ_{0}} ) ⊂ D

folgt dann:

(i) die Iteration x^{k+1} = x^{k} + Δx^{k}, Δx^{k} := −J(x^{k})^{−}F(x^{k}) ist wohldefiniert und bleibt in D_{0}, (ii) x^{k} ⟶ x^{∗} ∈ D_{0} mit J(x^{∗})^{−}F(x^{∗}) = 0, und es gilt die

(iii) a priori Abschätzung

∥x^{k} − x^{∗}∥ ≤ δ_{0}^{k} \frac{α_{0}}{1 − δ_{0}} .

(iv) Die Konvergenz ist linear

∥Δx^{k+1}∥ ≤ (\frac{α_{k}ω}{2} + κ) ∥Δx^{k}∥ .

Bemerkung 16 • Die Bedingung κ < 1 kennzeichnet im statistischen Sinne vernünftig gestellte Quadratmittelprobleme. Bock [Boc87] zeigt, daß unter der Voraussetzung κ < 1 der berechnete Fixpunkt x^{∗} nicht nur stationärer Punkt des Problems, sondern auch lokales Minimum und stabil gegen lokale Störungen ist. Umgekehrt wird gezeigt, daß andernfalls Störungen existieren, für die der Fixpunkt zwar stationär, aber kein Minimum ist.

• Alle Annahmen und Voraussetzungen des Kontraktionssatzes sind affin invariant.

4.2.3 Globalisierung

Zur Globalisierung des soeben beschriebenen Verallgemeinerten Gauß-Newton-Ver-fahrens führen wir eine Dämpfung aus:

x^{k+1} := x^{k} + α_{k}s^{k}.  (4.30)

Der Schrittweitenparameter α_{k} wird dabei so bestimmt, daß eine sogenannte Niveaufunktion einen ausreichenden Abstieg signalisiert. Neben dem Wert der Zielfunktion muß bei restringierten Quadratmittelproblemen auch berücksichtigt werden, wie stark die Nebenbedingungen verletzt werden. Bock [Boc87] betrachtet die Niveaufunktion

T (x) := \frac{1}{2} ∥F_{1}(x)∥^{2} + ∑_{i=1}^{m_{2}} β_{i} |F_{2i}(x)| + ∑_{i=1}^{m_{3}} γ_{i} |min (0, F_{3i}(x))| .

Die Bestimmung des Schrittweitenparameters α_{k} erfolgt über exakte Strahlsuche.

Gedämpftes Verallgemeinertes Gauß-Newton-Verfahren

S1: Wähle Startpunkt x^{0} ∈ ℝ^{n}; k := 0;

S2: Löse LP (x^{k});

S3: Bestimme α_{k}, so daß

T (x^{k} + α_{k}s^{k}) = min {T (x^{k} + αs) : α ∈ [0, 1]}

S3: Setze x^{k+1} := x^{k} + α_{k}s^{k}; k := k + 1; Goto S2;

Bemerkung 17 Die in Schritt S3 angenommene exakte Strahlsuche kann durch eine näherungsweise ersetzt werden (etwa Armijo-Goldstein-Schrittweitenstrategie).

Mit dieser Niveaufunktion beweist Bock globale Konvergenz des Verallgemeinerten Gauß-Newton-Verfahrens. Auf dem Grundgebiet D seien dabei die Regularitätsannahmen (4.24), (4.25) erfüllt.

Satz 11 (Globale Konvergenz des Verallgemeinerten GN-Verfahrens) (i) Sei x^{0} ∈ D beliebig, und die zugehörige Niveaumenge

D_{0} := {x | T (x) ≤ T (x^{0})} ⊂ D

sei kompakt.

(ii) Für alle x ∈ D_{0} gelte für den eindeutigen Kuhn-Tucker-Punkt (s, λ, μ) des linearisierten restringierten Problems LP (x)

|λ_{i}| < β_{i} |μ_{i}| < γ_{i}

Dann erzeugt das gedämpfte Verallgemeinerte Gauß-Newton-Verfahren eine Folge {α_{k}}, so daß die relaxierte Iteration (4.30) mit Startwert x^{0} einen Kuhn-Tucker-Punkt des nichtlinearen Problems (4.13) als Häufungspunkt besitzt.

4.2.4 Nichtlineare Quadratmittelprobleme mit linearen Ungleichheitsnebenbedingungen

In diesem Abschnitt schränken wir uns auf nichtlineare Quadratmittelprobleme mit linearen Ungleichheitsnebenbedingungen ein.

min \{\frac{1}{2} ∥F(x)∥^{2} : Cx ≥ h, x ∈ ℝ^{n}\}  (4.31)

mit F : ℝ^{n} ⟶ ℝ^{m}, C ∈ ℝ^{ncstr,n}, h ∈ ℝ^{ncstr}. Wright [Wri92] unterscheidet die aktuell verwendeten Algorithmen zur Lösung von (4.31) in zwei Klassen. Auf der einen Seite gibt es die active set Zugänge. Bei diesen Algorithmen wird in jeder äußeren Iteration eine Modellfunktion von g(x) = \frac{1}{2} ∥F(x)∥^{2} gebildet. In jeder inneren Iteration wird diese Modellfunktion dann über dem zulässigen Bereich {x ∈ ℝ^{n} |Cx ≥ h} minimiert. Bei der zweiten Klasse, den sogenannten gradient projection methods, wird in jeder Iteration eine Suchrichtung p gewählt. Anschließend wird auf dem Strahl P(x − αp) gesucht, wobei α > 0 und P die Projektion auf den zulässigen Bereich ist. Die von uns verwendete Strategie folgt dem im vorigen Abschnitt dargelegten Verallgemeinerten Gauß-Newton-Verfahren und läßt sich damit der ersten Klasse zuordnen. Einen allgemeinen Überblick über linear restringierte Optimierungsprobleme gibt [Fle87, Chapter 21]. Bei der Anwendung der Methoden aus Abschnitt 4.2.2 auf das Problem (4.31) ergeben sich einige Vereinfachungen. So ist die Regularitätsforderung (constraint qualification) automatisch erfüllt. Als Spezialisierung des linearisierten restringierten Problems LS auf (4.31) erhalten wir

Linearisiertes restringiertes Quadratmittelproblem LSI(x):

min \{\frac{1}{2} ∥F(x) + J(x)s∥^{2} : C(x + s) ≥ h, s ∈ ℝ^{n}\}  (4.32)

Dieses Problem kann unter der Voraussetzung Rang J(x) = n mit den in Abschnitt 3.2 entwickelten Methoden gelöst werden.

Bemerkung 18 Die Vollrangbedingung Rang J(x) = n ist analog zu der Vollrangbedingung für das Gauß-Newton-Verfahren für freie Quadratmittelprobleme. Sie sichert zugleich die

Rangbedingung (4.25) Rang [ \begin{matrix}J(x) \\ C_{a} \end{matrix} ] = n ( C_{a} aktive Nebenbedingungen) und

damit die Existenz eines strikten lokalen Minimums s.

Die Globalisierung erfolgt wieder über Dämpfung. Es ergeben sich jedoch auch hier einige wichtige Vereinfachungen. Die Terme

∑_{i=1}^{m_{2}} β_{i} |F_{2i}(x)| + ∑_{i=1}^{m_{3}} γ_{i} |min (0, F_{3i}(x))|

in der Definition der Niveaufunktion entfallen, da das Verfahren nur mit zulässigen Punkten arbeitet. Wir arbeiten stattdessen wieder mit der Armijo-Goldstein-Schrittweitenstrategie für die Funktion g(x) = \frac{1}{2} ∥F(x)∥^{2}. Abschließend formulieren wir das verwendete Verallgemeinerte Gauß-Newton-Ver-fahren zur Lösung von (4.31).

Gedämpftes Verallgemeinertes Gauß-Newton-Verfahren

S1: Wähle (zulässigen) Startwert x^{0} ∈ ℝ^{n}, α ∈ (0, 1), δ ∈ (0, \frac{1}{2})

Setze k := 0; γ_{−1} := 1;

S2: Iteration Repeat

S2.1: F := F(x^{k}); Berechnung des Residiums

J := F′(x^{k}); Berechnung der Jacobi-Matrix S2.2: Bestimme Abstiegsrichtung s^{k}(= −p^{k}) aus

∥F + Js∥^{2} ⟶ min
mit  C (x^{k} + s) ≥ h

S2.3: γ̃_{k} := γ_{k}; x^{+} := x^{k} + γ̃_{k}s^{k}; F^{+} := F(x^{+}); S2.4: Abstiegstest

if \frac{1}{2} ∥F∥^{2} − \frac{1}{2} ∥F^{+}∥^{2} ≥ δγ̃_{k}F^{T} Jp^{k}

then

x^{k+1} := x^{+};
γ_{k+1} := min(γ_{k}/α, 1);
k := k + 1;

else

γ_{k} := γ_{k} ∗ α; goto S2.3;

until (k > kmax) or Konvergenz

Bemerkung 19 • Die Berechnung der Abstiegsrichtung s^{k} in Schritt S2.2 erfolgt mittels des Algorithmus LSI aus Abschnit 3.2.

• Die Übertragung der Armijo-Goldstein-Schrittweitenstrategie auf linear restringierte Probleme bereitet keinerlei Schwierigkeiten, da alle Punkte längs des Strahls x^{k} − αp^{k}, α ∈ [0, 1] zulässig sind. Bei nichtlinearen Restriktionen ist die Verfolgung eines zulässigen Kurvenbogens x(α) mittels eines Korrekturterms notwendig, was die ”Strahl“suche i.a. sehr verteuert. Eine Armijo-Goldstein-Schrittweitenstrategie für solche nichtlinear restringierten Probleme beschreibt Spellucci [Spe93].

• Die Bemerkungen zum gedämpften Gauß-Newton-Verfahren mit akzeptiertem Schätzwert gelten sinngemäß.

4.2.5 Regularisierung

Die Durchführbarkeit und Konvergenz von (Verallgemeinerten) Gauß-Newton-Ver-fahren benötigt als wesentliche Voraussetzung die Vollrangeigenschaft Rang J(x) = n der Jacobi-Matrix. Wenn die Matrix J (fast) rangdefizient ist, so ist der Algorithmus LSI zur Lösung des linearisierten restringierten Teilproblems nicht durchführbar und bricht ab. Im Fall der freien Quadratmittelprobleme erfolgt eine Regularisierung meist über ein Trust-Region-Verfahren, d.h.

min \{μ_{GN}(x + s) = g(x) + ∇g(x)^{T} s + \frac{1}{2}s^{T} H(x)s : s ∈ ℝ^{n} bei ∥s∥^{2} ≤ Δ^{2}\}

(4.33)

mit dem Trust-Region-Radius Δ = Δ_{k} > 0. Bis auf Ausnahmefälle ergibt sich die Lösung von (4.33) als Lösung von

(H + λI) s = −∇g,

wobei der Lagrange-Parameter λ ≥ 0 so gewählt werden muß, daß H + λI positiv definit ist und außerdem ∥s(λ)∥ = Δ gilt. Für das Gauß-Newton-Modell erhalten wir

min \{∥[ \begin{matrix}J \\ \sqrt{λ}I \end{matrix} ] s + ( \begin{matrix}F \\ 0 \end{matrix} )∥^{2} = ∥F + Js∥^{2} + λ ∥s∥^{2} : s ∈ ℝ^{n}\}  (4.34)

Steuert man den Trust-Region-Radius über den Lagrange-Parameter λ spricht man vom Levenberg-Marquardt-Verfahren. Für Trust-Region-Verfahren kann man unter schwachen Voraussetzungen globale Konvergenz beweisen, d.h. der Trust-Region-Ansatz dient nicht nur zur Regularisierung, sondern auch zur Globalisierung. Dennis/Schnabel [DS83, S.151] beschreiben eine Strategie zur Regularisierung des Newton-Verfahrens bei Verwendung eines gedämpften Verfahrens. Falls die Hesse-Matrix H = ∇^{2}g(x) des quadratischen Modells schlecht konditioniert ist, wird das quadratische Modell

min \{μ_{N}(x + s) = g(x) + ∇g(x)^{T} s + \frac{1}{2}s^{T} H(x)s : s ∈ ℝ^{n}\}  (4.35)

gestört:

H_{reg} := H + \sqrt{n × macheps} ∥H∥_{1} × I

Diese Strategie ist in dem in [DS83] enthaltenen Programmpaket realisiert. Dort findet sich auch die theoretische Fundierung der Wahl der Größe des Lagrange-Parameters. Wir verwenden diese Strategie zur Regularisierung des Verallgemeinerten Gauß-Newton-Verfahrens. Der Algorithmus LSI führt eine QR-Faktorisierung von J durch und schätzt die Kondition der oberen Dreicksmatrix R. Falls die Konditionszahl größer als \frac{1}{\sqrt{macheps}} ist, bricht der Algorithmus ab und teilt dies dem aufrufenden

Algorithmus mit. Wir berechnen dann

λ := \sqrt{n × macheps} ∥J^{T} J∥_{1}

und lösen das Problem

min \{∥[ \begin{matrix}J \\ \sqrt{λ}I \end{matrix} ] s + ( \begin{matrix}F \\ 0 \end{matrix} )∥^{2} : C (x^{k} + s) ≥ h, s ∈ ℝ^{n}\} .  (4.36)

Bemerkung 20 • Die soeben beschriebene Regularisierung wurde in den verwendeten Algorithmen implementiert und liefert gute numerische Ergebnisse.

• Die Verwendung eines Trust-Region-Verfahrens zur Lösung von nichtlinearen restringierten Quadratmittelproblemen beschreiben [HF79]. Die Steuerung des Trust-Region-Radius erfolgt dabei nicht über den Lagrange-Parameter λ, sondern über Δ, da Δ nicht so empfindlich bei der Änderung der aktiven Menge ist. Trust-Region-Verfahren für Probleme mit linearen Ungleichheitsnebenbedingungen findet man in [Gay84] und [WH85].

5 Separable Quadratmittelprobleme

5.1 Problemformulierung

Ein nichtlineares Quadratmittelproblem, bei denen ein Teil der Unbekannten linear auftritt, wird als (linear) separabel bezeichnet.

min \{ψ(t̃, α) := \frac{1}{2} ∥F (t̃, α)∥^{2} = \frac{1}{2}∥y − B(t̃)α∥^{2} : t̃ ∈ ℝ^{l}, α ∈ ℝ^{n}\}  (5.1)

y ∈ ℝ^{m} sei dabei ein vorgegebener Vektor und B(.) : ℝ^{l} ⟶ ℝ^{m×n} eine hinreichend glatte Matrixfunktion. Es gelte m ≥ n + l, d.h. das System (5.1) sei überbestimmt. Separable Quadratmittelprobleme sind ein Spezialfall sog. reduzibler nichtlinearer Optimierungsprobleme, bei denen eine natürliche Unterscheidung der Variablen in zwei Gruppen auftritt (siehe [Par85] und [RW80]).

min {f(t̃, α) : t̃ ∈ ℝ^{l}, α ∈ ℝ^{n}}  (5.2)

Die Variablen, welche den Vektor α darstellen, sind dabei so gewählt, daß das reduzierte (n-dimensionale) Problem

min {f(t̃_{c}, α) : t̃_{c} ∈ ℝ^{l} fest, α ∈ ℝ^{n}}  (5.3)

effektiv und akkurat gelöst werden kann. Typische Beispiele für separable Quadratmittelprobleme sind:

• Parameterschätzung bei Spektren

\frac{1}{2} ∑_{i=1}^{m} [y_{i} − r (x_{i}, t̃, α)]^{2} ⟶ Min!

r (x, t̃, α) = ∑_{j=1}^{k} α_{j} exp (− (\frac{x − t̃_{j}}{t̃_{j+k}} )^{2}) + α_{k+1} + α_{k+2}x + α_{k+3}x^{2}

mit den k + 3 linearen Parametern (α_{1}, …, α_{k+3})^{T} und den 2k nichtlinearen Parametern (t̃_{1}, …, t̃_{2k})^{T}

• B-Splines mit freien Knoten

\frac{1}{2} ∑_{i=1}^{m} [y_{i} − r (x_{i}, t̃, α)]^{2} ⟶ Min!

r (x, t̃, α) = ∑_{j=1}^{n} B_{j,k}(x, t̃)α_{j}

wobei B_{j,k}(x, t̃) der j-te normalisierte B-Spline der Ordnung k an der Stelle x zur Knotenfolge t̃ ist.

5.2 Gauß-Newton-ähnliche Verfahren für separable Quadratmittelprobleme

Wir benutzen das Gauß-Newton-Verfahren, um die Ideen zur effektiven Lösung von (5.1) darzulegen. Der Gauß-Newton-Schritt s ∈ ℝ^{n+l} wird aus der Gleichung

J^{T} Js = −J^{T} F  (5.4)

berechnet. Die Jacobi-Matrix J ∈ ℝ^{m,n+l} der Residiuenfunktion F ist

J = [J_{α} Jt̃]
mit J_{α} = B \binom{}{t̃}
und Jt̃ = −DB \binom{}{t̃} α.

DB bezeichne die Fréchet-Ableitung der Matrixfunktion B(.).

Eine m × n Matrix-Funktion B \binom{}{t̃} ist eine nichtlineare Abbildung zwischen dem linearen Raum der Parameter ℝ^{l} und dem Raum der linearen Abbildungen 𝔏 (ℝ^{n}, ℝ^{m}). DB \binom{}{t̃} ist also ein Element von 𝔏 (ℝ^{l}, 𝔏 (ℝ^{n}, ℝ^{m})) und kann somit als dreidimensionaler Tensor interpretiert werden.

Wir können nun den Gauß-Newton-Schritt s ∈ ℝ^{n+l} in einen Schritt s_{α} bezüglich der Variablen α und einen Schritt st̃ bezüglich t̃ aufspalten. Schreiben wir die Normalgleichungen (5.4) aus, so erhalten wir folgende Gleichungen für s_{α} und st̃.

[ \begin{matrix}J_{α}^{T}J_{α} & J_{α}^{T}Jt̃ \\ JT ; t̃ J_{α} & JT ; t̃ Jt̃ \end{matrix} ] ( \begin{matrix}s_{α} \\ st̃ \end{matrix} ) = ( \begin{matrix}−J_{α}^{T}F \\ −JT ; t̃ F \end{matrix} )

B^{T} Bs_{α} + B^{T} Jt̃st̃ = −B^{T} F  (5.5)
JT ; t̃ Bs_{α} + JT ; t̃ Jt̃st̃ = −JT ; t̃ F  (5.6)

Indem wir die Gleichung (5.5) von links mit JT ; t̃ (B^{+})^{T} multiplizieren und das Resultat von (5.6) subtrahieren, erhalten wir wegen der Beziehung B^{+}^{T} B^{T} B =

(BB^{+})^{T} B =^{(P 3)} BB^{+}B =^{(P 1)} B die Gleichung

JT ; t̃ (I − BB^{+}) Jt̃st̃ = −JT ; t̃ (I − BB^{+}) F.

Definieren wir den orthogonalen Projektor P_{B}^{⊥} := I − BB^{+}, so ergibt sich

JT ; t̃ P_{B}^{⊥}Jt̃st̃ = −JT ; t̃ P_{B}^{⊥}F.

Wegen der Symmetrie und Idempotenz eines Projektors ist dazu äquivalent

(P_{B}^{⊥}Jt̃)^{T} (P_{B}^{⊥}Jt̃) st̃ = − (P_{B}^{⊥}Jt̃)^{T} P_{B}^{⊥}F.

Dies sind aber genau die Normalgleichungen zu

P_{B}^{⊥}Jt̃st̃ ≅ −P_{B}^{⊥}F.

Die eindeutige Normallösung ist

st̃ = − (P_{B}^{⊥}Jt̃)^{+} P_{B}^{⊥}F.

Aus der Gleichung (5.5) erhalten wir dann

s_{α} = −B^{+} (F + Jt̃st̃) .

Jetzt sind wir in der Lage, das Gauß-Newton-Verfahren für das unseparierte Problem (5.1) zu formulieren:

Algorithmus G: Gauß-Newton-Verfahren für das unseparierte Problem

S1: Wähle Startpunkt t̃^{(0)} ∈ ℝ^{l}, α^{(0)} ∈ ℝ^{n}

S2: For k = 0, 1, … do

S2.1: Konvergenztest

S2.2: Berechne

s(k) ; t̃ := − (P_{B}^{⊥}Jt̃)^{+} P_{B}^{⊥}F
s_{α}^{(k)} := −B^{+} (F + Jt̃s(k) ; t̃)

S2.3: Setze

t̃^{(k+1)} := t̃^{(k)} + s(k) ; t̃
α^{(k+1)} := α^{(k)} + s_{α}^{(k)}

Das separable Quadratmittelproblem (5.1) besitzt offensichtlich mehr Struktur als ein beliebiges nichtlineares Quadratmittelproblem. Das wesentliche Merkmal separabler Quadratmittelprobleme ist die Tatsache, daß für feste t̃_{c} ∈ ℝ^{l} die Lösung des reduzierten Problems

min \{\frac{1}{2}∥y − B(t̃_{c})α∥^{2} : t̃_{c} ∈ ℝ^{l} fest, α ∈ ℝ^{n}\}  (5.7)

sehr effizient berechnet werden kann. Das reduzierte Problem (5.7) ist ein lineares Quadratmittelproblem der Dimension n, welches mit Standardorthogonalisierungstechniken gelöst werden kann. Eine naheliegende Veränderung des Gauß-Newton-Verfahrens besteht darin, abwechselnd über einer Variablengruppe zu minimieren, während die Variablen der anderen Gruppe konstant gehalten werden. Diese Methode ist in der Literatur als Methode der alternierenden Variablen bekannt.

Algorithmus 1: Methode der alternierenden Variablen

S1: Wähle Startpunkt t̃^{(0)} ∈ ℝ^{l}

S2: For k = 0, 1, … do

S2.1: Konvergenztest

S2.2: Löse das lineare Quadratmittelproblem

min \{\frac{1}{2}∥y − B \binom{}{t̃^{(k)}} α∥^{2} : α ∈ ℝ^{n}\} ,

um α^{(k)} zu erhalten

S2.3: Berechne

s(k) ; t̃ := −J+ ; t̃F

S2.4: Setze

t̃^{(k+1)} := t̃^{(k)} + s(k) ; t̃

Bemerkung 21 Algorithmus 1 ist der Algorithmus NIPALS, welcher von Wold und Lyttkens [WL69] stammt und unter Statistikern verbreitet ist.

Betrachten wir nun das lineare Quadratmittelproblem in Schritt S2.2 in Algorithmus 1. Dieses lineare Quadratmittelproblem hat im allgemeinen, rankdefizienten Fall viele Lösungen. Wir wählen die eindeutige Normallösung

α_{opt} = α_{opt}(t̃) = B(t̃)^{+}y  (5.8)

Falls die Matrix B Vollrang l hat, so erhalten wir

B(t̃)^{+} = [B(t̃)^{T} B(t̃)]^{−1} B(t̃)^{T}

Für jedes t̃ ∈ ℝ^{l} erhalten wir den optimalen Parameter α_{opt}(t̃) durch (5.8). Setzen wir diesen optimalen Parameter in das Ausgangsproblem

min \{ψ(t̃, α) := \frac{1}{2} ∥F (t̃, α)∥^{2} = \frac{1}{2}∥y − B(t̃)α∥^{2} : t̃ ∈ ℝ^{l}, α ∈ ℝ^{n}\} (5.9)

ein, so erhalten wir folgendes reduzierte Problem:

min\{ϕ(t̃) := \frac{1}{2}∥\underbrace{[I − B(t̃)B(t̃)^{+}] y}_{=:G(t̃)}∥^{2} = \frac{1}{2}∥G(t̃)∥^{2} : t̃ ∈ ℝ^{l}\}  (5.10)

Da P_{B}^{⊥} := I − B(t̃)B(t̃)^{+} der orthogonale Projektor auf im B(t̃)^{⊥} ist, wird das Funktional ϕ(t̃) auch als variable projection functional und die beschriebene Reduktionstechnik als Verfahren der variablen Projektion bezeichnet (siehe [GP73]). Das reduzierte Problem (5.10) enthält nur noch die l nichtlinear auftretenden Parameter t̃, ist dafür aber von einer größeren Komplexität als das Ausgangsproblem. Für die Anwendung des Verfahrens der variablen Projektion bleibt die Zulässigkeit der Transformation zu zeigen. Außerdem muß für die Anwendung von Gauß-Newton-ähnlichen Verfahren die Jacobi-Matrix G′ von G(t̃) berechnet werden.

5.3 Der Übergang zum reduzierten Funktional

Das folgende Theorem von Golub/Pereyra [GP73] rechtfertigt den Übergang vom Ausgangsfunktional ψ(t̃, α) zum reduzierten Funktional ϕ(t̃).

Theorem 2 (Golub/Pereyra [GP73]) Seien ψ(t̃, α) und ϕ(t̃) wie oben definiert. B(t̃) habe in der offenen Menge Ω ⊂ ℝ^{l} konstanten Rang r ≤ min(m, n).

(a) Wenn t̃^{∗} ein kritischer Punkt (oder eine globale Minimumstelle in Ω) von ϕ(t̃) ist und

α^{∗} = B(t̃^{∗})^{+}y,  (5.11)

so ist (t̃^{∗}, α^{∗}) ein kritischer Punkt (oder eine globale Minimumstelle für t̃ ∈ Ω) von ψ(t̃, α), und es gilt

ϕ(t̃^{∗}) = ψ(t̃^{∗}, α^{∗}).

(b) Wenn (t̃^{∗}, α^{∗}) eine globale Minimumstelle von ψ(t̃, α) für t̃ ∈ Ω ist, so ist t̃^{∗} eine globale Minimumstelle von ϕ(t̃) in Ω, und es gilt

ϕ(t̃^{∗}) = ψ(t̃^{∗}, α^{∗}).

Wenn es ein eindeutiges α^{∗} unter den Minimumstellen von ψ(t̃, α) gibt, so muß α^{∗} die Beziehung (5.11) erfüllen.

Bemerkung 22 Da die Moore-Penrose-Inverse B^{+} bei rangerhöhenden Störungen nicht stetig ist (siehe z.B. [KS88]), ist die Forderung nach einem lokal konstanten Rang der Matrix B natürlich. Wir machen daher im weiteren die Voraussetzung

Die Matrix B(t̃) habe lokal konstanten Rang r ≤ min(m, n) für alle t̃ ∈ Ω ⊂ ℝ^{l}, wobei Ω eine offene Menge ist, welche die gesuchte Lösung enthält.

5.4 Die Bestimmung der Jacobi-Matrix

Zur Berechnung der Jacobi-Matrix G′ benötigen wir zunächst die Fréchet-Ableitung eines orthogonalen Projektors P_{A(t̃)} einer differenzierbaren m × n Matrix-Funktion A(t̃) von lokal konstantem Rang r. Wir werden dabei zunächst nicht die Moore-Penrose-Inverse A^{+} und den zugehörigen Projektor P_{A} = AA^{+} betrachten, sondern eine verallgemeinerte Inverse A^{−}, welche nicht alle vier Moore-Penrose-Bedin-gungen erfüllt. Der zugehörige Projektor P_{A} = AA^{−} wird dabei nicht geändert.

Lemma 2 (Golub/Pereyra [GP73]) Sei A^{−}(t̃) eine n×m Matrixfunktion, so daß (P1) AA^{−}A = A und (P3) (AA^{−})^{T} = AA^{−}. Dann gilt:

DP_{A} = P_{A}^{⊥} (DA) A^{−} + (P_{A}^{⊥} (DA) A^{−})^{T}

Beweis: Wegen P_{A}A = A gilt

D (P_{A}A) = (DP_{A}) A + P_{A} (DA) = DA,

d.h.

(DP_{A}) A = DA − P_{A} (DA) = P_{A}^{⊥} (DA)

und daher wegen P_{A} = AA^{−} schließlich

(i)  (DP_{A}) P_{A} = (DP_{A}) AA^{−} = P_{A}^{⊥} (DA) A^{−}.

Wegen (ii) ((DP_{A}) P_{A})^{T} = P_{A} (DP_{A}) erhalten wir aus (i) und (ii) letztendlich

DP_{A} = D (P_{A}^{2}) = (DP_{A}) P_{A} + P_{A} (DP_{A})
= P_{A}^{⊥} (DA) A^{−} + (P_{A}^{⊥} (DA) A^{−})^{T}

□

Bemerkung 23 • Das Lemma gilt natürlich auch, wenn die verallgemeinerte Inverse A^{−} durch die Moore-Penrose-Inverse A^{+} ersetzt wird.

• Wegen P_{A}^{⊥} = I − P_{A} gilt DP_{A}^{⊥} = −DP_{A}.

Wir sind nun in der Lage, den Gradienten des reduzierten Funktionals

ϕ(t̃) = \frac{1}{2} ∥P_{B(t̃)}^{⊥}y∥^{2}
= \frac{1}{2} ⟨P_{B}^{⊥}y, P_{B}^{⊥}y⟩

zu berechnen. Es gilt:

∇ϕ(t̃) = ⟨P_{B}^{⊥}y, D (P_{B}^{⊥}y)⟩
= − ⟨P_{B}^{⊥}y, (DP_{B}) y⟩
= − ⟨P_{B}^{⊥}y, [P_{B}^{⊥}(DB)B^{−} + (B^{−})^{T} (DB)^{T} P_{B}^{⊥}] y⟩
= −y^{T} P_{B}^{⊥} [P_{B}^{⊥}(DB)B^{−} + (B^{−})^{T} (DB)^{T} P_{B}^{⊥}] y

Beachte: Diese Darstellung für den Gradienten des reduzierten Funktionals gilt für beliebige orthogonale Projektoren P_{B}^{⊥} = I − BB^{−}, sofern B^{−} die Moore-Penrose-Bedingungen (P1) und (P3) erfüllt.

Wir nehmen jetzt an, daß B^{−} zusätzlich die Bedingung (P2) B^{−}BB^{−} = B^{−} erfüllt. Dann gilt:

P_{B}^{⊥} (B^{−})^{T} = [I − BB^{−}] (B^{−})^{T}
= (B^{−})^{T} − BB^{−} (B^{−})^{T}
= (B^{−})^{T} − (BB^{−})^{T} (B^{−})^{T}
= (B^{−})^{T} − (B^{−}BB^{−})^{T} = 0

Falls B^{−} die Moore-Penrose-Bedingungen (P1), (P2) und (P3) erfüllt, erhalten wir

für den Gradienten des reduzierten Funktionals ϕ(t̃) = \frac{1}{2} ∥P_{B(t̃)}^{⊥}y∥^{2} den Ausdruck

∇ϕ(t̃) = −y^{T} P_{B(t̃)}^{⊥} (DB(t̃)) B(t̃)^{−}y.  (5.12)

Kommen wir nun zur angekündigten Berechnung der Jacobi-Matrix G′ der reduzierten Funktion G.

G(t̃) = [I − B(t̃)B(t̃)^{+}] y = P_{B(t̃)}^{⊥}y
G′(t̃) = D (P_{B}^{⊥}y) = − (DP_{B}y)
= −P_{B}^{⊥}(DB)B^{+}y − (P_{B}^{⊥}(DB)B^{+})^{T} y

Wir erhalten daher für beliebige s ∈ ℝ^{l}

G′(t̃)s = {A(t̃)[s] + (A(t̃)[s])^{T} } y  (5.13)
mit A(t̃)[s] := −P_{B}^{⊥} (B′(t̃)[s]) B(t̃)^{+}  (5.14)
= − [I − B(t̃)B(t̃)^{+}] (B′(t̃)[s]) B(t̃)^{+}  (5.15)

Nun können wir das Gauß-Newton-Modell für das reduzierte Problem (5.10) bilden:

μ_{GP} (t̃ + s) := \frac{1}{2} ∥G + G′s∥^{2}  (5.16)

Bemerkung 24 • Im weiteren werden wir dieses Modell als Golub/Pereyra-Modell bezeichnen.

• In der Originalarbeit von Golub/Pereyra [GP73] wird anstatt mit der Pseudoinversen B^{+} mit einer symmetrischen verallgemeinerten Inversen B^{−}, welche nur die Bedingungen (P1), (P2) und (P3) erfüllt, gearbeitet. Diese verallgemeinerte Inverse wird über eine QR-Transformation auf Trapezform berechnet (siehe etwa [GvL89]).

Jetzt sind wir in der Lage, das Gauß-Newton-Verfahren für das reduzierte Problem (5.10) zu formulieren.

Algorithmus 2: Gauß-Newton-Verfahren für das reduzierte Problem

S1: Wähle Startpunkt t̃^{(0)} ∈ ℝ^{l}

S2: For k = 0, 1, … do

S2.1: Konvergenztest, Falls Konvergenz goto S3; S2.2: Berechne

s(k) ; t̃ := − {− [P_{B}^{⊥}(DB)B^{+} + (P_{B}^{⊥}(DB)B^{+})^{T} ] y}^{+} P_{B}^{⊥}y

S2.3: Setze

t̃^{(k+1)} := t̃^{(k)} + s(k) ; t̃

S3: Setze α^{∗} := B(t̃)y

Bemerkung 25 Eine effiziente Implementierung dieses Verfahrens findet sich in Krogh [Kro74]. Dort wird angenommen, daß B Vollrang hat.

Sowohl in den Implementierungen von Golub/Pereyra als auch in der von Krogh muß der Tensor DB gebildet werden, da er einmal als DB als auch als (DB)^{T} auftritt. Wie oben (5.12) gezeigt wurde, trägt der Term

(A(t̃))^{T} y := − (P_{B}^{⊥}(DB)B^{+})^{T} y

jedoch gar nicht zum Gradienten ∇ϕ(t̃) bei. Eine genaue Analyse zeigt, daß für das Golub/Pereyra-Modell (5.16) gilt

∥G + G′s∥^{2} = ∥G + A[s]y∥^{2} + s^{T} T_{G}s,  (5.17)

wobei die Matrix T_{G} gemäß

u^{T} T_{G}v := {(A[u])^{T} y}^{T} {(A[u])^{T} y}  (5.18)
:= G^{T} (B′[u]) B^{+}B^{+}^{T} (B′[v])^{T} G  (5.19)

definiert ist. Der Beitrag des Terms (A(t̃)[s])^{T} y zur Hesse-Matrix ∇^{2}ϕ(t̃) des re-

duzierten Funktionals ist T_{G}, dieser ist aber von der Größenordnung O (∥G∥^{2}). Er ist somit für kleines ∥G∥ kleiner als der beim Gauß-Newton-Verfahren ohnehin weggelassene Term S = G ∘ G′′. Entsprechend der Gauß-Newton-Philosophie wird der Term T_{G} daher ebenfalls vernachlässigt. Dies führt auf das vereinfachte Modell

μ_{K}(t̃ + s) := \frac{1}{2}∥G + Js∥^{2}  (5.20)
mit Js := A(t̃)[s]y = − [I − B(t̃)B(t̃)^{+}] (B′(t̃)[s]) B^{+}y, (5.21)

welches nach Kaufman [Kau75] benannt ist.

Bemerkung 26 Die soeben beschriebene Vereinfachung wurde erstmals von Kaufman [Kau75] vorgenommen. Sie wird dort unter Benutzung einer QR-Transformation auf Trapezform hergeleitet. Eine direkte Herleitung, an welche sich die obige Darstellung anlehnt, findet sich in [RW80]. Der Kaufman-Zugang eignet sich nach [Bjö 90] besser zur Verallgemeinerung auf restringierte Probleme (siehe auch [KP78]).

Algorithmus 3: Gauß-Newton-Verfahren für das red. Problem (Kaufman-Modell)

S1: Wähle Startpunkt t̃^{(0)} ∈ ℝ^{l}

S2: For k = 0, 1, … do

S2.1: Konvergenztest, Falls Konvergenz goto S3;

S2.2: Berechne

s(k) ; t̃ := − {− [P_{B}^{⊥}(DB)B^{+}] y}^{+} P_{B}^{⊥}y

S2.3: Setze

t̃^{(k+1)} := t̃^{(k)} + s(k) ; t̃

S3: Setze α^{∗} := B(t̃)y

5.5 Konvergenzraten und Aufwand

Ruhe/Wedin [RW80] haben die asymptotischen Konvergenzraten der vorgestellten Algorithmen untersucht und geben hinreichende Bedingungen für die lokale Konvergenz der ungedämpften Verfahren an. Sei

ϵ := \begin{matrix}∥y − B \binom{}{t̃} α∥ & 2 \\ ∥B \binom{}{t̃} α∥ & 2 \end{matrix} = ∥[I − B \binom{}{t̃} B \binom{}{t̃}^{+}] y∥_{2} ; ∥B \binom{}{t̃} B \binom{}{t̃}^{+} y∥_{2}

ein Maß für die relative Größe des Residiums im Lösungspunkt. Das Residium wird als klein betrachtet, wenn ϵ klein ist. Für jede Fixpunktiteration

z_{k+1} := h(z_{k}), lim_{k→∞} z_{k} = z^{∗},

ist die R-Konvergenzrate als

R = ρ (h′(z^{∗}))

definiert (ρ – Spektralradius).

Satz 12 (asymptotische Konvergenzraten) Die R-Konvergenzraten der Algorithmen hängen in folgender Weise von ϵ ab:

R_{G} = R_{II} + O(ϵ^{2}) Gauß-Newton unsepariertes Problem
R_{I} = O(1)  Alternierende Variablen
R_{II} = O(ϵ)  Golub/Pereyra-Modell
R_{III} = R_{II} + O(ϵ^{2}) Kaufman-Modell

Jeder der vier Algorithmen ist also im allgemeinen R-linear konvergent. Wenn der Gauß-Newton-Algorithmus für das unseparierte Problem superlinear konvergiert, so konvergieren auch die Algorithmen II und III superlinear. Die Methode der alternierenden Variablen dagegen konvergiert stets nur R-linear. In Kaufman [Kau75] und Ruhe/Wedin [RW80] werden die Algorithmen für freie Optimierungsprobleme verglichen und festgestellt, daß der Kaufman-Algorithmus mit weniger arithmetischen Operationen nicht mehr Iterationen als der ursprüngliche Algorithmus von Golub/Pereyra benötigt. Sie konvergieren beide schneller als der Gauß-Newton-Algorithmus für das unseparierte Problem. Dabei ist jedoch zu beachten, daß stets nur bestimmte Implementierungen der Algorithmen verglichen werden, und daß das Abschneiden der Algorithmen stark vom Aufwand zur Berechnung der Residiumfunktion abhängig ist. Ein Vorteil der Algorithmen für separable Quadratmittelprobleme liegt in der Tatsache, daß keine Anfangswerte für die linearen Parameter benötigt werden. Eine eventuelle spezielle Struktur (Schwachbesetztheit) von B′ (dies ist insbesondere bei der Approximation mit B-Splines der Fall) kann in den Algorithmen von Golub-Pereyra und Kaufman gut ausgenutzt werden.

6 Splineapproximation mit freien Knoten

6.1 Problemstellung

Gegeben seien die Meßdaten (x_{i}, y_{i}), i = 1, …, m, mit

a ≤ x_{1} < … < x_{m} ≤ b

einer unbekannten glatten Funktion f. Die mit Fehlern behafteten Meßwerte

y_{i} = f(x_{i}) + ϵ_{i},  i = 1, …, m
ϵ = (ϵ_{1}, …, ϵ_{m})^{T} = N (0, σ^{2}I)

sollen durch eine einfach zu berechnende Funktion s approximiert werden. Wir wählen als Ansatzfunktionen die B-Splines B_{1,k}, …, B_{n,k} der Ordnung k bezüglich einer Knotenfolge t = (t_{1}, …, t_{n+k}). Dabei gelte m ≥ n. Für die Knotenfolge t gelte

t_{1} = … = t_{k} = a < t_{k+1} ≤ … ≤ t_{n} < b = t_{n+1} = … = t_{n+k}  (6.1)
und  t_{j+k} > t_{j}  j = 1, …, n.  (6.2)

Die Ansatzfunktion s lautet dann

s(x) := ∑_{j=1}^{n} B_{j,k}(x, t)α_{j} , α_{j} ∈ ℝ.

Die Bezeichnung B_{j,k}(x, t) deutet dabei an, daß die B-Splines von der Knotenfolge t abhängig sind. Unter der Voraussetzung (6.2) hat der Splineraum

S_{k,t} := {∑_{j=1}^{n} B_{j,k}(x, t)α_{j} , α_{j} ∈ ℝ}

die Dimension n. Für die ”wahren“ Meßwerte ỹ_{i} erhalten wir den Ansatz

ỹ_{i} = ∑_{j=1}^{n} B_{j,k}(x_{i}, t)α_{j}

und mittels der diskreten l_{2}-Norm das Minimierungsproblem

∑_{i=1}^{m} [y_{i} − ỹ_{i}]^{2} ⟶ min

∑_{i=1}^{m} [y_{i} − ∑_{j=1}^{n} B_{j,k}(x_{i}, t)α_{j}]^{2} ⟶ min .  (6.3)

Wir definieren die Beobachtungsmatrix

B(t) = [ \begin{matrix}B_{1,k}(x_{1}, t) & … & B_{n,k}(x_{1}, t) \\ ⋮ & & ⋮ \\ B_{1,k}(x_{m}, t) & … & B_{n,k}(x_{m}, t) \end{matrix} ] ∈ ℝ^{m,n}

und erhalten mit den Splinekoeffizienten α = (α_{1}, …, α_{n})^{T} ∈ ℝ^{n} die zu (6.3) äquivalente Matrixformulierung

∥y − B(t)α∥^{2} ⟶ min .  (6.4)

Die Matrix B(t) hat unter der Bedingung (6.2) genau dann Vollrang n, wenn die Whitney-Schoenberg-Bedingung (2.6)

t_{j} < x_{i_{j}} < t_{j+k}

erfüllt ist. Die Splineknoten t und die Splinekoeffizienten α sollen jetzt so bestimmt werden, daß ∥y − B(t)α∥^{2} möglichst klein wird. Um größere Flexibilität zu erreichen, wählen wir nur eine Teilmenge t̃, die sogenannten freien Knoten, der Splineknoten t aus. l sei die Anzahl der freien Knoten und

t̃ = (t_{p(1)}, …, t_{p(l)})^{T} ∈ ℝ^{l}

sei der Vektor der freien Knoten. p = (p(1), …, p(l))^{T} ist ein Indexvektor, welcher auf die freien Knoten zeigt. Es gelte

k < p(1) < p(2) < … p(l) < n + 1.  (6.5)

Bemerkung 27 • Die Bedingung (6.5) besagt, daß nur die inneren Knoten t_{k+1}, …, t_{n} freie Knoten sein können. Diese Beschränkung erscheint wegen der besonderen Wahl der Knotenfolge (6.1) natürlich, da Informationen über das Meßintervall [a, b] zumeist vorliegen.

• Wegen (6.5) gilt l ≤ n − k.

An Stelle des Minimierungsproblems (6.4) erhalten wir jetzt

∥y − B(t̃)α∥^{2} ⟶ min .  (6.6)

Die freien Splineknoten t̃ können jedoch nicht beliebig gewählt werden. Es muß stets gelten

t_{j+k} > t_{j}  j = 1, …, n,  (6.7)

oder bei stärkeren Glattheitsforderungen sogar

t_{j+k−r} > t_{j}  j = r + 1, …, n.  (6.8)

Diese Nebenbedingungen für die Splineknoten werden wir im nächsten Abschnitt formulieren.

6.2 Anordnungsbedingung für die Splineknoten

Für eine Knotenfolge t muß gelten

t_{j+1} − t_{j} ≥ 0  j = 1, …, n + k − 1.

Da bei der Minimierung von (6.6) nur die freien Knoten t̃ verändert werden, lautet diese Bedingung

t_{p(j)} − t_{p(j)−1} ≥ 0
t_{p(j)+1} − t_{p(j)} ≥ 0  j = 1, …, l.

Für die numerische Praxis ist diese Bedingung ungünstig ( (6.7) bzw. (6.8) wird nicht erfüllt!). Wir fordern daher

t_{p(j)} − t_{p(j)−1} ≥ ϵ  (6.9)
t_{p(j)+1} − t_{p(j)} ≥ ϵ  j = 1, …, l.  (6.10)

mit einer vorgegebenen Konstanten 0 < ϵ ≪ 1. Die Bedingungen (6.9), (6.10) können äquivalent in Matrixform angegeben werden

Ct̃ ≥ h  (6.11)

mit C ∈ ℝ^{ncstr,l}, t̃ ∈ ℝ^{l}, h ∈ ℝ^{ncstr}. ncstr ist die Anzahl der linear unabhängigen Nebenbedingungen. Die Matrix C und der Vektor h hängen von der Lagebeziehung der Knoten ab und können leicht algorithmisch erzeugt werden. Zur Illustration möge ein kleines Beispiel dienen:

Beispiel: k = 4, n = 8, l = 3, p(1) = 5, p(2) = 6, p(3) = 8

\begin{matrix}t_{1} \\ t_{2} \\ t_{3} \\ t_{4} \end{matrix}  t_{5}  t_{6}  t_{7}  \begin{matrix}t_{9} \\ t_{10} \\ t_{11} \\ t_{12} \end{matrix} t_{8}

[ \begin{matrix}1 & 0 & 0 \\ −1 & 1 & 0 \\ 0 & −1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & −1 \end{matrix} ] ( \begin{matrix}t_{5} \\ t_{6} \\ t_{8} \end{matrix} ) ≥ ( \begin{matrix}t_{4} + ϵ \\ ϵ \\ −t_{7} + ϵ \\ t_{7} + ϵ \\ −t_{9} + ϵ \end{matrix} )

Ct̃ ≥ h  C ∈ ℝ^{5,3}, h ∈ ℝ^{5}

Bemerkung 28 • Es gilt l + 1 ≤ ncstr ≤ 2l.

• Prinzipiell läßt sich auch

t_{p(j)} − t_{p(j)−1} ≥ \underline{ϵ}_{j}
t_{p(j)+1} − t_{p(j)} ≥ \overline{ϵ}_{j}  j = 1, …, l

mit 2l Konstanten \underline{ϵ}_{j}, \overline{ϵ}_{j} fordern. Das Problem liegt dann in der Wahl geeigneter Konstanten. Deshalb haben wir im vorliegenden Programmpaket lediglich (6.9) und (6.10) mit einer, in Abhängigkeit von der durchschnittlichen Größe der Knoten zu wählenden Konstante ϵ realisiert.

• Die Matrix C enthält maximal zwei Nichtnullelemente in jeder Zeile.

• Mit der Wahl eines beliebig kleinen, aber positiven ϵ ist – zumindest theoretisch in rundungsfehlerfreier Arithmetik – die höchstmögliche Glattheit des Splines an den Stellen t_{p(j)}, j = 1, …, l, sichergestellt, da die freien Knoten dann einfache Knoten sind.

6.3 Untersuchung der Struktur des Problems

Jetzt können wir das vollständige Minimierungsproblem formulieren.

min \{ψ(t̃, α) := \frac{1}{2}∥y − B(t̃)α∥^{2} : Ct̃ ≥ h , t̃ ∈ ℝ^{l}, α ∈ ℝ^{n}\}  (6.12)

Bei diesem Problem handelt es sich um ein separables Quadratmittelproblem mit linearen Ungleichheitsnebenbedingungen. Die Splinekoeffizienten α ∈ ℝ^{n} treten linear auf, die freien Splineknoten t̃ ∈ ℝ^{l} treten nichtlinear auf. Wenn die Whitney-Schoenberg-Bedingung erfüllt ist, hat die Beobachtungsmatrix B Vollrang n und (6.12) liefert für jedes feste zulässige t̃ ∈ ℝ^{l} die eindeutige Lösung α_{opt}(t̃) = B(t̃)^{+}y.

6.3.1 Der Übergang zum reduzierten Funktional beim freien Problem

Wir betrachten zunächst das freie Optimierungsproblem

min \{ψ(t̃, α) := \frac{1}{2}∥y − B(t̃)α∥^{2} : t̃ ∈ ℝ^{l}, α ∈ ℝ^{n}\} .  (6.13)

Wir setzen voraus, daß die Whitney-Schoenberg-Bedingung überall erfüllt ist. Theorem 2 von Golub/Pereyra rechtfertigt dann den Übergang zum reduzierten Funktional ϕ(t̃)

min\{ϕ(t̃) := \frac{1}{2}∥\underbrace{[I − B(t̃)B(t̃)^{+}] y}_{=:G(t̃)}∥^{2} = \frac{1}{2}∥G(t̃)∥^{2} : t̃ ∈ ℝ^{l}\}.  (6.14)

6.3.2 Der Übergang zum reduzierten Funktional beim restringierten Problem

Die Zulässigkeit des Überganges zum reduzierten Funktional bezieht sich bis jetzt nur auf freie Quadratmittelprobleme. Bei dem zu lösenden Problem (6.12) handelt es sich aber um ein separables Quadratmittelproblem mit linearen Ungleichheitsnebenbedingungen an die nichtlinearen Variable t̃. Kaufman hat in [GP76] gezeigt, daß dieser Übergang ebenfalls zulässig ist.

Satz 13 (Kaufman, zitiert nach [Par85]) Ein Problem der Form

min \{ψ(t̃, α) := \frac{1}{2}∥y − B(t̃)α∥^{2} : Ct̃ ≥ h , t̃ ∈ ℝ^{l}, α ∈ ℝ^{n}\}

ist äquivalent mit der Lösung von

min \{ϕ(t̃) := \frac{1}{2}∥[I − B(t̃)B(t̃)^{+}] y∥^{2} = \frac{1}{2}∥P_{B}^{⊥}(t̃)∥^{2} : Ct̃ ≥ h , t̃ ∈ ℝ^{l}\}

gefolgt von

α^{∗} = B(t̃)^{+}y.

Bemerkung 29 • Die Aussage des Satzes ist anschaulich klar, da die variable Projektion ja bezüglich der linearen Variablen α erfolgt, und nicht bezüglich der restringierten nichtlinearen Variablen t̃.

• Den Fall von nichtlinearen Gleichungsnebenbedingungen an t̃ untersuchen [KP78]. Das Problem läßt sich dann in ein verallgemeinertes separables Quadratmittelproblem überführen und mit den Methoden aus Kapitel 5 lösen.

Zusätzliche Nebenbedingungen auch an die linearen Variablen erschweren die algorithmische Lösung sehr. Diese allgemeinen Probleme wurden in [Par85] untersucht. Dort findet sich auch eine Klassifizierung nach den verschiedenen Typen von Nebenbedingungen.

Wir können nun an Stelle von (6.12) das restringierte, reduzierte Problem

min\{ϕ(t̃) := \frac{1}{2}∥\underbrace{[I − B(t̃)B(t̃)^{+}] y}_{=:G(t̃)}∥^{2} = \frac{1}{2}∥G(t̃)∥^{2} : Ct̃ ≥ h , t̃ ∈ ℝ^{l}\} (6.15)

betrachten. Die effektive Lösung von (6.15) ist Gegenstand des restlichen Kapitels. Dabei sind die folgenden Besonderheiten zu beachten und – soweit möglich – auszunutzen.

• B hat stets Vollrang n und besitzt eine Bandbreite kleiner oder gleich k

• B′ hat ebenfalls eine Bandbreite kleiner oder gleich k

• B(t̃)^{+}y kann sehr effektiv mittels zeilenweiser Givens-Drehungen berechnet werden

• C hat maximal zwei Nichtnullelemente in jeder Zeile

Grundlage für die Lösung von (6.15) ist ein Verallgemeinertes Gauß-Newton-Ver-fahren.

6.4 Existenz einer Lösung des reduzierten Problems

Satz 14 (Existenz einer Lösung des reduzierten Problems) Die Knotenfolge t erfülle überall die Whitney-Schoenberg-Bedingung und es gelten (6.9), (6.10). Dann besitzt das reduzierte Problem (6.15) eine Lösung t̃^{∗}.

Beweis: Die Beobachtungsmatrix und die Moore-Penrose-Inverse sind unter den gestellten Voraussetzungen stetige – sogar differenzierbare – Funktionen der Splineknoten. Das stetige Funktional ϕ(t̃) nimmt daher über der abgeschlossenen und beschränkten (a ≤ t_{p(1)}, t_{p(l)} ≤ b) Menge {t̃ ∈ ℝ^{l} : Ct̃ ≥ h} nach dem Satz von Weierstraß ihr Minimum an. □

Die Charakterisierung und Berechnung einer eindeutigen Bestapproximation t̃^{∗} ist ein schwieriges Problem. Wir beschränken uns aus praktischen Gesichtspunkten daher auf die Bestimmung von lokal guten Approximationen. Unser Ziel ist es, die Lage der Knoten t̃ so zu verbessern, daß wir mit weniger, aber gut plazierten Knoten die notwendige Approximationsgüte erreichen. Es können also solange Knoten entfernt werden, bis die Approximationsgüte in der Größenordnung der Standardabweichung der Meßdaten ist.

6.5 Die Berechnung der Jacobi-Matrix

Zur Anwendung eines Verallgemeinerten Gauß-Newton-Verfahrens zur Lösung von (6.15) benötigen wir die Jacobi-Matrix G′ ∈ ℝ^{m,l} der reduzierten Funktion G :

ℝ^{l} ⟶ ℝ^{m} mit G(t̃) := [I − B(t̃)B(t̃)^{+}] y. Die reduzierte Funktion G(t̃) ist unter den gestellten Vorausssetzungen – Whitney-Schoenberg-Bedingung erfüllt, Knotenfolge genügt (6.9),(6.10) – differenzierbar und es gilt

G′(t̃)s = {A(t̃)[s] + (A(t̃)[s])^{T} } y
mit A(t̃)[s] := − [I − B(t̃)B(t̃)^{+}] (B′(t̃)[s]) B(t̃)^{+}.

Die Jacobi-Matrix G′ wird in dem entwickelten Programmpaket auf zwei unterschiedliche Arten approximiert.

6.5.1 Die Berechnung der Jacobi-Matrix im Golub/Pereyra-Modell

Im Golub/Pereyra-Modell wird das reduzierte Funktional ϕ(t̃ + s) durch die quadratische Modellfunktion

μ_{GP} (t̃ + s) := \frac{1}{2} ∥G(t̃) + G′(t̃)s∥^{2}

ersetzt. Wegen des Auftretens des Tensors A(t̃)^{T} eignet sich dieses Modell nur bedingt für die numerische Approximation der Ableitungen durch Differenzenquotienten. Wir approximieren daher direkt

G′(t̃)e^{k} ≈ \frac{G(t̃ + h_{k}e^{k}) − G(t̃)}{h_{k}}

Diese Approximation benötigt l + 1 Funktionsberechnungen von G.

6.5.2 Die Berechnung der Jacobi-Matrix im Kaufman-Modell

Im Kaufman-Modell wird das reduzierte Funktional ϕ(t̃+s) durch die quadratische Modellfunktion

μ_{K}(t̃ + s) := \frac{1}{2} ∥G(t̃) + J′(t̃)s∥^{2}
mit G′(t̃)s ≈ J(t̃)s := {A(t̃)[s]} y
= − [I − B(t̃)B(t̃)^{+}] (B′(t̃)[s]) B(t̃)^{+}y

ersetzt. Die Approximation J ∈ ℝ^{m,l} der Jacobi-Matrix G′ wird spaltenweise berechnet. Wir erhalten

Je^{k} = − [I − B(t̃)B(t̃)^{+}] (B′(t̃)[e^{k}]) B(t̃)^{+}y
B(t̃) = B_{j} (x_{i}, t̃) e^{i} ∘ e^{j}  i = 1, …, m; j = 1, …, n
B′(t̃) = \frac{∂B_{j} (x_{i}, t̃)}{∂t_{p(k)}} e^{i} ∘ e^{j} ∘ e^{k}  k = 1, …, l

B′(t̃) [e^{k}] = \frac{∂B_{j} (x_{i}, t̃)}{∂t_{p(k)}} δ_{k,k} e^{i} ∘ e^{j}.

Hierbei haben wir von der Einsteinschen Summenkonvention und der Darstellung eines Tensors als dyadisches Produkt Gebrauch gemacht. Die Matrix B′(t̃) [e^{k}] wird numerisch über Differenzenquotienten berechnet.

(B′(t̃) [e^{k}])_{i,j} ≈ \frac{B_{j} (x_{i}, t̃ + h_{k}e^{k}) − B_{j} (x_{i}, t̃)}{h_{k}}

Man beachte dabei die Blockbandstruktur von B!

Bemerkung 30 • Ausdrücke der Form B(t̃)^{+}u werden stets über die QR-Faktorisierung von B(t̃) berechnet und niemals über die direkte Bildung von B(t̃)^{+}.

• B(t̃) wird nur einmal faktorisiert. Die Drehungsparameter werden dann abgespeichert. Die wiederholte Berechnung von B(t̃)u für verschiedene u kann dann sehr effektiv erfolgen.

• Die Approximation von B′(t̃) [e^{k}] wird sehr effektiv durchgeführt, da der verwendete Algorithmus zur Berechnung von B-Splines simultan alle nichtverschwindenen B-Splines an der Stelle x_{i} berechnet.

6.6 Ein Verfahren zur Lösung des Problems

Jetzt können wir ein Verfahren zur Lösung des restringierten, reduzierten Problems

min\{ϕ(t̃) := \frac{1}{2}∥\underbrace{[I − B(t̃)B(t̃)^{+}] y}_{=:G(t̃)}∥^{2} = \frac{1}{2}∥G(t̃)∥^{2} : Ct̃ ≥ h , t̃ ∈ ℝ^{l}\} (6.16)

angeben. Wir verwenden ein gedämpftes Verallgemeinertes Gauß-Newton-Verfahren mit der in 4.2.5 beschriebenen Regularisierung. Wir formulieren das Verfahren noch einmal mit den aktuellen Bezeichnungen.

Gedämpftes Verallgemeinertes Gauß-Newton-Verfahren mit Regularisierung

S1: Wähle (zulässigen) Startwert t̃^{0} ∈ ℝ^{l}, α ∈ (0, 1), δ ∈ (0, \frac{1}{2})

Setze k := 0; γ_{−1} := 1;

S2: Iteration

Repeat

S2.1: G := G(t̃^{k}) = [I − B(t̃^{k})B(t̃^{k})^{+}] y; Berechnung des Residiums

J :≈ G′(t̃^{k}); Berechnung einer Approximation der Jacobi-Matrix (Golub/Pereyra- oder Kaufman-Modell)

S2.2: Bestimme Abstiegsrichtung s^{k} aus

∥G + Js∥^{2} ⟶ min
mit  C (t̃^{k} + s) ≥ h

Falls das Problem schlecht konditioniert ist, löse das regularisierte Problem

∥( \begin{matrix}G \\ 0 \end{matrix} ) + [ \begin{matrix}J \\ \sqrt{λ}I \end{matrix} ] s∥^{2} ⟶ min
mit  C (t̃^{k} + s) ≥ h

S2.3: γ̃_{k} := γ_{k}; t̃^{+} := t̃^{k} + γ̃_{k}s^{k}; G^{+} := G(t̃^{+});

S2.4: Abstiegstest

if \frac{1}{2} ∥G∥^{2} − \frac{1}{2} ∥G + γ̃_{k}Js^{k}∥^{2} ≥ −δγ̃_{k}G^{T} Js^{k}

then

t̃^{k+1} := t̃^{+};

γ_{k+1} := min(γ_{k}/α, 1);

k := k + 1;

else

γ_{k} := γ_{k} ∗ α; goto S2.3;

until (k > kmax) or Konvergenz

Bemerkung 31 • Die Lösung des linearen Quadratmittelproblems mit linearen Ungleichheitsnebenbedingungen in Schritt S2.2 erfolgt mittels des Algorithmus LSI (siehe 3.2 ).

• Die Feststellung, wann ein Problem in S2.2 in unserem Sinne schlecht konditioniert ist, sowie die Wahl des Regularisierungsparameters λ ist in 4.2.5 beschrieben.

• Für das Abbruchkriterium gilt das in 4.1.2 zum Gauß-Newton-Verfahren gesagte. Insbesondere ist bei restringierten Problemen der, bei freien Problemen mit ∥G′(t̃^{k})G(t̃^{k})∥ ≤ ε_{3}^{a} äquivalente, Test (|G(t̃^{k})^{T} G′(t̃^{k})s| ≤ ε_{4}^{a}) wichtig. Er stellt fest, ob noch eine zulässige Abstiegsrichtung s gefunden werden kann.

• Für die numerischen Tests wählen wir δ = 0.2 und α = 0.2.

6.7 Testbeispiele

Dieser Abschnitt erhält einige Beispiele aus der Literatur sowie ein selbstgewähltes Beispiel, welche die Leistungsfähigkeit des entwickelten Verfahrens zeigen sollen. Alle Testrechnungen wurden auf einem 486-PC in IEEE double Arithmetik mit einer relativen Maschinengenauigkeit von macheps ≈ 2.22E − 16 durchgeführt. Die Abbruchkonstanten^{1} (siehe 4.1.2) wurden wie folgt gewählt:

\begin{matrix}ε_{0}^{a}= & 1.0E − 06 \\ ε_{1}^{a}= & 1.0E − 04 \\ ε_{3}^{a}= & 1.0E − 10 \\ ε_{4}^{a}= & 1.0E − 10 \end{matrix}

Zum Vergleich des entwickelten Verfahrens mit Standardoptimierungsverfahren führen wir noch eine Minimierung des reduzierten Funktionals mit der Matlab-Pro-zedur constr durch. Die Abbruchkonstanten wurden auf die vorgegebenen Standardwerte gesetzt.

1 Die entwickelten Verallgemeinerten Gauß-Newton-Verfahren liefern folgende Return Codes:

1  ∥G∥ ≤ ε_{0}^{a} 2  ∥J^{T} G∥ ≤ ε_{3}^{a} 3  |G^{T} Jp| ≤ ε_{4}^{a} 4 ∥x^{k+1} − x^{k}∥ ≤ ε_{1}^{a} (∥x^{k}∥ + ε_{2}^{a})

5  k > kmax 6  failed otherwise

6.7.1 Titan-Daten

Unser erstes Testbeispiel sind die, in der Spline-Literatur bekannten Titanium Heat Data (siehe [dB78]). Sie beschreiben eine Eigenschaft von Titan in Abhängigkeit von der Temperatur. m = 49 Meßpunkte sind äquidistant im Intervall [595, 1075] verteilt. Die Meßwerte verlaufen relativ flach bis auf den typischen Peak im Bereich x = 900. Wir wollen nun die m = 49 Datenpunkte durch n = 11 B-Splines der Ordnung k = 4 approximieren. Zunächst wählen wir die l = 7 inneren Splineknoten äquidistant im Intervall [595, 1075] verteilt, was angesichts des Peaks sicher keine sehr gute Wahl ist. ^{2} Die Konstante ϵ, welche den minimalen Abstand zwischen den freien Splineknoten angibt, wurde auf ϵ = 0.0 für das Kaufman-Modell und ϵ = 0.1 für die Matlab-Routine gesetzt.

Bemerkung 32 Die theoretische Durchführbarkeit ist zwar nur für den Fall ϵ > 0.0 gesichert. Die von uns entwickelte Routine ist aber relativ unempfindlich gegenüber der Wahl von ϵ, während die Matlab-Routine bei zu kleinem ϵ wegen der Rundungsfehler abbricht.

Startpunkt  Matlab-constr Kaufman-Modell
t_{5}  655  656.648608  603.122987
t_{6}  715  722.545037  805.055831
t_{7}  775  843.889596  866.128923
t_{8}  835  874.837667  866.128923
t_{9}  895  897.730378  897.867208
t_{10}  955  916.990605  916.774224
t_{11}  1015  972.317384  973.773457
∥r∥  7.924657E − 01 7.370720E − 02 5.651648E − 02
Schritte  484 Func. calls  24
Zeit  12.037 MFlops  0.56 sec
|G^{T} Js|  2.30011E − 11
∥J^{T} G∥  3.98777E − 08
Return Code  erfolgreich  3

Ein erfahrener Nutzer würde natürlich die Splineknoten geschickter wählen. Wir beginnen daher die Iteration noch einmal mit einer guten Knotenfolge.

Startpunkt  Matlab-constr Kaufman-Modell
t_{5}  750  802.354053  802.345327
t_{6}  825  865.366264  859.461095
t_{7}  885  865.366264  879.534951
t_{8}  890  905.535024  879.534951
t_{9}  920  905.635024  910.350046
t_{10}  935  934.584027  910.350046
t_{11}  975  969.940811  974.653448
∥r∥  1.054318E − 01 5.875918E − 02 3.934182E − 02
Schritte  704 Func. calls  17
Zeit  17.706 MFlops  0.40 sec
|G^{T} Js|  5.93399E − 11
∥J^{T} G∥  1.49747E − 08
Return Code  Abbruch  3

2Die Splineknoten t_{j} werden in den folgenden Abbildungen durch Rechtecke, die Meßdaten (x_{i}, y_{i}) durch Kreise dargestellt.

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.68; 0.85; 1.02; 1.18; 1.35; 1.52; 1.68; 1.85; 2.02; E 0]

Abbildung 1: Titanium Heat Data, k = 4, n = 11, equidistant knots

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.75; 0.92; 1.08; 1.24; 1.41; 1.57; 1.73; 1.89; 2.06; E 0]

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.75; 0.92; 1.08; 1.24; 1.40; 1.57; 1.73; 1.89; 2.05; E 0]

Abbildung 3: Titanium Heat Data, k = 4, n = 11, l = 7, good startpoint

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.75; 0.91; 1.07; 1.23; 1.39; 1.55; 1.72; 1.88; 2.04; E 0]

6.7.2 Titan-Daten II

Ein Motiv für die Entwicklung des vorgestellten Algorithmus war die bessere Plazierung der freien Knoten, um anschließend gegebenfalls die Anzahl der Knoten zu reduzieren. Diese Knotenreduktion möchten wir mit dem folgenden Beispiel demonstrieren. Wir approximieren erneut die Titan-Daten. Diesmal wählen wir allerdings n = 17 B-Splines der Ordnung k = 4. Wir starten die Approximation wieder mit äquidistanten inneren Knoten (l = 13).

Startpunkt  Matlab-constr Golub/Pereyra-Modell
t_{5}  629.285714  626.746755  595.017590
t_{6}  663.571428  662.988340  595.020012
t_{7}  697.857142  697.976877  694.798680
t_{8}  732.142857  732.369775  740.838863
t_{9}  766.428571  766.966956  826.246672
t_{10}  800.714285  803.464577  863.392120
t_{11}  835.000000  866.069225  863.459128
t_{12}  869.285714  866.169225  863.597505
t_{13}  903.571428  897.564015  897.564233
t_{14}  937.857142  917.369575  917.109542
t_{15}  972.142857  967.743774  971.160283
t_{16}  1006.428571  1005.407904  994.053308
t_{17}  1040.714285  1040.315564  1066.665752
∥r∥  3.237391E − 01 5.468452E − 02  5.235772E − 02
Schritte  924 Func. calls  39
Zeit  26.912 MFlops  1.65 sec
|G^{T} Js|  9.25177E − 11
∥J^{T} G∥  4.89054E − 08
Return Code  erfolgreich  3

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.72; 0.89; 1.05; 1.21; 1.37; 1.54; 1.70; 1.86; 2.02; E 0]

Abbildung 5: Titanium Heat Data, k = 4, n = 17, equidistant knots

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.75; 0.92; 1.08; 1.24; 1.40; 1.57; 1.73; 1.89; 2.06; E 0]

Wir modifizieren jetzt die Knotenfolge, welche wir als Ergebnis der Optimierung mit dem Golub/Pereyra-Modell erhalten haben.

Knoten  Modifikation
t_{5}, t_{6}  entfernen
t_{10}, t_{11}, t_{12} ein einfacher Knoten \frac{t_{10}+t_{11}+t_{12}}{3}
t_{15}, t_{17}  entfernen

Statt n = 17 erhalten wir n = 11 Knoten. Auf diese neue Knotenfolge wenden wir erneut das Minimierungsverfahren an.

Startpunkt  Golub/Pereyra-Modell
t_{5}  694.798680  604.372097
t_{6}  740.838863  805.055833
t_{7}  826.246672  866.128905
t_{8}  863.482918  866.128905
t_{9}  897.564233  897.867279
t_{10}  917.109542  916.774162
t_{11}  994.053308  973.773496
∥r∥  2.863939E − 01  5.651648E − 02
Schritte  18
Zeit  0.50 sec
|G^{T} Js|  2.86444E − 11
∥J^{T} G∥  4.44755E − 08
Return Code  3

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.74; 0.90; 1.06; 1.22; 1.38; 1.54; 1.70; 1.86; 2.03; E 0]

Abbildung 7: Titanium Heat Data, k = 4, n = 11, Number of knots reduced

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.75; 0.92; 1.08; 1.24; 1.41; 1.57; 1.73; 1.89; 2.06; E 0]

6.7.3 Spektraldaten

Wir betrachten Meßwerte, welche durch die Spektralfunktion

f(x) = 7 × exp (−\frac{(x − 5)^{2}}{0.25} ) − 5 × exp (−(x + 3)^{2}) + 0.1 × x

beschrieben werden. Mittels Pseudozufallszahlen ϵ_{i}, −0.05 ≤ ϵ_{i} ≤ 0.05, erzeugen wir die Werte

y_{i} = f (x_{i}) + ϵ_{i}

an m = 90 äquidistanten Punkten im Intervall [−10, 10]. Für den verwendeten Meßwertsatz gilt E(ϵ) = 8.3E − 04 und D^{2}(ϵ) = 3.0E − 03.

Bemerkung 33 Wenn man den physikalischen Hintergrund der Daten kennt, ist natürlich das ebenfalls separable Parameterschätzproblem mit der Ansatzfunktion

r(x, t̃, α) = ∑_{j=1}^{k} α_{j} exp (− (\frac{x − t̃_{j}}{t̃_{j+k}} )^{2}) + α_{k+1} + α_{k+2}x + α_{k+3}x^{2}

der Problemstellung besser angepaßt. Wir stellen uns jedoch auf den Standpunkt, daß nichts über den physikalischen Hintergrund bekannt ist.

Diese Daten wollen wir durch n = 25 B-Splines der Ordnung k = 5 approximieren. Wir wählen l = 20 äquidistante innere Splineknoten als Ausgangspunkt der Iteration. Wir verzichten diesmal auf die numerische Angabe der berechneten Knoten.

Startpunkt  Kaufman-Modell Golub/Pereyra-Modell
∥r∥  3.796930E + 00 2.143237E − 01  1.445334E + 00
Schritte  34  50
Zeit  5.01 sec  8.94 sec
|G^{T} Js|  8.77864E − 11  1.06799E − 03
∥J^{T} G∥  4.22577E − 08  4.53193E + 00
Return Code  3  5

[figure: -8.0; -6.0; -4.0; -2.0; 0.00; 2.00; 4.00; 6.00; 8.00; E 0; -4.0; -2.7; -1.4; -0.2; 1.08; 2.37; 3.66; 4.95; 6.25; E 0]

Abbildung 9: Spectral Data, k = 5, n = 25, equidistant knots

[figure: -8.0; -6.0; -4.0; -2.0; 0.00; 2.00; 4.00; 6.00; 8.00; E 0; -4.0; -2.7; -1.4; -0.1; 1.13; 2.43; 3.73; 5.03; 6.33; E 0]

Abbildung 10: Spectral Data, k = 5, n = 25, l = 20, Kaufman-Model

[figure: -8.0; -6.0; -4.0; -2.0; 0.00; 2.00; 4.00; 6.00; 8.00; E 0; -4.0; -2.6; -1.3; -0.0; 1.30; 2.64; 3.98; 5.31; 6.65; E 0]

Abbildung 11: Spectral Data, k = 5, n = 25, l = 20, Golub/Pereyra-Model

7 Splineglättung mit freien Knoten

Bei der praktischen Anwendung von Approximationsschemata ist oftmals ein Kompromiß zwischen der Güte der Approximation und der Glattheit gefragt. Neben dem Approximationsterm

φ(t̃, α) := ∥y − B(t̃)α∥^{2}  (7.1)

betrachten wir jetzt zusätzlich einen sogenannten Glättungsterm

ρ̄(t̃, α) := ∥s^{(r)}∥_{L_{2}}^{2} = ∫_{a}^{b} [s^{(r)}(x)]^{2} dx  (7.2)

mit festem r < k und das zugehörige Minimierungsproblem

φ(t̃, α) + μρ̄(t̃, α) = ∥y − B(t̃)α∥^{2} + μ ∫_{a}^{b} [s^{(r)}(x)]^{2} dx ⟶ min . (7.3)

Der Glättungsparameter μ > 0 liefert im Grenzfall μ → +0 die Approximation mit kleinster Seminorm φ und im Grenzfall μ → +∞ die beste polynomiale Approximation der Ordnung r. Die Bestimmung des Glättungsparameters μ kann in den Optimierungsprozeß einbezogen werden ( siehe [Wah90] für die entsprechenden Verfahren und den statistischen Hintergrund). Wir werden dagegen in dieser Arbeit den Glättungsparameter interaktiv durch den Nutzer vorgeben lassen. Wir werden zunächst den Glättungsterm (7.2) in Abhängigkeit von den Splinekoeffizienten α und den Knoten t ausdrücken. Anschließend wird dieser exakte Glättungsterm durch eine einfacher zu berechnende Approximation ersetzt. Die Darstellung in den beiden folgenden Abschnitten sowie die Algorithmen zur Berechnung eines glättenden Splines folgen [Kun89] und [SK93].

7.1 Darstellung des Glättungsterms

Für die Knotenfolge t gelte

t_{j+k−r} > t_{j}  j = r + 1, …, n.

Dann existiert s^{(r)}(x) und besitzt die Darstellung

s^{(r)}(x) = ∑_{j=r+1}^{n} B_{j,k−r}(x)α_{j}^{(r)}
mit α^{(r)} := D_{r}α
D_{r} := H_{r}L_{r} … H_{1}L_{1} ∈ ℝ^{n−r,n}
H_{ν} := (k − ν) diag (\frac{1}{t_{k+j−ν} − t_{j}} )_{j=1+ν,...,n} ∈ ℝ^{n−ν,n−ν}

L_{ν} := [ \begin{matrix}−1 & 1 \\ & −1 & 1 \\ & & ... & ... \\ & & & −1 & 1 \\ & & & & −1 & 1 \end{matrix} ] ∈ ℝ^{n−ν,n−ν+1} ν = 1, …, r.

Setzt man dies in den Glättungsterm ρ̄ ein, so erhält man

ρ̄ = ∫_{a}^{b} [ ∑_{i=r+1}^{n} B_{i,k−r}(x)α_{i}^{(r)}][ ∑_{j=r+1}^{n} B_{j,k−r}(x)α_{j}^{(r)}] dx

= ∑_{i,j=r+1}^{n} α_{i}^{(r)}α_{j}^{(r)} ∫_{a}^{b} B_{i,k−r}(x)B_{j,k−r}(x) dx
= α^{(r)^{T}} M̄_{r}α^{(r)}
mit  M̄_{r} = (m̄_{ij}^{(r)}) ∈ ℝ^{n−r,n−r}
m̄_{ij}^{(r)} := ∫_{a}^{b} B_{i,k−r}(x)B_{j,k−r}(x) dx i, j = r + 1, …, n.

Bemerkung 34 • M̄_{r} ist eine Bandmatrix der Bandbreite 2(k − r) − 1.

• M̄_{r} ist als Gramsche Matrix linear unabhängiger Funktionen ( t_{j+k−r} > t_{j}, j = r + 1, …, n!) positiv definit.

Wir führen die Cholesky-Faktorisierung der symmetrischen, positiv definiten Matrix M̄_{r} durch.

M̄_{r} = F̄_{r}^{T} F̄_{r}
F̄_{r} ∈ ℝ^{n−r,n−r}  reguläre obere Dreiecksmatrix der Bandbreite k − r

Damit erhält man für den Glättungsterm

ρ̄ = α^{(r)^{T}} M̄_{r}α^{(r)}
= α^{T} D_{r}^{T} F̄_{r}^{T} F̄_{r}D_{r}α
= ∥S̄_{r}α∥^{2}

mit der Matrix

S̄_{r} = S̄_{r}(t̃) := F̄_{r}D_{r} ∈ ℝ^{n−r,n}.

S̄_{r} ist eine obere Dreiecksmatrix und besitzt die Bandbreite k. Die Berechnung der Elemente m̄_{ij} der Gramschen Matrix ist für k − r > 2 sehr teuer. Sie kann etwa, wie von [Sch81] empfohlen, über Gaußsche Quadraturformeln erfolgen. Wir verwenden eine billigere Approximation des Glättungstermes, welche sich in den Anwendungen (siehe [SK93]) als günstig erwiesen hat.

7.2 Approximation des Glättungsterms

Wir ersetzen die L_{2}-Norm in der Definition des Glättungsterms durch ihr diskretes Äquivalent im Splineraum S_{k−r,t}.

ρ̃(t̃, α) := ∥s^{(r)}∥_{l_{2}}^{2} = ∑_{j=r+1}^{n} (α_{j}^{(r)})^{2} \frac{t_{j+k−r} − t_{j}}{k − r}  (7.4)

Die Normen sind in folgendem Sinne äquivalent:

Es existiert eine Konstante d_{k−r}, welche nicht von der Knotenfolge t abhängt, so daß gilt:

d_{k−r}^{−1} ∥s^{(r)}∥_{l_{2}} ≤ ∥s^{(r)}∥_{L_{2}} ≤ ∥s^{(r)}∥_{l_{2}} ∀s ∈ S_{k,t}.

Die Existenz der Konstanten d_{k−r} für beliebige L_{p}-Normen und ihr diskretes Äquivalent l_{p} geht auf de Boor [dB76] zurück. Die Ersetzung der L_{2}-Norm durch die diskrete Norm l_{p} wird ebenfalls in [LM88] und [ADLM90] vorgeschlagen. Für diese Approximation erhalten wir

ρ̃ = α^{(r)^{T}} F̃_{r}^{T} F̃_{r}α^{(r)}
= α^{T} D_{r}^{T} F̃_{r}^{T} F̃_{r}D_{r}α

mit F̃_{r} := diag (\sqrt{\frac{t_{j+k−r}−t_{j}}{k−r}} )_{j=r+1,...,n}, und schließlich

ρ̃ = ∥S̃_{r}α∥^{2}

mit der oberen Dreiecksmatrix

S̃_{r} = S̃_{r}(t̃) := F̃_{r}D_{r} ∈ ℝ^{n−r,n}

der Bandbreite r + 1 (im Gegensatz zu k bei S̄_{r}!). Wegen der Äquivalenz der Normen besitzen die Matrizen S̄_{r} und S̃_{r} beide vollen Rang n − r und haben denselben Nullraum

ker (S̃_{r}) = ker (S̄_{r}) = {α ∈ ℝ^{n} : α_{j}^{(r)} = 0 , j = r + 1, …, n} ,

welcher durch alle Polynome der Ordnung r aus S_{k,t} gekennzeichnet ist.

7.3 Formulierung des kompletten Funktionals

Mit ρ = ρ̄ = ∥S̄_{r}α∥^{2}, S = S̄_{r} (exakter Glättungsterm) oder ρ = ρ̃ = ∥S̃_{r}α∥^{2}, S = S̃_{r} (Approximation des Glättungsterms) erhalten wir das zu minimierende Funktional

φ(t̃, α) + μρ(t̃, α) = ∥y − B(t̃)α∥^{2} + μ ∥S_{r}α]^{2}  (7.5)
= ∥( \begin{matrix}y \\ 0 \end{matrix} ) − [ \begin{matrix}B(t̃) \\ \sqrt{μ}S(t̃) \end{matrix} ] α∥^{2} ⟶ min .  (7.6)

Zusammen mit den Anordnungsbedingungen Ct̃ ≥ h ergibt sich das separable Quadratmittelproblem

min \{ψ(t̃, α) := \frac{1}{2} ∥( \begin{matrix}y \\ 0 \end{matrix} ) − [ \begin{matrix}B(t̃) \\ \sqrt{μ}S(t̃) \end{matrix} ] α∥^{2} : Ct̃ ≥ h , t̃ ∈ ℝ^{l}, α ∈ ℝ^{n}\} .

(7.7)

Wir setzen im folgenden voraus, daß

m ≥ r.  (7.8)

Dann folgt aus Bα = 0 und S = \{ S̄_{r} ; S̃_{r} \} = 0, daß α = 0, d.h. die Matrix [ \begin{matrix}B \\ \sqrt{μ}S \end{matrix} ]

hat für jedes μ > 0 Vollrang und (7.5) besitzt eine eindeutige Lösung α_{opt}(t̃).

Bemerkung 35 Die Bedingung (7.8) ersetzt gewissermaßen die Schoenberg-Whitney-Bedingung. Die, im Minimierungsprozeß schwer zu realisierende, Forderung nach einer Vollrangmatrix B entfällt durch die Regularisierung (7.8) und wird durch die, in unseren Anwendungen immer erfüllte, Bedingung m ≥ r ersetzt.

Unter der Voraussetzung (7.8) können wir zum reduzierten, restringierten Problem

min \{ϕ(t̃) := \frac{1}{2} ∥\{[ \begin{matrix}I_{m} \\ I_{n−r} \end{matrix} ] − [ \begin{matrix}B(t̃) \\ \sqrt{μ}S(t̃) \end{matrix} ] [ \begin{matrix}B(t̃) \\ \sqrt{μ}S(t̃) \end{matrix} ]^{+}\} ( \begin{matrix}y \\ 0 \end{matrix} )∥^{2}

= \frac{1}{2} ∥G(t̃)∥^{2} : Ct̃ ≥ h , t̃ ∈ ℝ^{l}\}  (7.9)

mit der differenzierbaren Abbildung G : ℝ^{l} ⟶ ℝ^{m+n−r}

G(t̃) := \{[ \begin{matrix}I_{m} \\ I_{n−r} \end{matrix} ] − [ \begin{matrix}B(t̃) \\ \sqrt{μ}S(t̃) \end{matrix} ] [ \begin{matrix}B(t̃) \\ \sqrt{μ}S(t̃) \end{matrix} ]^{+}\} ( \begin{matrix}y \\ 0 \end{matrix} )  (7.10)

übergehen. Dieses Problem kann mit dem gedämpften Verallgemeinerten Gauß-Newton-Verfah-

ren aus Abschnitt 6 gelöst werden, indem man dort B(t̃) durch [ \begin{matrix}B(t̃) \\ \sqrt{μ}S(t̃) \end{matrix} ] und y

durch ( \begin{matrix}y \\ 0 \end{matrix} ) ersetzt.

Die Berechnung der Jacobi-Matrix im Kaufman-Modell ist wegen der komplizierten Gestalt von S(t̃) schwierig. Da die Testrechnungen im Fall der Quadratmittelapproximation gezeigt haben, daß Golub/Pereyra-Modell und Kaufman-Modell nahezu gleichwertig sind, verwenden wir die einfacher zu berechnende Approximation

G′(t̃)e^{k} ≈ \frac{G(t̃ + h_{k}e^{k}) − G(t̃)}{h_{k}}

nach dem Golub/Pereyra-Modell. Diese Approximation benötigt l+1 Funktionsberechnungen von G. Es ist also lediglich die Berechnung eines glättenden Splines sicherzustellen. Dazu verwenden wir ein existierendes Verfahren aus [Kun89]. Bei diesem Verfahren wird zuerst B mittels zeilenweiser GIVENS-Drehungen in eine obere Dreiecksmatrix R_{B} transformiert.

Im nächsten Schritt wird die erweiterte Matrix [ \begin{matrix}R_{B} \\ \sqrt{μ}S \end{matrix} ] auf obere Dreiecksform

gebracht. Dazu wird die Matrix mit Permutationen auf Standardform gebracht und mittels GIVENS-Drehungen faktorisiert. Die Bandbreite k bleibt bei diesem Prozeß erhalten (siehe auch [Eld84]).

7.4 Testbeispiele

Ein Verfahren zur Berechnung eines glättenden Splines mit freien Knoten wurde unter Verwendung der Approximation des Glättungsterms mittels des Golub/Pereyra-Modells implementiert. Die Aussagen zur Wahl der Abbruchkonstanten des Verallgemeinerten Gauß-Newton-Verfahrens aus Abschnitt 6 gelten weiterhin. Bei der Berechnung eines glättenden Splines ist jedoch darauf zu achten, daß t_{j} < t_{j+k−r}, j = r +1, …, n gilt. Wir wählen deshalb für die Konstante ϵ, welche den minimalen Abstand zwischen den freien Knoten angibt, den Wert ϵ = 0.1. Für kleine Werte des Glättungsparameters μ > 0 erhält man annähernd die gleichen Ergebnisse wie in Abschnitt 6.7. Der Aufwand ist dabei natürlich wegen der größeren Dimension von G entsprechend größer.

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.69; 0.86; 1.03; 1.19; 1.36; 1.52; 1.69; 1.85; 2.02; E 0]

Abbildung 12: Titanium Heat Data, r = 2, μ = 1E3, equidistant knots

7.4.1 Titan-Daten III

Wir wenden uns ein letztes Mal den Titan-Daten zu. Erneut wollen wir die m = 49 Meßpunkte durch n = 11 B-Splines der Ordnung k = 4 approximieren. Für den Glättungsparameter μ wählen wir μ = 1E3. Wir beginnen den Minimierungsprozeß wieder mit äquidistanten inneren Knoten. Man vergleiche die folgenden Ergebnisse mit den Ergebnissen aus Abschnitt 6!

Startpunkt  Golub/Pereyra-Modell
t_{5}  655  627.239439
t_{6}  715  657.590926
t_{7}  775  801.398888
t_{8}  835  869.748225
t_{9}  895  897.358543
t_{10}  955  928.206791
t_{11}  1015  994.289363
∥r∥ = ∥y − Bα∥  8.151829E − 01  4.318360E − 01
(∥y − Bα∥^{2} + μ ∥Sα∥^{2})^{1/2} 9.219330E − 01  7.045600E − 01
Schritte  63
Zeit  2.15 sec
|G^{T} Js|  5.85493E − 11
∥J^{T} G∥  3.45285E − 08
Return Code  3

Zum Vergleich ist an der ermittelten Optimalstelle der Spline ohne Glättungsterm (μ = 0.0) berechnet worden.

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.73; 0.89; 1.05; 1.21; 1.38; 1.54; 1.70; 1.86; 2.02; E 0]

Abbildung 13: Titanium Heat Data, l = 7, r = 2, μ = 1E3, Golub/Pereyra-Model

[figure: 64.3; 69.1; 73.9; 78.7; 83.5; 88.3; 93.1; 97.9; 102.; E 1; 0.73; 0.89; 1.05; 1.21; 1.38; 1.54; 1.70; 1.86; 2.02; E 0]

Abbildung 14: Titanium Heat Data, μ = 0.0 7.4.2 RPN-14-Daten

Zum Abschluß wollen wir die RPN 14 Daten von Fritsch/Carlson [FC80] approximieren. Diese Daten stammen von einem radioaktiven Prozeß und zeigen eine für physikalisch-chemische Prozesse typische Sättigungserscheinung.

\begin{matrix}i & x_{i} & y_{i} \\ 1 & 7.99 & 0.00000E + 00 \\ 2 & 8.09 & 2.76429E − 05 \\ 3 & 8.19 & 4.37498E − 02 \\ 4 & 8.70 & 1.69183E − 01 \\ 5 & 9.20 & 4.69428E − 01 \\ 6 & 10.00 & 9.43740E − 01 \\ 7 & 12.00 & 9.98636E − 01 \\ 8 & 15.00 & 9.99919E − 01 \\ 9 & 20.00 & 9.99994E − 01 \end{matrix}

Die m = 9 Meßpunkte wollen wir durch quadratische B-Splines (k = 3) approximieren. Für die l = 2 inneren Splineknoten wählen wir t_{4} = 9 und t_{5} = 12. Es sei r = 1 und μ = 1E − 02.

Startpunkt  Golub/Pereyra-Modell
t_{4}  9.000000  9.464841
t_{5}  12.000000  10.341722
∥r∥ = ∥y − Bα∥  1.725484E − 01  2.752452E − 02
(∥y − Bα∥^{2} + μ ∥Sα∥^{2})^{1/2} 1.917395E − 01  9.099074E − 02
Schritte  8
Zeit  0.02 sec
|G^{T} Js|  9.42421E − 12
∥J^{T} G∥  3.70994E − 07
Return Code  3

[figure: 9.19; 10.3; 11.5; 12.7; 13.9; 15.1; 16.3; 17.5; 18.7; E 0; 0.10; 0.22; 0.33; 0.44; 0.56; 0.67; 0.78; 0.90; 1.01; E 0]

Abbildung 15: RPN14Data, k = 3, n = 5, r = 1, μ = 1E − 2, start point

[figure: 9.19; 10.3; 11.5; 12.7; 13.9; 15.1; 16.3; 17.5; 18.7; E 0; 0.10; 0.20; 0.30; 0.40; 0.50; 0.60; 0.70; 0.80; 0.90; E 0]

Abbildung 16: RPN14Data, l = 2, r = 1, μ = 1E − 2, Golub/Pereyra-Model

Literatur

[ADLM90] E. Arge, M. Dæhlen, T. Lyche, and K. Mørken. Constrained spline approximation of functions and data based on constrained knot removal. In J. C. Mason and M. G. Cox, editors, Algorithms for Approximation II, pages 4–20. Chapman and Hall, 1990.

[Bjö90] Å. Björck. Least squares methods. In P. G. Ciarlet and J. L. Lions, editors, Handbook of Numerical Analysis, volume I, Solution of equations in R^{n}, pages 589–652. North Holland, 1990.

[Boc87] H. G. Bock. Randwertproblemmethoden zur Parameteridentifizierung in Systemen nichtlinearer Differentialgleichungen (Habilitationschrift). Bonner Mathematische Schriften Nr.183, Universität Bonn, 1987.

[CMSW79] A. K. Cline, C. B. Moler, G. W. Stewart, and J. H. Wilkenson. An estimate for the condition of a matrix. SIAM J. Numer. Anal., 16:368– 375, 1979.

[Cox81] M. G. Cox. The least squares solution of overdetermined linear equations having band or augmented band structure. IMA Journal Numerical Analysis, 1:3–22, 1981.

[Cox82] M. G. Cox. Practical spline approximation. In P. R. Turner, editor, Topics in Numerical Analysis, pages 79–112. Springer Verlag, 1982.

[Cox89] M. G. Cox. Linear algebra support modules for approximation and other software. Technical Report DITC 142/89, National Physical Laboratory, Teddington, 1989.

[dB76] C. de Boor. Splines as linear combinations of B-splines. In G. G. Lorentz, C. K. Chui, and L. L. Schumaker, editors, Approximation Theory II, pages 1–47. Academic Press, New York, 1976.

[dB78] C. de Boor. A practical guide to splines. Springer-Verlag, New York, Heidelberg, Berlin, 1978.

[dBH85] C. de Boor and K. Hoellig. B-Splines without divided differences. TSR 2896, University of Wisconsin at Madison, 1985.

[DBMS79] J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. Stewart. LIN-PACK. Users guide. SIAM Publications, 1979.

[dBP77] C. de Boor and A. Pinkus. Backward error analysis for totally positive linear systems. Numer. Math., 27:485–490, 1977.

[DH91] P. Deuflhard and A. Hohmann. Numerische Mathematik, Eine algorithmisch orientierte Einführung. Walter de Gruyter, Berlin. New York, 1991.

[DS83] J. E. Dennis and R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice Hall, Inc., 1983.

[Eld84] L. Elden. An algorithm for the regularization of ill-conditioned banded least squares problems. SIAM J. Sci. Stat. Comput., 5:237–254, 1984.

[FC80] F. N. Fritsch and R. E. Carlson. Monotone piecewise cubic interpolation. SIAM J. Numer. Anal., 17(2):238–246, April 1980.

[Fle87] R. Fletcher. Practical Methods of Optimization. John Wiley & Sons, 1987.

[Gay84] D. M. Gay. A trust-region approach to linearly constrained optimization. In D. F. Griffiths, editor, Proceedings 1983 Dundee Conference on Numerical Analysis, pages 72–105, Berlin, 1984. Lecture Notes in Mathematics 1066, Springer.

[GMW81] P. E. Gill, W. Murray, and M. H. Wright. Practical Optimization. Academic Press Inc., 1981.

[GP73] G. H. Golub and V. Pereyra. The differentiation of pseudoinverses and nonlinear least squares problems whose variables separate. SIAM J. Numer. Anal., 10:413–432, 1973.

[GP76] G. H. Golub and V. Pereyra. Differentiation of pseudoinverses, separable nonlinear least squares problems and other tales. In M. Nashed, editor, Generalized Inverses and Applications, pages 303–324. Academic Press, New York, 1976.

[GvL89] G. H. Golub and C. F. van Loan. Matrix Computations. The John Hopkins University Press, Oxford, Baltimore, 1989.

[HF79] J. N. Holt and R. Fletcher. An algorithm for constrained nonlinear least squares. J. Inst. Math. Appl., 23:449–463, 1979.

[Kau75] L. Kaufman. A variable projection method for solving separable nonlinear least squares problems. BIT, 15(4):49–57, 1975.

[KP78] L. Kaufman and V. Pereyra. A method for separable nonlinear least squares problems with separably nonlinear equality constraints. SIAM J. Numer. Anal., 15:12–20, 1978.

[Kro74] F. T. Krogh. Efficient implementation of a variable projection algorithm for nonlinear least squares problems. Comm. ACM, 17(3):167– 169, March 1974.

[KS88] A. Kiełbasinski and H. Schwetlick. Numerische Lineare Algebra. Deutscher Verlag der Wissenschaften, Berlin, 1988.

[Kun89] V. Kunert. Splineglättung mit Ungleichheitsbedingungen an Ableitungen. Master’s thesis, MLU Halle-Wittenberg, 1989.

[LH74] C. L. Lawson and R. J. Hanson. Solving Least Squares Problems. Prentice Hall, Inc., 1974.

[LM88] T. Lyche and K. Mørken. A data reduction strategy for splines with applications to the approximation of functions and data. IMA J. Numer. Anal., 8:185–208, 1988.

[Par85] T. A. Parks. Reducible Nonlinear Programming Problems. PhD thesis, Houston Univ., Dept. of Mathematics, Houston, 1985.

[RW80] A. Ruhe and P. Å. Wedin. Algorithms for separable nonlinear least squares problems. SIAM Review, 22(3):318–337, July 1980.

[Sch79] H. Schwetlick. Numerische Lösung nichtlinearer Gleichungen. Deutscher Verlag der Wissenschaften, 1979.

[Sch81] L. L. Schumaker. Spline Functions: Basic Theory. John Wiley and Sons, New York, 1981.

[SK93] H. Schwetlick and V. Kunert. Spline smoothing under constraints on derivatives. BIT, 1993. to appear.

[Spe93] P. Spellucci. Numerische Verfahren der nichtlinearen Optimierung. Birkhäuser, 1993.

[Wah90] G. Wahba. Spline Models for Observational Data. SIAM, Philadelphia, Pennsylvania, 1990.

[WH85] S. J. Wright and J. N. Holt. Algorithms for nonlinear least squares with linear inequality constraints. SIAM J. Sci. Statist. Comp., 6:1033–1048, 1985.

[WL69] H. Wold and E. Lyttkens. Nonlinear iterative partial least squares (NIPALS) estimation procedures. Bull. ISI, 43:29–51, 1969.

[Wri92] S. J. Wright. An interior point algorithm for linearly constrained optimization. SIAM J. Optimization, 1992. to appear.

Erklärung

Hiermit erkläre ich, daß ich die vorliegende Arbeit selbständig und nur unter Verwendung der angegebenen Literatur angefertigt habe.

Halle/Saale, den . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .