�

Westfälische Wilhems-Universität Münster Institut für Numerische und instrumentelle Mathematik

Iterative Verfahren zur Lösung linearer Komplementaritätsprobleme,

angewendet auf ein Marktgleichgewichtsproblem

Diplomarbeit vorgelegt von Michael Pusmentirer im Februar 1995

nach einer Aufgabenstellung von Prof.Dr.C.W.Cryer

Inhaltsverzeichnis

1 Einführung 4

2 Mathematische Grundlagen 6

2.1 Lineare Ungleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Trennung konvexer Mengen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Die Alternativtheoreme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Konvexe Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Q-Matrizen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Iterative Lösung linearer Gleichungssysteme . . . . . . . . . . . . . . . . . . . . . . . 18

2.4.1 Die Standard-Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4.2 Konvergente Matrizen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.3 Konvergenz des SOR-Verfahrens . . . . . . . . . . . . . . . . . . . . . . . . . 22

3 Iterative Verfahren zur Lösung eines LKPs 24

3.1 Basis-Splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.1 Die BS-Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.2 Das projezierte Jacobi-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.2 Ein klassischer Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.1 Das LKP (q, M) als duales quadratisches Programm . . . . . . . . . . . . . . 28 3.2.2 Das Hildreth-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2.3 Erweiterter Hildreth-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.3 Das projezierte SOR-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3.1 Konvergenzgeschwindigkeit des Verfahrens . . . . . . . . . . . . . . . . . . . . 42

4 Marktgleichgewicht (Economic Equilibrium) 44

4.1 Was ist Marktgleichgewicht ? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1.1 Ein kurzer Exkurs in die Makroökonomik . . . . . . . . . . . . . . . . . . . . 44 4.1.2 Netzwerkoptimierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.1.3 Der Baum-Wald-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Verfahren zur Netzwerkoptimierung . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.2.1 “Quadratisches Netzwerk” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2.2 Nicht-lineares Netzwerk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

5 Der Konvergenzbeweis von Luo und Tseng (1991) 57 5.1 Eigenschaften regulärer Splittings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.2 Konvergenz des Verfahrens z^{(k+1)} = F_{B}(z^{(k)}) . . . . . . . . . . . . . . . . . . . . . . 64

5.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2.2 Separation der Koordinaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5.3 Die SOR-Erweiterung des Konvergenzbeweises . . . . . . . . . . . . . . . . . . . . . 72

6 Numerik 75

6.1 Lösung eines quadratischen Netzwerkes mit dem erweiterten Hildreth-Verfahren … 75 6.1.1 Der numerische Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.1.2 Numerische Ergebnisse und Bemerkungen . . . . . . . . . . . . . . . . . . . . 77

6.2 Der PSOR-Algorithmus als Lösungsverfahren eines nicht-linearen Netzwerkproblems 79 6.2.1 Der numerische Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 6.2.2 Numerische Ergebnisse und Bemerkungen . . . . . . . . . . . . . . . . . . . . 80

6.3 Programmierhinweise zum Baum-Wald-Algorithmus . . . . . . . . . . . . . . . . . . 81

7 Schlußbemerkung 84

Kapitel 1

Einführung

Als lineares Komplementaritätsproblem bezeichnet man die folgende Aufgabe: Gegeben ist eine quadratische Matrix M ∈ IR^{n×n} sowie ein Vektor q ∈ IR^{n}, und gesucht ist ein Vektor x ∈ IR^{n}, für den gilt:

x ≥ 0  (1.1)
Mx + q ≥ 0  (1.2)
x^{t}(Mx + q) = 0  (1.3)

Die Abschätzungen (1.1) und (1.2) sind komponentenweise gemeint, d.h. {x ≥ 0 ⇔ x_{i} ≥ 0 ∀ i = 1, ..., n} und {Mx ≥ −q ⇔ (Mx)_{i} ≥ −q_{i} ∀ i = 1, ..., n}. Solche Vektorabschätzungen werden im Verlauf dieser Arbeit immer wieder auftauchen. Ein Vektor, der (1.1) und (1.2) erfüllt, heißt zulässig. Ein Vektor, der (1.3) erfüllt, heißt komplementär. Es ist klar, daß die Lösbarkeit des linearen Komplementaritätsproblems (1.1) bis (1.3) im wesentlichen von den vorgegebenen Werten in q und M abhängt. Wir bezeichnen (1.1)-(1.3) als LKP (q, M) . Um ein LKP (q, M) zu lösen, bieten sich Pivotmethoden und iterative Verfahren an. Pivotmethoden beruhen auf dem Austauschritt, der auch bei linearer Programmierung verwendet wird. Hochdimensionale LKP s (d.h. besonders große n) führen bei Pivotmethoden zu einem enormen Rechenaufwand. In diesem Fall ist es effizienter, iterative Verfahren zu benutzen. Als eigenständige mathematische Aufgabe betrachtet man das LKP (q, M) seit Mitte der 60er Jahre. Bis dahin ergab es sich nur als Umformulierung eines Quadratischen Programms, wie etwa bei [Hildreth ]. Lineare Komplementaritätsprobleme tauchen in der Volkswirtschaftslehre auf, wenn Gleichgewichtspreise bestimmt werden sollen. Eine Reihe von Autoren ( i.B. [Glassey ], [Jones ], [Matth 1 ] ) nennen diese Modelle, erklären jedoch nicht ihren ökonomischen Hintergrund. So erwähnt etwa [Cottle ], daß die Lösung eines Netzwerkproblems auf die Lösung eines LKP (q, M) mit M = ABA^{t} und q = Aa + t hinausläuft – die Matrix A ( im folgenden Beschränkungsmatrix genannt ) wird aber nicht definiert.

KAPITEL 1. EINFÜHRUNG

Wesentliches Anliegen in Kapitel 4 ist es darum, “Quadratisches Netzwerk” und nicht-lineares Netzwerk vor dem volkswirtschaftlichen Hintergrund darzustellen. Die Ansätze von [Cottle ] und [Jones ] werden aufgenommen und erläutert.

Zur Lösung eines LKP (q, M) werden die Iterationsverfahren von [Hildreth ], [Lent ] und [Cryer ] vorgestellt. Dazu gehört die vollständige Ausführung ihrer Konvergenzbeweise. Interpretiert man die 3 Algorithmen als Basis-Splitting-Methoden, dann läßt sich die Konvergenz nach einer Idee von Luo und Tseng (siehe [Luo ] ) zeigen. Dieser Beweis ist weniger restriktiv als die klassischen Ansätze, die allesamt eine streng konvexe Zielfunktion voraussetzen.

[Lent ] benutzt das erweiterte Hildreth-Verfahren zur Lösung eines speziellen LKP (q, M) (bzw. eines speziellen quadratischen Programms, was auf dasselbe hinausläuft). Die Transformation des LKP (q, M) auf diese besondere Form paßt haargenau zu den Datenvorgaben im “Quadratischen Netzwerk”. Ein solches läßt sich darum sehr elegant mit dem erweiterten Hildreth-Verfahren lösen.

Im Modell des nichtlinearen Netzwerk wird eine Folge von LKP -Lösungen berechnet. [Jones ] nennt kein besonderes Verfahren, mit dem die Folge-LKP s gelöst werden sollen. [Matth 2 ] benutzt dazu den Lemke-Algorithmus, eine Pivotmethode. Er spricht sich kritisch hinsichtlich der Anwendung iterativer Verfahren aus. Verwendet man an dieser Stelle z.B. das PSOR-Verfahren, so muß man durch die Diagonalelemente der Nullmatrix dividieren. Wir umgehen diese Schwierigkeit, indem wir anstelle der Nullmatrix eine Diagonalmatrix mit kleinen Diagonalelementen benutzen ( vergleiche Abschnitt 6.2.1 ). Mathematisch ist dieser Ansatz korrekt, er widerspricht jedoch der Modelltheorie des nichtlinearen Netzwerkes. Es bleibt Gegenstand der Diskussion, beides gegeneinander abzuwägen.

In dieser Arbeit werden z.T. sehr neue Ideen ( z.B. Kapitel 5 oder Abschnitt 4.2.2 ) angesprochen, aber auch Quellen zitiert, die wesentlich älter als 30 Jahre sind. Obligatorische Numerik ( z.B. Abschnitt 2.4 ) hat seinen Platz, ebenso wie eher ungewöhnliche Theorien ( z.B. Abschnitt 2.1.2 ). Diese “mathematischen Mosaiksteine” ergeben zusammen eine Theorie, aus der ersichtlich wird, was ein Netzwerkproblem ist, und wie man es iterativ mit verschiedenen LKP - Verfahren lösen kann.

Mathematische Grundlagen

2.1 Lineare Ungleichungen

Die Tatsache, daß ein LKP (q, M) unter bestimmten Voraussetzungen äquivalent ist zu einem konvexen quadratischen Programm, wird bei vielen Beweisen ausgenutzt. Ein konvexes quadratisches Programm besteht aus einer Zielfunktion und einem Zulässigkeitsbereich; dieser hat stets die Gestalt {z ∈ IR^{n} : Az ≥ c, A ∈ IR^{m×n}, c ∈ IR^{m}}. Es ist daher unerlässlich, die Eigenschaften dieser konvexen Menge zu untersuchen und Kriterien für die Zulässigkeit von Vektoren herzuleiten. Solche Kriterien sind die Alternativtheoreme.

2.1.1 Trennung konvexer Mengen

Definition 2.1 : Eine Menge D ⊆ IR^{n} ist konvex, wenn für alle x, y ∈ D und λ ∈ [0, 1] gilt: λx + (1 − λ)y ∈ D.

Lemma 2.2 : D ⊆ IR^{n} sei offen und konvex. Wähle δ > 0 so, daß für d ∈ D B_{δ}(d) := {x ∈ D : ∥x − d∥ < δ} in D enthalten ist. Dann ist B_{δ}(d) konvex.

Beweis: Es seien λ ∈ [0, 1] und y, z ∈ B_{δ}(d). Zu zeigen ist: λz + (1 − λ)y ∈ B_{δ}(d). y ∈ B_{δ}(d) ⇒ ∥y − d∥ < δ ⇒ ∥(1 − λ)y − (1 − λ)d∥ < (1 − λ)δ. z ∈ B_{δ}(d) ⇒ ∥z − d∥ < δ ⇒ ∥λz − λd∥ < λδ. Addiere beide Ungleichungen und wende die △-Ungleichung an: ∥λz + (1 − λ)y − d∥ ≤ ∥(1 − λ)y − (1 − λ)d∥ + ∥λz − λd∥ < δ. Das ist die Behauptung. □

Satz 2.3 : Für A ∈ IR^{n×m} und d_{i} ∈ IR^{n}, i = 1, 2 definiere man: P_{1} = {y ∈ IR^{m} : Ay ≤ d_{1}} und P_{2} = {y ∈ IR^{m} : Ay ≤ d_{2}} , P_{2} ≠ ∅. dann existiert eine Konstante α(A) , so daß es zu jedem y_{1} ∈ P_{1} ein y_{2} ∈ P_{2} gibt mit ∥y_{1} − y_{2}∥ ≤ α ∥ [d_{1} − d_{2}]^{+}∥.

KAPITEL 2. MATHEMATISCHE GRUNDLAGEN

Beweis: Dieser Satz ist dem Artikel von [Daniel ] entnommen. □ Es seien a ∈ IR^{n} \ {0} und b ∈ IR. Dann ist mit H := {x ∈ IR^{n} : a^{t}x = b} eine Hyperebene im IR^{n} definiert. Weiter seien S_{1}, S_{2} konvexe Punktmengen im IR^{n}, S_{i} ≠ ∅, i = 1, 2. H trennt S_{1} und S_{2}, wenn a^{t}x ≥ b ∀ x ∈ S_{1} und a^{t}x ≤ b ∀ x ∈ S_{2} erfüllt ist. H trennt S_{1} und S_{2} strikt, wenn a^{t}x > b ∀ x ∈ S_{1} und a^{t}x < b ∀ x ∈ S_{2} erfüllt ist. Eine S_{1} und S_{2} trennende Hyperebene H existiert, wenn S_{1} ∩ S_{2} = ∅.

Satz 2.4 : K ≠ ∅ sei eine abgeschlossene und konvexe Teilmenge des IR^{n}. Sei p ∈ IR^{n} \ K. Dann folgt die Existenz eines eindeutigen x_{p} ∈ K mit ∥x_{p} − p∥ ≤ ∥x − p∥ ∀ x ∈ K. x_{p} ist der am nächsten bei p liegende Punkt aus K. ⇔ (x_{p} − p)^{t}(x − x_{p}) ≥ 0 ∀ x ∈ K.

Beweis: Der Satz ist aus der FA als “Projektionssatz” bekannt. Ein Beweis findet sich z.B. in ( [Alt ], S.82ff. ). □ Den Vektor x_{p} bezeichnet man als orthogonale Projektion von p auf die Menge K.

Satz 2.5 : K ≠ ∅ sei eine abgeschlossene und konvexe Teilmenge des IR^{n}. Zu jedem p ∉ K existiert ein a ∈ IR^{n} und ein b ∈ IR mit a^{t}p > b ≥ a^{t}x ∀ x ∈ K.

Beweis: Nach Satz 2.4 gibt es ein eindeutiges x_{p} ∈ K mit x_{p}^{t}x − p^{t}x − x_{p}^{t}x_{p} + p^{t}x_{p} = (x_{p} − p)^{t}(x − x_{p}) ≥ 0 ∀ x ∈ K.

Wir setzen b := −x_{p}^{t}x_{p} + p^{t}x_{p} und a := p − x_{p}. Damit erhalten wir sofort die 2.Abschätzung:

b = −x_{p}^{t}x_{p} + p^{t}x_{p} ≥ −x_{p}^{t}x + p^{t}x = a^{t}x ∀ x ∈ K.

Beweis der 1.Abschätzung: a ≠ 0, da p ≠ x_{p}, d.h. 0 > −a^{t}a = −a^{t}(p − x_{p}) ⇒ a^{t}p > a^{t}x_{p} = x_{p}^{t}(p − x_{p}) = −x_{p}^{t}x_{p} + p^{t}x_{p} = b. □

Definition 2.6 : X ≠ ∅, X ⊆ IR^{n}. X ist ein Kegel, wenn tx ∈ X ∀ x ∈ X und t ≥ 0. Ist X außerdem konvex, dann ist X ein konvexer Kegel. keg A = {y ∈ IR^{m} : ∃ x ≥ 0 : y = Ax} ist der von der Matrix A ∈ IR^{m×n} erzeugte Kegel.

Ein LKP (q, M) läßt sich mit Hilfe des Kegel-Begriffes folgendermaßen ausdrücken: Sind q ∈ IR^{n} und M ∈ IR^{n×n} gegeben, so sucht man einen Vektor z ∈ IR^{2n}, der q als Element von keg (I, −M) darstellt, und für den gilt: {z_{j} = 0 ∨ z_{j+n} = 0 ∀ j ∈ {1, ..., n}}.

Korollar 2.7 : K sei abgeschlossener konvexer Kegel, p ∉ K ⇒ ∃ a ∈ IR^{n} : a^{t}p > 0 ≥ a^{t}x ∀ x ∈ K.

Beweis: Nach Satz 2.5 gibt es ein b ∈ IR und ein a ∈ IR^{n} mit a^{t}p > b ≥ a^{t}x ∀ x ∈ K. K ist ein Kegel, also folgt:

1. 0 ∈ K ⇒ b ≥ 0.

2. a^{t}x ≤ 0 ∀ x ∈ K, denn wäre dies nicht der Fall, dann gäbe es ein x̃ ∈ K : a^{t}x̃ > 0 ⇒ ∃ λ > 0 : a^{t}(λx̃) > b und λx̃ ∈ K. Dies ist ein Widerspruch.

Also: a^{t}p > b ≥ 0 ≥ a^{t}x ∀ x ∈ K. □

2.1.2 Die Alternativtheoreme

Ein Alternativtheorem besteht aus zwei Aussagen I und II. Entweder ist nur I oder nur II erfüllt. I^{C} (II^{C}) bedeutet: I (II) ist nicht erfüllt. Der Beweis eines Alternativtheorems erfolgt in zwei Schritten: Zuerst zeige, es gelten nicht beide Aussagen gleichzeitig: {I ∧ II} geht nicht. Dann zeige: Aus der Ungültigkeit von I folgt die Gültigkeit von II: {I^{C} ⇒ II}. Damit ist dann auch bewiesen: {II^{C} ⇒ I} (denn Annahme: {II^{C} ⇒ I^{C}}, dann führt {I^{C} ⇒ II} zum Widerspruch). Wir halten fest, daß für ein Alternativtheorem mit den Aussagen I und II stets gilt:

{I ⇔ II^{C}} ⇔ {I^{C} ⇔ II}  (2.1)

Satz 2.8 : (Satz von Farkas) A ∈ IR^{m×n}, b ∈ IR^{m}, b ≠ 0. Dann gilt entweder I. Ax = b und x ≥ 0 hat ein Lösung x ∈ IR^{n}. oder II. y^{t}A ≤ 0 und y^{t}b > 0 hat eine Lösung y ∈ IR^{m}.

Beweis: I und II können nicht gleichzeitig erfüllt sein, denn dann wäre 0 < y^{t}b = y^{t}(Ax) = (y^{t}A)x ≤ 0, was ein Widerspruch ist. Setze voraus, daß I nicht gilt und wende Korollar 2.7 an mit K =keg A.

I^{C} ⇒ b ∉ keg A ⟹^{Korollar 2.7} ∃ y ∈ IR^{n} : y^{t}b > 0 ≥ y^{t}z ∀ z ∈ keg A.

Alle Spalten von A sind in keg A enthalten ⇒ y^{t}A ≤ 0. □

Farkas formulierte 1902 eine Theorie der Ungleichungen; dabei ging es ihm darum, das Fouriersche Ungleichheitsprinzip der virtuellen Verschiebungen auf die analytische Mechanik anzuwenden ( Vergleiche [Farkas] ). Der Satz von Farkas ist ein wichtiges Alternativtheorem (es ist nicht das historisch älteste !), weil man auf seiner Grundlage eine Reihe weiterer Alternativtheoreme beweisen kann.

Korollar 2.9 : Es seien A ∈ IR^{m×n}, b ∈ IR^{m}, b ≠ 0. Gilt y^{t}b ≤ 0 ∀ y mit y^{t}A ≤ 0. ⇒ ∃ x ≥ 0 : Ax = b.

Beweis: Dieses Korollar ist eine Umformulierung des Satzes von Farkas. Nehme an: es gibt kein x ≥ 0 mit Ax = b.

⟹^{Satz 2.8} ∃ y : A^{t}y ≤ 0 und b^{t}y > 0. Widerspruch ! □

Satz 2.10 : Es seien A ∈ IR^{m×n}, b ∈ IR^{m} mit b ≠ 0. Dann gilt entweder: I. Ax ≥ b und x ≥ 0 hat eine Lösung oder II. y^{t}A ≤ 0 und y ≥ 0 und y^{t}b > 0 hat eine Lösung.

Beweis: I. und II. geht nicht, denn gelte beides ⇒ 0 ≤ y^{t}(Ax − b) ⇒ 0 < y^{t}b ≤ y^{t}(Ax) = (y^{t}A)x ≤ 0 ⇒ Widerspruch! Sei also I. nicht lösbar, d.h. es gibt kein x ∈ IR^{n}, x ≥ 0 und kein ξ ∈ IR^{m}, ξ ≥ 0 mit:

∑_{i=1}^{n} a_{i}x_{i} − ∑_{j=1}^{m} ξ_{j}e_{j} = b,

e_{j} ist der j-te Einheitsvektor, a_{i} ist die i-te Spalte von A. Anders formuliert ist (A, −I)(x, ξ) = b und (x, ξ) ≥ 0 nicht lösbar.

⟹^{Satz 2.8} y^{t}(A, −I) ≤ 0 und y^{t}b > 0 lösbar. ⇔ y^{t}A ≤ 0 und y^{t}b > 0 und y ≥ 0 lösbar. □

Satz 2.11 : (Satz von Ville) A ∈ IR^{m×n}. Dann gilt entweder: I. Es gibt x ∈ IR^{n} : Ax > 0 und x > 0 oder II. Es gibt y ∈ IR^{m} : y^{t}A ≤ 0 und y ≥ 0 und y ≠ 0.

Beweis: Beide Aussagen gelten nicht simultan, denn dann gälte (y^{t}A)x ≤ 0 und y^{t}(Ax) ≥ 0 und daher y^{t}Ax = 0. Dies ist ein Widerspruch zu y ≠ 0 und Ax > 0 und y ≥ 0. Aussage I sei nicht erfüllt. Sei b > 0 ⇒ es gibt kein x̃ mit Ax̃ ≥ b und x̃ ≥ 0. Denn da die Funktion z ↦ Az stetig ist, könnte man dann folgern, daß es ein λ > 0 gibt mit A(x̃ + λe) > 0 und x̃ + λe > 0, e ∈ IR^{n} hat 1 in jeder Komponente.

⟹^{Satz 2.10} y^{t}A ≤ 0 und y ≥ 0 und y^{t}b > 0 ⇒ y^{t}A ≤ 0 und y ≥ 0 und y ≠ 0.

⇒ Aussage II gilt. □

2.2 Konvexe Analysis

Daß das LKP (q, M) äquivalent zu einem quadratischen konvexen Programm ist, wurde eingangs bereits erwähnt. Auf der Grundlage von Abschnitt 2.1 läßt sich diese Eigenschaft nun genau erläutern.

Definition 2.12 : D ⊆ IR^{n} sei konvex. Eine Funktion f : D → IR heißt konvex,wenn f(λx + (1 − λ)y) ≤ λf(x) + (1 − λ)f(y) ∀ x, y ∈ D und λ ∈ [0, 1]. f heißt strikt konvex, wenn unter der Voraussetzung x ≠ y und λ ∈ (0, 1) die Ungleichung strikt erfüllt ist.

Eine Folgerung aus der mehrdimensionalen Taylorformel ist folgendes Korollar:

Korollar 2.13 : D ⊆ IR^{n} sei offen, f : D ⟶ IR sei stetig differenzierbar auf D. Dann ist f stetig auf D, für alle x ∈ D existiert ∇f(x) und für alle ξ mit x + ξ ∈ D gilt: f(x + ξ) = f(x) + [∇f(x)]^{t}ξ + φ(x, ξ)∥ξ∥, Dabei ist φ eine Funktion mit lim_{ξ⟶0} φ(x, ξ) = 0.

Beweis: Standardbeweis aus der Analysis 2. □

Satz 2.14 : (Satz über konvexe Funktionen) D ⊆ IR^{n} sei konvex und offen. f : D ⟶ IR sei konvex und stetig differenzierbar. Dann gilt: f(x) − f(y) ≥ [∇f(y)]^{t}(x − y) ∀ x, y ∈ D.

Beweis: Seien x, y ∈ D, x ≠ y. Wähle δ > 0, so daß 0 < μ < 1 und μ < \frac{δ}{∥y − x∥}. ⇒ x̂ := y + μ(x − y) = (1 − μ)y + μx ∈ B_{δ}(y) ⊂ D. Da B_{δ}(y) nach Lemma 2.2 konvex ist, ist f auf B_{δ}(y) konvex; folglich ist für 0 < λ < 1 :

(1 − λ)f(y) + λf(x̂) ≥ f((1 − λ)y + λx̂). Definiere φ wie in Korollar 2.13, dann ergibt sich

f(x̂) − f(y) ≥ \frac{f((1 − λ)y + λx̂) − f(y)}{λ} = \frac{λ[∇f(y)]^{t}(x̂ − y) + φ(y, λ(x̂ − y))λ∥y − x̂∥}{λ} =

= [∇f(y)]^{t}(x̂ − y) + φ(y, λ(x̂ − y))∥y − x̂∥.

Mit λ ⟶ 0 folgt

f(x̂) − f(y) ≥ [∇f(y)]^{t}(x̂ − y).  (2.2)

Aus x̂ ∈ D, f konvex in y und x̂ = (1 − μ)y + μx zieht man die folgenden Schlüsse:

f(x̂) − f(y) ≤ μ(f(x) − f(y))  (2.3)
x̂ − y = μ(x − y).  (2.4)

(2.2), (2.3), (2.4) ⇒

f(x̂) − f(y) ≥ [∇f(y)]^{t}μ(x − y) ≥ [∇f(y)]^{t} \frac{f(x̂) − f(y)}{f(x) − f(y)}(x − y).

⇒ f(x) − f(y) ≥ [∇f(y)]^{t}(x − y). □

Korollar 2.15 : Als direkte Folgerung aus Satz 2.14 halten wir fest: f : IR^{n} ⟶ IR sei konvex. a) Gilt für x, y ∈ IR^{n} : [∇f(x)]^{t}y > 0. ⇒ f(x + λy) > f(x) ∀ λ > 0. b) Gilt für x, y ∈ IR^{n} : [∇f(x)]^{t}y < 0. ⇒ ∃ λ_{0} > 0. mit f(x + λy) < f(x) ∀ 0 < λ < λ_{0}.

Beweis: zu a) f(x + λy) − f(x) ≥ λ[∇f(x)]^{t}y > 0. zu b) Definiere g(ξ) := y^{t}∇f(x + ξy). g(0) < 0 gilt nach Voraussetzung. g ist stetig in ξ. ⇒ ∃ λ_{0} mit g(ξ) < 0 ∀ 0 ≤ ξ < λ_{0}. Sei λ mit λ > 0 und λ < λ_{0}. ⇒ f(x) − f(x + λy) ≥ −λ[∇f(x + λy)]^{t}y = −λg(λ) > 0.

□

Lemma 2.16 : D ⊆ IR^{n} sei offen und konvex. f : D → IR, f(z) := q^{t}z + \frac{1}{2}z^{t}Mz, wobei q ∈ IR^{n} beliebig und M ∈ IR^{n×n} eine symmetrische und positiv semidefinite Matrix ist. Dann ist f konvex.

Beweis: Seien x, y ∈ D, λ ∈ [0, 1]. Es ist λz^{t}Mz ≥ λ^{2}z^{t}Mz ∀ z ∈ IR^{n}. Weiter:

λf(x) + (1 − λ)f(y) = λq^{t}x + \frac{1}{2}λx^{t}Mx + q^{t}(1 − λ)y + \frac{1}{2}y^{t}My − \frac{λ}{2}y^{t}My =

= q^{t}(λx + (1 − λ)y) + \frac{1}{2}y^{t}My + \frac{λ}{2}(x^{t}Mx − y^{t}My) =

= q^{t}(λx + (1 − λ)y) + \frac{1}{2}y^{t}My + \frac{λ}{2}(2(x − y)^{t}My + (x − y)^{t}M(x − y)) ≥^{s.o.}

≥ q^{t}(λx + (1 − λ)y) + \frac{1}{2}y^{t}My + λ(x − y)^{t}My + \frac{λ^{2}}{2} (x − y)^{t}M(x − y) =

= q^{t}(λx + (1 − λ)y) + \frac{1}{2}(λ(x − y) + y)^{t}M(λ(x − y) + y) =

= q^{t}(λx + (1 − λ)y) + \frac{1}{2}(λx + (1 − λ)y)^{t}M(λx + (1 − λ)y) = f(λx + (1 − λ)y). □

Ist die Matrix M in Lemma 2.16 positiv definit, dann ist die entsprechende Funktion f streng konvex.

Satz 2.17 : f sei wie in Lemma 2.16, allerdings sei M ∈ IR^{n×n} positiv definit. Für alle z ≥ 0 ist f dann nach unten beschränkt.

Beweis: Wähle \overline{z} ≥ 0 mit M\overline{z} ≥ −q. Ein solches \overline{z} existiert, denn angenommen, es gibt keins

⟹^{Satz 2.10} ∃ y ∈ IR^{n} : My ≤ 0 und y ≥ 0 und −y^{t}q > 0 ⇒ y^{t}My ≤ 0.

Das ist ein Widerspruch zu M positiv definit. Für z ≥ 0 gilt z^{t}M\overline{z} ≥ −z^{t}q, mit dem Satz über konvexe Funktionen folgert man:

f(z) − f(\overline{z}) ≥ (z − \overline{z})^{t}∇f(\overline{z}) = (z − \overline{z})^{t}(q + M\overline{z}) = z^{t}q + z^{t}M\overline{z} + K ≥^{s.o.} z^{t}q − z^{t}q + K = K,

mit K := −\overline{z}^{t}q − \overline{z}^{t}M\overline{z}. □

Konvexe quadratische Programmierung

Unter konvexer quadratischer Programmierung versteht man folgendes Problem: Seien A ∈ IR^{m×n}; q ∈ IR^{n}; b ∈ IR^{m}. Die Matrix M ∈ IR^{n×n} sei symmetrisch und positiv semi-definit. Dann lautet das QP-Problem:

minimiere f(z) = q^{t}z + \frac{1}{2}z^{t}Mz
unter den Nebenbedingungen: Az ≥ b  (QP)
z ≥ 0

Satz 2.18 : Ist z ∈ IR^{n} lokale Lösung des (QP), dann existiert ein y ∈ IR^{m}, so daß (z, y) ∈ IR^{n+m} die Karush-Kuhn-Tucker-Bedingungen erfüllt: \begin{aligned}u := & q + Mz − A^{t}y & ≥ 0 , z ≥ 0 , z^{t}u = 0, \\ v := & −b + Az & ≥ 0 , y ≥ 0 , y^{t}v = 0. \end{aligned} Die KKT-Bedingungen sind hinreichend dafür, daß z globale Lösung des (QP) ist.

Beweis: Die zweite Aussage läßt sich verhältnismäßig kurz beweisen: Sei z ∈ IR^{n}, so daß z die KKT-Bedingungen erfüllt. Definiere P := {y ∈ IR^{n} : Ay ≥ b, y ≥ 0}. z ∈ P ist klar. Zu zeigen: f(z) ≤ f(x) ∀ x ∈ P. Sei also x ∈ P, d.h. Ax ≥ b und x ≥ 0. Aus den KKT-Bedingungen folgt: 0 = y^{t}Az − y^{t}b ≥ y^{t}Az − y^{t}Ax ⇒ y^{t}Ax ≥ y^{t}Az. Weiterhin: x^{t}(q + Mz − A^{t}y) ≥ 0, d.h. q^{t}x + x^{t}Mz ≥ x^{t}A^{t}y ≥ y^{t}Az = q^{t}z + z^{t}Mz. Addiert man \frac{1}{2}x^{t}Mx − x^{t}Mz + \frac{1}{2}z^{t}Mz = \frac{1}{2}(x − z)^{t}M(x − z) ≥ 0, so folgt: f(x) ≥ f(z).

Der Beweis der ersten Aussage:

Definiere die verallgemeinerte Lagrange-Funktion

L(z, y) = f(z) + ∑_{j=1}^{m} y_{j}(b − Az)_{j} = f(z) + y^{t}(b − Az) = f(z) − y^{t}(Az − b).

Es gibt, so setzen wir voraus, einen zulässigen Vektor \overline{z} ∈ IR^{n} mit \overline{z} ≥ 0 und A\overline{z} > b. z^{∗} sei Lösung des (QP). Wir zeigen:

∃ y^{∗} ≥ 0 : L(z^{∗}, y) ≤ L(z^{∗}, y^{∗}) ≤ L(z, y^{∗}) ∀ z ≥ 0 ∀ y ≥ 0.  (2.5)

Die Vektoren des IR^{m+1} bezeichne mit (ξ_{0}, ..., ξ_{m})^{t}. Definiere die beiden Punktmengen

K_{1} := {ξ ∈ IR^{m+1} : ξ_{0} ≥ f(z), ξ_{i} ≥ (b − Az)_{i}, i = 1, ..., m für ein z ≥ 0}. K_{2} := {ξ ∈ IR^{m+1} : ξ_{0} < f(z^{∗}), ξ_{i} < 0, i = 1, ..., m}.

Die Funktionen f(z) und (b − Az) sind konvex, folglich sind K_{1} und K_{2} konvexe Mengen. z^{∗} ist Optimalvektor ⇒ K_{1} ∩ K_{2} = ∅. Nach den Bemerkungen aus Abschnitt 2.1.1 gibt es eine K_{1} und K_{2} trennende Hyperebene H = {x ∈ IR^{m+1} : a^{t}x = b, a ≠ 0}, für die gilt: a^{t}y_{1} ≥ a^{t}y_{2} ∀ y_{1} ∈ K_{1} ∀ y_{2} ∈ K_{2}.

Die y_{2} können beliebig klein werden im negativen Bereich, d.h. aus der Ungleichung folgt a ≥ 0. Definiere für z ≥ 0 : y_{1} = (f(z), (b − Az)_{1}, ..., (b − Az)_{m})^{t} und y_{2} = (f(z^{∗}), 0, .., 0)^{t}. Die Ungleichung ist auch erfüllt für Punkte auf dem Rand von K_{2}. ⇒ a_{0}f(z) + ∑_{j=1}^{m} a_{j}(b − Az)_{j} ≥ a_{0}f(z^{∗}) ∀ z ≥ 0. Es ist a_{0} > 0 ( Wäre dies nicht der Fall, dann gälte für jedes z ≥ 0 : ∑_{j=1}^{m} a_{j}(b − Az)_{j} ≥ 0 und a_{j} ≥ 0 ∀ j, und es gäbe ein k mit a_{k} > 0. Für das oben vorausgesetzte \overline{z} ist dies jedoch falsch! ).

Setze y_{j}^{∗} = \frac{a_{j}}{a_{0}} , j = 1, ..., m ⇒ f(z) + ∑_{j=1}^{m} y_{j}^{∗}(b − Az)_{j} ≥ f(z^{∗}) bzw. L(z, y^{∗}) ≥ f(z^{∗}).

Setze z = z^{∗} ⇒ (b − Az^{∗})^{t}y^{∗} ≥ 0. z^{∗} ist zulässig, folglich b − Az^{∗} ≤ 0 ⇒ (b − Az^{∗})^{t}y^{∗} = 0 ⇒ L(z^{∗}, y^{∗}) = f(z^{∗}).

Ferner gilt: f(z^{∗}) ≥ f(z^{∗}) + y^{t}(b − Az^{∗}) = L(z^{∗}, y) ∀ y ≥ 0.

Insgesamt kommt man auf (2.5).

Zu zeigen bleibt, daß (2.5) auf die KKT-Bedingungen führt. L ist differenzierbar in z und y. Die Abschätzungskette (2.5) ist äquivalent zu den Bedingungen

\frac{∂L}{∂z_{i}} (z^{∗}, y^{∗}) ≥ 0, z_{i}^{∗} ≥ 0, z_{i}^{∗} (\frac{∂L}{∂z_{i}} (z^{∗}, y^{∗})) = 0  (2.6)
\frac{∂L}{∂y_{j}} (z^{∗}, y^{∗}) ≤ 0, y_{j}^{∗} ≥ 0, y_{j}^{∗} (\frac{∂L}{∂y_{j}} (z^{∗}, y^{∗})) = 0  (2.7)

∀ i = 1, ..., n ∀ j = 1, ..., m. (2.6) und (2.7) sind notwendige Bedingungen, denn: setze y^{∗} ≥ 0 fest, dann ist L(z, y^{∗}) konvex in z_{i}. Gilt (2.6) für z_{i}^{∗} ≥ 0 nicht,

d.h. es ist \{z_{i}^{∗} > 0 ∧ \frac{∂L(z^{∗}, y^{∗})}{∂z_{i}} ≠ 0\} ∨ \{z_{i}^{∗} = 0 ∧ \frac{∂L(z^{∗}, y^{∗})}{∂z_{i}} < 0\}

⇒ es gibt Punkte z_{i} ≥ 0, in denen L einen geringeren Wert annimmt als in z_{i}^{∗}. Dies ist eine Folgerung aus Korollar 2.15 b). Analog führt die Ungültigkeit von (2.7) zu dem Schluß, daß es ein y_{j} ≥ 0 gibt, in dem L größer als in y_{j}^{∗} ist. (2.6) und (2.7) sind hinreichende Bedingungen, denn L(z, y^{∗}) ist konvex in z. Aus Satz 2.14 folgt

L(z, y^{∗}) ≥ L(z^{∗}, y^{∗}) + (z − z^{∗})^{t} \frac{∂L(z^{∗}, y^{∗})}{∂z} = L(z^{∗}, y^{∗}) + z^{t} \frac{∂L(z^{∗}, y^{∗})}{∂z} ≥ L(z^{∗}, y^{∗}) ∀ z ≥ 0.

L(z^{∗}, y) ist linear in y

⇒ L(z^{∗}, y) = L(z^{∗}, y^{∗}) + (y − y^{∗})^{t} \frac{∂L(z^{∗}, y^{∗})}{∂y} = L(z^{∗}, y^{∗}) + y^{t} \frac{∂L(z^{∗}, y^{∗})}{∂y} ≤ L(z^{∗}, y^{∗}) ∀ y ≥ 0.

(2.6) und (2.7) lassen sich schreiben als

\frac{∂L(z^{∗}, y^{∗})}{∂z} ≥ 0, [z^{∗}]^{t} \frac{∂L(z^{∗}, y^{∗})}{∂z} = 0, z^{∗} ≥ 0,

\frac{∂L(z^{∗}, y^{∗})}{∂y} ≤ 0, [y^{∗}]^{t} \frac{∂L(z^{∗}, y^{∗})}{∂y} = 0, y^{∗} ≥ 0.

Und das sind die KKT-Bedingungen. □

Lemma 2.19 : z ∈ IR^{n} löst das (QP) ohne die Restriktionen z ≥ 0 ⇔ es existiert ein y ∈ IR^{m} mit v := Az − b  (2.8) v ≥ 0  (2.9) q + Mz − A^{t}y = 0  (2.10) y ≥ 0  (2.11) y^{t}v = 0  (2.12)

Beweis: ′′ ⟹′′ Unter Berücksichtigung folgender Bemerkung verläuft der Beweis analog zum Beweis der ersten Aussage in Satz 2.18. L ist die verallgemeinerte Lagrange-Funktion. Fällt die Restriktion z ≥ 0 in (2.5) weg, so gilt L(z^{∗}, y^{∗}) ≤ L(z, y^{∗}) ∀ z ∈ IR^{n}. ⇔ ∇_{z}L(z^{∗}, y^{∗}) = 0. Bedingung (2.6) reduziert sich folglich auf \frac{∂L(z^{∗}, y^{∗})}{∂z_{i}} = 0. ′′ ⟸′′ Der Beweis der Rückrichtung wird geführt wie der Beweis der zweiten Aussage in Satz 2.18. □

Korollar 2.20 : Das Problem minimiere f(z) = q^{t}z + \frac{1}{2}z^{t}Mz unter den N.B. z ≥ 0 ist äquivalent zum LKP (q, M) , wenn M positiv definit und symmetrisch ist.

Beweis: Nach Satz 2.18 sind die KKT- Bedingungen notwendig und hinreichend für die Lösung des (QP)- Problems. Für das angegebene (QP)-Problem lauten die KKT-Bedingungen (d.h. A = 0 und b = 0):

q + Mz ≥ 0, z ≥ 0, z^{t}(q + Mz) = 0.

Dies ist das LKP (q, M) . □

2.3 Q-Matrizen

Das LKP (q, M) ist eindeutig lösbar für q ∈ IR^{n} beliebig, wenn die Matrix M positiv definit ist. Ist M positiv semi-definit, dann sind m.E. Aussagen über die Lösbarkeit des LKP (q, M) möglich. Wir betrachten an dieser Stelle die Matrix-Klassen Q und Q_{0}: Q_{0} ist die Klasse der Matrizen M, für die jedes zulässige LKP (q, M) lösbar ist. Q ist die Klasse der Matrizen M, für die das LKP (q, M) lösbar ist. Da die Verfahren in Kapitel 3 sich vornehmlich auf Probleme mit positiv-semidefiniten, bzw. positiv definiten Matrizen beziehen, soll gezeigt werden, daß diese Matrizen Unterklassen von Q_{0}, bzw. Q sind.

Definition 2.21 : Ein Polyhedron ist eine Menge der Form P = {x ∈ IR^{n} : a_{i}^{t}x ≥ b_{i}, a_{i} ∈ IR^{n}, b_{i} ∈ IR, i = 1, ..., n}. P ist der Durchschnitt endlich vieler abgeschlossener Halbräume im IR^{n}.

Satz 2.22 : (Frank-Wolfe) Seien M ∈ IR^{n×n}, q ∈ IR^{n}. Ist die quadratische Funktion f(z) := q^{t}z + z^{t}Mz nach unten beschränkt auf einem nicht-leeren Polyhedron P, dann nimmt f ihr Infimum auf P an.

Beweis: Dieser Satz findet sich in seiner ursprünglichen Form in [Frank ]. □

Definition 2.23 : M ∈ IR^{n×n} ist positiv semidefinit : ⇔ z^{t}Mz ≥ 0 ∀ z ∈ IR^{n}. M ∈ IR^{n×n} ist positiv definit : ⇔ z^{t}Mz > 0 ∀ z ∈ IR^{n}.

Eine positiv (semi-) definite Matrix ist in diesem Text nicht notwendig symmetrisch. Die Aufgabe

minimiere z^{t}(q + Mz)
unter den Nebenbedingungen q + Mz ≥ 0  (AQP)
z ≥ 0

bezeichnet man als alternatives quadratisches Programm. Sehr leicht einzusehen ist die Äquivalenz: {z^{∗} löst das (AQP) ∧ (q + Mz^{∗})^{t}z^{∗} = 0} ⇔ z^{∗} löst das LKP (q, M) .

Lemma 2.24 : Ist das LKP (q, M) zulässig, dann hat das (AQP) eine optimale Lösung z^{∗}. Darüberhinaus existiert ein u ∈ IR^{n} mit q + (M + M^{t})z^{∗} − M^{t}u ≥ 0  (2.13) (z^{∗})^{t}(q + (M + M^{t})z^{∗} − M^{t}u) = 0  (2.14) u ≥ 0  (2.15) u^{t}(q + Mz^{∗}) = 0  (2.16) (z^{∗} − u)_{i}(M^{t}(z^{∗} − u))_{i} ≤ 0 ∀ i = 1, ..., n  (2.17)

Beweis: Die Existenz einer optimalen Lösung z^{∗} folgt aus Satz 2.22. z^{∗} ist zulässiger Vektor des LKP (q, M) . Da z^{t}(q+Mz) = z^{t}q+ \frac{1}{2}z^{t}(M +M^{t})z sind die Bedingungen (2.13)-(2.16) die KKT-Bedingungen des (AQP).

(2.14) ⇒ \underbrace{(z^{∗})^{t}q + (z^{∗})^{t}Mz^{∗}}_{≥0} +(z^{∗})^{t}M^{t}(z^{∗} − u) = 0.

⇒ z_{i}^{∗}(M^{t}(z^{∗} − u))_{i} ≤ 0 ∀ i = 1, ..., n.

(2.13) ⇒ \underbrace{u^{t}(q + Mz^{∗})}_{0, wg.(2.16)} +u^{t}M^{t}z^{∗} − u^{t}M^{t}u ≥ 0

⇒ −u_{i}(M^{t}(z^{∗} − u))_{i} ≤ 0 ∀ i = 1, ..., n. Addiere beide Ungleichungen ⇒ (2.17). □

Satz 2.25 : Ist M ∈ IR^{n×n} positiv semidefinit und das LKP (q, M) zulässig, dann ist das LKP (q, M) lösbar.

Beweis: Das LKP (q, M) ist zulässig, aus Lemma 2.24 folgt dann die Existenz von Vektoren z^{∗} und u, so daß z^{∗} zulässiger Vektor des LKP (q, M) ist und die Bedingungen (2.13)-(2.17) gelten. (2.17) ⇒ (z^{∗} − u)^{t}M^{t}(z^{∗} − u) ≤ 0. M^{t} ist positiv semidefinit, d.h. (z^{∗} − u)^{t}M^{t}(z^{∗} − u) = 0. ⇒

\underbrace{(z^{∗})^{t}M^{t}(z^{∗} − u)}_{≤0} + \underbrace{(−u^{t}M^{t}(z^{∗} − u))}_{≤0} = 0.

Daß die beiden Summanden nicht-positiv sind, wurde im Beweis von Lemma 2.24 verwendet. Es ergibt sich (z^{∗})^{t}M^{t}(z^{∗} − u) = 0. Setze dies in (2.14) ein: (z^{∗})^{t}(q + Mz^{∗}) = 0, also löst z^{∗} das LKP (q, M) . □

Bei den sich anschließenden Sätzen greifen wir zurück auf den Abschnitt 2.1.2 Alternativtheoreme; vor allen Dingen benötigen wir den Satz 2.11 (Satz von Ville).

Definition 2.26 : Eine Matrix M ∈ IR^{n×n}, zu der ein z ∈ IR^{n} existiert mit Mz > 0 und z > 0, heißt S-Matrix.

Lemma 2.27 : Ist M ∈ IR^{n×n} positiv definit ⇒ M ist S-Matrix.

Beweis: Wir nehmen an, M ist keine S-Matrix, d.h. es gibt kein z mit Mz > 0 und z > 0

⟹^{Satz 2.11} ∃ u ∈ IR^{n} : u ≥ 0 und u ≠ 0 und M^{t}u ≤ 0

⇒ u^{t}M^{t}u ≤ 0. Dies ist ein Widerspruch zur Voraussetzung, daß M positiv definit ist. □

Korollar 2.28 : M ist S-Matrix ⇔ zu M ∈ IR^{n×n} existiert ein z ∈ IR^{n} mit Mz > 0 und z ≥ 0.

Beweis: ′′ ⟹′′ ist offensichtlich

′′ ⟸′′ Mz ist stetig in z, d.h. für λ > 0 hinreichend klein gilt: M(z + λe) > 0 und z + λe > 0. Der Vektor e ∈ IR^{n} hat 1 in jeder Komponente. □

Lemma 2.29 : M ∈ IR^{n×n} ist S-Matrix ⇔ das LKP (q, M) ist zulässig ∀ q ∈ IR^{n}.

Beweis: ′′ ⟹′′ Gegeben sei ein LKP (q, M) mit M ∈ S. z erfülle Mz > 0 und z > 0. λ sei so groß, daß: λMz = Mλz ≥ −q. λz ≥ 0 ist klar. ⇒ λz ist zulässiger Vektor. ′′ ⟸′′ Das LKP (q, M) sei zulässig für q beliebig. Wähle q < 0. Jede zulässige Lösung z erfüllt Mz ≥ −q > 0 sowie z ≥ 0. Folglich ist M S-Matrix. □

Satz 2.30 : Ist M ∈ IR^{n×n} positiv definit, dann hat das LKP (q, M) eine eindeutige Lösung für alle q ∈ IR^{n}.

Beweis: Aus Lemma 2.26 ergibt sich, daß M S-Matrix, und darum nach Lemma 2.29 zulässig ist. Aus Satz 2.25 folgt damit die Lösbarkeit des LKP (q, M) .

Die Lösung z^{∗} des LKP (q, M) minimiert die Zielfunktion f(z) = z^{t}(Mz + q) des (AQP). Die Eindeutigkeit von z^{∗} ergibt sich daraus, daß M positiv definit und f somit streng konvex ist. □

Klar ist: Q = Q_{0} ∩ S. Ob M eine S-Matrix ist, läßt sich z.B. mit linearer Programmierung feststellen. Könnte man numerisch untersuchen, ob M zur Klasse der Q_{0}-Matrizen gehört, so hätte man damit überprüft, ob M eine Q-Matrix ist. Allerdings existieren dazu noch keine effizienten Algorithmen. Aus Satz 2.30 folgt, daß alle positiv definiten Matrizen zur Klasse der Q-Matrizen gehören. Aus Satz 2.25 folgt, daß alle positiv semi-definiten Matrizen zur Klasse der Q_{0}-Matrizen gehören.

2.4 Iterative Lösung linearer Gleichungssysteme

Der Inhalt dieses Abschnitts ist Numerik-Grundwissen. Allerdings rechtfertigt der direkte Bezug der Numerik linearer Gleichungssysteme zur LKP -Numerik diese ausführliche Darstellung.

A ∈ IR^{n×n} sei regulär mit nicht-negativen Diagonalelementen. Für b ∈ IR^{n} suchen wir eine Lösung des Gleichungssystems Ax = b. Das Gleichungssystem numerisch exakt zu lösen, bietet sich nur an, wenn n relativ klein ist – sonst wird der Rechenaufwand zu groß. Die numerischen Iterationsverfahren basieren darauf, Ax = b als Fixpunktproblem zu formulieren; d.h. wir bestimmen eine Matrix B(A) ∈ IR^{n×n} und einen Vektor d ∈ IR^{n} derart, daß

Ax = b ⇔ x = B(A)x + d.

Die einzelnen Methoden haben alle die Form

x^{(k+1)} = B(A)x^{(k)} + d.  (2.18)

Sie sind formal mit den iterativen Verfahren zur Lösung des LKP (q, M) verwandt. Z.T. lassen sich die Konvergenzaussagen vom SOR- zum PSOR-Verfahren übertragen.

Definition 2.31 : Sei A ∈ IR^{n×n}. Eine Zerlegung der Form A = A_{1} + A_{2}, A_{i} ∈ IR^{n×n} heißt Matrix- Splitting und wird bezeichnet mit (A_{1}, A_{2}).

2.4.1 Die Standard-Methoden

Gesamtschrittverfahren (Jacobi-Verfahren)

Ax = b ⇔ ∑_{j=1}^{n} a_{ij}x_{j} − b_{i} = 0, i = 1, ..., n. Man löst nach x_{i} auf:

x_{i} = −\frac{1}{a_{ii}} [∑_{{j=1 \atop j≠i}}^{n} a_{ij}x_{j} + b_{i}], i = 1, ..., n.

Das Verfahren lautet demzufolge:

x_{i}^{(k+1)} = −\frac{1}{a_{ii}} [∑_{{j=1 \atop j≠i}}^{n} a_{ij}x_{j}^{(k)} + b_{i}]  (2.19)

i = 1, ..., n. Bezüglich A wird das Splitting (M, −N) durchgeführt mit M := D und N := L + R. D ist die Diagonale von A, L und R sind die linken/rechten Dreiecksmatrizen von A. Ax = b ⇒ Dx = b + (L + R)x ⇒ x^{(k+1)} = D^{−1}(L + R)x^{(k)} + D^{−1}b. Die Iterationsmatrix lautet also: B(A) = D^{−1}(L + R).

Einzelschrittverfahren (Gauss-Seidel-Verfahren)

Im Unterschied zum Jacobi-Verfahren werden die neu berechneten Komponenten zur Berechnung der folgenden Komponenten des Iterationsvektors mit verarbeitet:

x_{i}^{(k+1)} = −\frac{1}{a_{ii}} [∑_{j=1}^{i−1} a_{ij}x_{j}^{(k+1)} + ∑_{j=i+1}^{n} a_{ij}x_{j}^{(k)} − b_{i}]

i = 1, ..., n und k = 1, 2, 3, ... Verwendet wird ebenfalls das Splitting (M, −N), M := D und N := L + R wie oben, allerdings lautet das Verfahren: Dx^{(k+1)} = Lx^{(k+1)} + Rx^{(k)} + b ⇒ x^{(k+1)} = (D − L)^{−1}Rx^{(k)} + (D − R)^{−1}b. B(A) = (D − L)^{−1}R ist hier die Iterationsmatrix.

Das SOR-Verfahren

Das Gauss-Seidel-Verfahren ist ein Sonderfall des SOR-Verfahrens. Anstatt den Iterationsvektor x_{i}^{(k+1)}, i = 1, ..., n des ES-Verfahrens direkt zu verwenden, bestimmt man ein gewichtetes Mittel zwischen der alten und der neuen Iteration:

△x_{i}^{(k+1)} := x_{i}^{(k+1)} − x_{i}^{(k)} = −\frac{1}{a_{ii}} (∑_{j=1}^{i−1} a_{ij}x_{j}^{(k+1)} + ∑_{j=i}^{n} a_{ij}x_{j}^{(k)} − b_{i}) .

⇒

x_{i}^{(k+1)} = x_{i}^{(k)} + ω △ x_{i}^{(k+1)} = (1 − ω)x_{i}^{(k)} − \frac{ω}{a_{ii}} (∑_{j=1}^{i−1} a_{ij}x_{j}^{(k+1)} + ∑_{j=i+1}^{n} a_{ij}x_{j}^{(k)} − b_{i}) .

Verwendet wird dabei das Matrix-Splitting (M, −N), also A = M − N mit

M = \frac{1}{ω} D − L, N = \frac{1 − ω}{ω} D + R, 0 < ω < 2.

Als Iterationsverfahren ergibt sich:

x^{(k+1)} = B_{ω}(A)x^{(k)} + d_{ω}, B_{ω}(A) := M^{−1}N, d_{ω} := M^{−1}b

. Ausgeschrieben: B_{ω}(A) = (D − ωL)^{−1}((1 − ω)D + ωR) und d_{ω} = (\frac{1}{ω} D − L)^{−1}.

2.4.2 Konvergente Matrizen

In diesem Abschnitt soll untersucht werden, welche Eigenschaften die Matrix B(A) ∈ IR^{n×n} erfüllen muß, damit das Verfahren (2.18) konvergiert. Es wird definiert, was unter der Konvergenzgeschwindigkeit des Verfahrens zu verstehen ist. Wir kommen in Abschnitt 3.3.1 auf diese Definition zurück, wenn es um die Konvergenz des PSOR-Verfahrens geht. Der Fehler von (2.18) beträgt: e^{(k)} = x^{(k)} − x. ⇒ e^{(k+1)} = B(A)e^{(k)}. ⇒ e^{(k)} = [B(A)]^{k}e^{(0)}, k ≥ 0. Es ist klar, daß

lim_{k⟶∞} e_{j}^{(k)} = 0 ⇔ lim_{k⟶∞} x_{j}^{(k)} = x

∀ j ∈ {1, ..., n}. Gilt also

lim_{k⟶∞}[B(A)]^{k}e^{(0)} = 0 ∈ IR^{n×n},  (2.20)

dann ist (2.18) konvergent.

Definition 2.32 : A ∈ IR^{n×n} habe die Eigenwerte λ_{1}, ..., λ_{m}, m ≤ n. Dann ist ρ(A) = max{ |λ_{i}|; 1 ≤ i ≤ m } der Spektralradius von A.

Definition 2.33 : A ∈ IR^{n×n} ⇒ ∥A∥ = sup_{x≠0} \frac{∥Ax∥}{∥x∥} wird als Spektralnorm von A bezeichnet.

Lemma 2.34 : Für A, B ∈ IR^{n×n} gelten die Aussagen 1. ∥A∥ > 0, wenn A ≠ 0 ∈ IR^{n×n} und ∥αA∥ = |α| ∥A∥ ∀ α ∈ IR. 2. ∥A + B∥ ≤ ∥A∥ + ∥B∥. ∥AB∥ ≤ ∥A∥ ∥B∥. 3. ∥Ax∥ ≤ ∥A∥ ∥x∥ ∀ x ∈ IR^{n}.

Beweis: Aussage 1. ergibt sich sofort aus den Eigenschaften von Vektornormen. zu Aussage 2.:

μ, η, ξ ∈ IR^{n} \ \{0\} seien Vektoren mit ∥A∥ = \frac{∥Aξ∥}{∥ξ∥}, ∥B∥ = \frac{∥Bη∥}{∥η∥} und ∥A + B∥ = \frac{∥(A + B)μ∥}{∥μ∥}.

Aus der △-Ungleichung für Vektornormen erhält man:

∥A∥ + ∥B∥ = \frac{∥Aξ∥}{∥ξ∥} + \frac{∥Bη∥}{∥η∥} ≥ \frac{1}{∥μ∥}(∥Aμ∥ + ∥Bμ∥) ≥ \frac{∥(A + B)μ∥}{∥μ∥} = ∥A + B∥.

Des weiteren:

∥A∥ ≥ \frac{∥ABξ∥}{∥Bξ∥} = \frac{∥ABξ∥}{∥Bξ∥} \frac{∥ξ∥}{∥ξ∥} = \frac{∥AB∥}{∥B∥}.

Aussage 3 ist klar für x = 0. Ist x ≠ 0, dann folgt ∥A∥ ≥ \frac{∥Ax∥}{∥x∥}. □

Korollar 2.35 : A ∈ IR^{n×n} ⇒ ∥A∥ ≥ ρ(A).

Beweis: Sei λ Eigenwert von A, v der dazugehörige Eigenvektor. Dann erhält man: |λ| ∥v∥ = ∥λv∥ = ∥Av∥ ≤ ∥A∥ ∥v∥. ⇒ |λ| ≤ ∥A∥ für jeden Eigenwert von A. □

Satz 2.36 : A ∈ IR^{n×n} ⇒ ∥A∥ = \sqrt{ρ(A^{t}A)}.

Beweis: A^{t}A ∈ IR^{n×n} ist symmetrisch und positiv semidefinit, denn x^{t}A^{t}Ax = ∥Ax∥^{2} ≥ 0 ∀ x ∈ IR^{n}. Folglich hat A^{t}A eine Orthonormalbasis von Eigenvektoren v_{1}, ..., v_{n} und nicht-negative Eigenwerte λ_{1}, ..., λ_{n}. Œ sei 0 ≤ λ_{1} ≤ ... ≤ λ_{n}.

Sei x ∈ IR^{n}, dann existieren ξ_{1}, ..., ξ_{n} mit x = ∑_{i=1}^{n} ξ_{i}v_{i}. ⇒

(\frac{∥Ax∥}{∥x∥})^{2} = \frac{x^{t}A^{t}Ax}{x^{t}x} = \frac{∑_{i=1}^{n} ξ_{i}^{2}λ_{i}}{∑_{i=1}^{n} ξ_{i}^{2}} ⇒ 0 ≤ λ_{1} ≤ (\frac{∥Ax∥}{∥x∥})^{2} ≤ λ_{n}.

Setzt man x = v_{n}, dann gilt die letzte Abschätzung als Gleichung und darum:

∥A∥^{2} = sup_{x≠0} (\frac{∥Ax∥}{∥x∥})^{2} = λ_{n} = ρ(A^{t}A). □

Korollar 2.37 : Ist A ∈ IR^{n×n} symmetrisch, dann folgt ∥A∥ = ρ(A).

Beweis:

A^{t} = A ⟹^{Satz 2.36} ∥A∥ = \sqrt{ρ(A^{t}A)} = \sqrt{ρ(A^{2})} = λ_{n} = ρ(A).

λ_{n} ist der betragsgrößte Eigenwert von A. □

Definition 2.38 : A ∈ IR^{n×n}. A heißt konvergent (gegen 0 ∈ IR^{n×n}), wenn A^{k}_{\overrightarrow{k⟶∞}}0.

Wegen Definition 2.38 und (2.20) ist das Verfahren (2.18) konvergent, wenn die Iterationsmatrix es ist.

Satz 2.39 : A ∈ C^{n×n}. A ist konvergent ⇔ ρ(A) < 1.

Beweis: Stelle von A die Jordansche Normalform her, d.h. es gibt eine reguläre Matrix S ∈ IR^{n×n} :

SAS^{−1} =: A_{J} = ( \begin{matrix}J_{1} & & & & & 0 \\ & J_{2} \\ & & . & . & . \\ 0 & & & & & J_{r} \end{matrix} ) mit J_{l} ∈ IR^{n_{l}×n_{l}} dergestalt J_{l} = ( \begin{matrix}λ_{l} & 1 & & 0 \\ & ... & ... \\ & & ... & 1 \\ 0 & & & λ_{l} \end{matrix} ) ,

1 ≤ l ≤ r und n = ∑_{l=1}^{r} n_{l}.

A_{J} ist obere trianguläre △-Matrix, die λ_{l} sind paarweise verschiedene Eigenwerte ∀ l. A und A_{J} sind ähnlich. Für m ≥ 1 erhält man:

A_{J}^{m} = ( \begin{matrix}J_{1}^{m} & & & & & 0 \\ & J_{2}^{m} \\ & & . & . & . \\ 0 & & & & & J_{r}^{m} \end{matrix} ).

Durch die besondere Struktur der J_{l} lassen sich die J_{l}^{m} auf folgende Weise berechnen:

(J_{l}^{m})_{i,j} = \{ \begin{matrix}0, & j < i \\ \binom{m}{j−i}λ_{l}^{m−j+i}, & i ≤ j ≤ min\{n_{l}, m + i\} \\ 0, & m + i < j ≤ n_{l} \end{matrix}  (2.21)

\binom{m}{k} = \frac{m!}{k!(m − k)!} ist der Binominalkoeffizient.

Nun zum eigentlichen Beweis: ′′ ⟹′′ Ist A konvergent ⇒ A^{k}_{\overrightarrow{k⟶∞}}0 ⇒ (A_{J})^{k}_{\overrightarrow{k⟶∞}}0, da A_{J}^{k} = SA^{k}S^{−1}.

⇒ (J_{l})^{k}_{\overrightarrow{k⟶∞}}0 ⇒ |λ_{l}| < 1 ∀ 1 ≤ l ≤ r.

Somit ist klar, daß ρ(A) = ρ(A_{J}) = max_{1≤l≤r} |λ_{l}| < 1.

′′ ⟸′′ Es sei ρ(A) = ρ(A_{J}) < 1 ⇒ |λ_{l}| < 1 ∀ 1 ≤ l ≤ r.

Mit (2.21) erhält man: lim_{k⟶∞}(J_{l}^{k})_{i,j} = 0 ∀ 1 ≤ i, j ≤ n_{l}.

Jede Matrix J_{l} ist also konvergent. Nach Definition ist A_{J} ebenfalls konvergent. Daß A konvergent ist, folgt schließlich aus A^{k} = S^{−1}A_{J}^{k}S. □

2.4.3 Konvergenz des SOR-Verfahrens

Das Verfahren (2.18) konvergiert, wenn ρ(B(A)) < 1. Es ist bekannt, daß das SOR-Verfahren nur für 0 < ω < 2 konvergiert ( siehe z.B. [Schwarz ] ).

Wie schnell das SOR-Verfahren in Abhängigkeit von ω konvergiert, soll i.f. skizziert werden.

Definition 2.40 : Es seien A, B ∈ IR^{n×n}, m > 0 und ∥A^{m}∥ < 1. Unter der durchschnittlichen Konvergenzgeschwindigkeit von (2.15) mit der Ite- rationsmatrix A nach m Schritten versteht man: R(A^{m}) = − ln[(∥A^{m}∥^{\frac{1}{m}}] = \frac{− ln ∥A^{m}∥}{m}. R(A^{m}) < R(B^{m}) bedeutet, daß bei m Iterationsschritten das Verfahren x^{(k+1)} = Bx^{(k)} + d schneller konvergiert als x^{(k+1)} = Ax^{(k)} + d.

Nutzt man die Eigenschaften von Vektor-, bzw. Matrizennormen aus, dann ergibt sich: ∥e^{(k)}∥ ≤ ∥[B(A)∥^{k}∥e^{(0)}∥, k ≥ 0.

Der durchschnittliche Reduzierungsfaktor beträgt R = (\frac{∥e^{(k)}∥}{∥e^{(0)}∥})^{\frac{1}{k}}.

Für ∥A^{k}∥ < 1 folgt dann: R ≤ (∥A^{k}∥)^{\frac{1}{k}} = e^{−R(A^{k})}.

R(A^{k}) ist demnach die exponentielle Abstiegsrate, mit der sich ρ(A) einer unteren Schranke nähert.

Satz 2.41 : A ∈ IR^{n×n} sei eine konvergente Matrix ⇒ lim_{k⟶∞} R(A^{k}) = − ln ρ(A) =: R_{∞}(A).

Beweis: Sei A ∈ IR^{n×n} mit ρ(A) > 0. J_{r} seien die Jordan-Kästchen der Jordanschen Normalform von A mit ρ(J_{r}) = ρ(A). p sei deren größte Ordnung. ν > 0 konstant. Dann läßt sich zeigen:

∥A^{k}∥ ∼ ν\binom{k}{p − 1}[ρ(A)]^{k−(p−1)}, für k ⟶ ∞.

Ein ausführlicher Beweis ist zu finden in ( [Varga ], S. 65 ff). Damit folgt sofort:

R(A^{k}) = \frac{− ln ∥A^{k}∥}{k} ∼ \frac{− ln ν}{k} − \frac{ln \binom{k}{p−1}}{k} − (\frac{k − p + 1}{k}) ln ρ(A).

Da \frac{ln \binom{k}{p−1}}{k}_{\overrightarrow{k⟶∞}}0 ergibt sich die Behauptung. □

R_{∞}(A) heißt asymptotische Konvergenzgeschwindigkeit. Äquivalent läßt sich auch schreiben:

R_{∞}(A) = − ln\{sup_{e^{(0)}} lim sup_{k⟶∞} ∥e^{(k)}∥^{\frac{1}{k}} \}  (2.22)

Man führt Relaxationsparameter bei Iterationsverfahren ein, um R_{∞}(A) zu verbessern. ω_{opt} ist der Relaxationsparameter, für den das Verfahren am schnellsten konvergiert. Ist B_{ω}(A) die Iterationsmatrix, dann ist klar: R_{∞}(B_{ω_{opt}}(A)) = max_{0<ω<2} R_{∞}(B_{ω}(A)).

Zur Bestimmung von ω_{opt} siehe [Varga ].

Iterative Verfahren zur Lösung eines LKPs

3.1 Basis-Splitting

Basis-Splitting-Methoden sind Iterationsverfahren zur Lösung eines LKP (q, M) . Bei der BS-Methode bestimmt man eine Vektorenfolge {z^{(k)}}_{k} ⊆ IR^{n}, so daß z^{(k+1)} das LKP (q + Cz^{(k)}, B) löst. Dabei wird M gesplittet in M = B + C. Unter bestimmten Voraussetzungen löst jeder Häufungspunkt dieser Folge das LKP (q, M) . Die BS-Methode

z^{(k+1)} = [z^{(k+1)} − (Bz^{(k+1)} + Cz^{(k)} + q)]^{+}

konvergiert, wenn die Matrizen B − C und M positiv semi-definit sind. Der Beweis dieser Aussage gelang Luo und Tseng 1991 – ihr Ansatz wird ausgeführt in Kapitel 5.

3.1.1 Die BS-Methode

M ∈ IR^{n×n} sei symmetrisch und positiv semidefinit. Das LKP (q, M) , q ∈ IR^{n} soll gelöst werden. Dazu wird zunächst das Matrix-Splitting (B, C) bzgl. M durchgeführt, d.h. M = B + C. Definiere q_{z} = q + Cz ∀ z ∈ IR^{n}. L sei Lösungsmenge des LKP (q_{z}, B). F sei eine Funktion, die einem Punkt des IR^{n} eine Teilmenge des IR^{n} zuordnet. ^{1}

Ist das LKP (q_{z}, B) eindeutig lösbar und gilt F (z) = L, dann ist klar:

F (z) = z ⇔ z löst das LKP (q, M)  (3.1)

Auf dieser Grundlage formuliert man die

1In der Literatur, so z.B. bei [Cottle ], wird diese Funktion als eine multivalued function bezeichnet.

KAPITEL 3. ITERATIVE VERFAHREN ZUR LÖSUNG EINES LKPS

Basis-Splitting-Methode (BS-Methode)

Schritt 0:

Initialisierung

Wähle z^{(0)} ≥ 0 beliebig, k = 0.

Schritt 1:

Iterationsschritt z^{(k)} ≥ 0 sei gegeben Löse das LKP (q_{k}, B) mit q_{k} = q + Cz^{(k)}. Bezeichne die Lösung als z^{(k+1)}.

Schritt 2: \underline{Abbruchtest}

Gilt für ε > 0 \{ \begin{matrix}||z^{(k)} − z^{(k+1)}|| < ε : & Abbruch \\ sonst : & gehe zu Schritt 1 \end{matrix}

Die Folge {z^{(k)}}_{k}, die durch die BS-Methode entsteht, nennen wir BS-Folge. Die Effektivität der BS-Methode hängt ab von dem Matrix-Splitting (B, C), insbesondere von B. B muß so gewählt sein, daß das Verfahren z^{(k+1)} = F (z^{(k)}) für alle z^{(0)} ≥ 0 konvergiert. Dies setzt voraus:

1. B ist Q-Matrix ⇔ L ≠ ∅.

2. F (z) ∈ L, wobei die Definition von F Konvergenzuntersuchungen möglich macht.

Definition 3.1 : M ∈ IR^{n×n} wird gesplittet in (B, C). Ist B eine Q-Matrix, dann ist (B, C) ein Q-Splitting.

Definition 3.2 : Ist (B − C) positiv semidefinit, dann heißt das Splitting (B, C) schwach regulär; ist (B − C) positiv definit, dann wird es als regulär bezeichnet.

In den sich anschließenden Sätzen dieses Abschnitts ist f(z) = q^{t}z + \frac{1}{2}z^{t}Mz, mit M ∈ IR^{n×n} wie jeweils bezeichnet.

Lemma 3.3 : M ∈ IR^{n×n} sei symmetrisch, (B, C) sei ein schwach reguläres Q-Splitting von M. Die Folge \{z^{(k)}\}_{k} entstehe durch die BS-Methode. ⇒ f(z^{(k)}) − f(z^{(k+1)}) ≥ \frac{1}{2}(z^{(k)} − z^{(k+1)})^{t}(B − C)(z^{(k)} − z^{(k+1)}) ≥ 0. Ist (B, C) regulär, so gilt: f(z^{(k)}) = f(z^{(k+1)}) ⇔ z^{(k)} = z^{(k+1)}.

Beweis: f(z^{(k)}) − f(z^{(k+1)}) = (z^{(k)} − z^{(k+1)})^{t}(q + Mz^{(k+1)}) + \frac{1}{2}(z^{(k)} − z^{(k+1)})^{t}M(z^{(k)} − z^{(k+1)}) = = (z^{(k)} − z^{(k+1)})^{t}(q + Cz^{(k)} + Bz^{(k+1)}) + \frac{1}{2}(z^{(k)} − z^{(k+1)})^{t}(B − C)(z^{(k)} − z^{(k+1)}) ≥ ≥ \frac{1}{2}(z^{(k)} − z^{(k+1)})^{t}(B − C)(z^{(k)} − z^{(k+1)}).

Die letzte Ungleichung ergibt sich dadurch, daß z^{(k)} ≥ 0 ist und z^{(k+1)} das LKP (q_{k}, B) löst: q + Cz^{(k)} + Bz^{(k+1)} ≥ 0 (q + Cz^{(k)} + Bz^{(k+1)})^{t}z^{(k+1)} = 0. □

Ist z^{(k)} = z^{(k+1)}, dann löst z^{(k)} das LKP (q, M) . Also können wir Œ im weiteren annehmen, daß z^{(k)} ≠ z^{(k+1)} für die Vektoren der BS-Methode gilt und {f(z^{(k)})}_{k} monoton fallend ist.

Satz 3.4 : (B, C) sei reguläres Q-Splitting der symmetrischen Matrix M. Dann löst jeder Häufungspunkt einer BS-Folge {z^{(k)}}_{k}, das LKP (q, M) .

Beweis: z^{HP} sei Häufungspunkt der BS-Folge. \{z^{(k_{i})}\}_{i} sei eine Teilfolge mit z^{(k_{i})}_{\overrightarrow{i⟶∞}}z^{HP}. Daraus folgt, daß f(z^{(k_{i})})_{\overrightarrow{i⟶∞}}f(z^{HP}). Da die Folge monoton fallend ist, ist sie nach unten beschränkt.

Folglich konvergiert {f(z^{(k)})}_{k}. Mit Lemma 3.3 schlußfolgert man aus der Regularität des Splittings (B, C): \{z^{(k)} − z^{(k+1)}\}_{k} konvergiert gegen 0 und daher z^{(k_{i}+1)}_{\overrightarrow{i⟶∞}}z^{HP}. Nach Definition erfüllt z^{(k_{i}+1)} die Bedingungen:

q + Cz^{(k_{i})} + Bz^{(k_{i}+1)} ≥ 0
z^{(k_{i}+1)} ≥ 0
(q + Cz^{(k_{i})} + Bz^{(k_{i}+1)})^{t}z^{(k_{i}+1)} = 0

Für k_{i} ⟶ ∞ folgt, daß z^{HP} das LKP (q, M) löst. □

3.1.2 Das projezierte Jacobi-Verfahren

Man programmiert den Iterationsschritt mit Hilfe einer Funktion F, für die gilt: F (z^{(k)}) löst das LKP (q_{k}, B). Die Bezeichnungen entsprechen denen im Algorithmus. Mit z^{(k+1)} = F (z^{(k)}) folgt aus der Abbruchbedingung, daß z^{(k)} mit der Genauigkeit ε ein Fixpunkt von F ist. Nach (3.1) ist z^{(k)} dann eine Lösung des LKP (q, M) . Indem bestimmte Bedingungen an die Matrix B im Splitting (B, C) gestellt werden, läßt sich eine Funktion F herleiten, die im Iterationsschritt der BS-Methode verwendet werden kann.

Definition 3.5 : Eine invertierbare Matrix B ∈ IR^{n×n} heißt inversstabil, wenn gilt: Bx ≥ y ⇔ x ≥ B^{−1}y ∀ x, y ∈ IR^{n}.

Satz 3.6 : Sei q_{x} := q + Cx ∀ x, B = M − C wie oben eingeführt. B sei inversstabil. Definiere F : IR^{n} ⟶ IR^{n} als F (x) := max{0, −B^{−1}q_{x}}. Dann gilt: F (z) = L, wobei L die Lösungsmenge des LKP (q_{z}, B) ist.

Beweis:

F (z) ⊆ L  F (z) ist zulässig und komplementär: (a) F (z) ≥ 0 (b) F (z) ≥ −B^{−1}q_{z} ⇒ q_{z} + BF (z) ≥ 0 (c) In jeder Komponente gilt entweder (a) oder (b) als Gleichung Aus (a), (b) und (c) folgt: F (z) ∈ L. .

F (z) ⊇ L Sei x ∈ L, dann gilt:

(a) x ≥ 0

(b) x ≥ −B^{−1}q_{z}

(c) In jeder Komponente gilt entweder (a) oder (b) als Gleichung.

Es ist klar, daß { x_{i} = 0 ⇒ (−B^{−1}q_{z})_{i} ≤ 0 } und { x_{i} = (−B^{−1}q_{z})_{i} ⇒ (−B^{−1}q_{z})_{i} ≥ 0 }. Insofern gilt x = max{0; −B^{−1}q_{z}} = F (z). □

Mit Satz 3.36 läßt sich folgende BS-Methode zur iterativen Lösung des LKP (q, M) herleiten:

\begin{aligned}z^{(0)} & ≥ & 0 \\ w^{(k+1)} & := & −B^{−1}(q + Cz^{(k)}) \\ z^{(k+1)} & := & max\{0, w^{(k+1)}\} \end{aligned}

Dabei ist (B, C) ein Q-Splitting von M, und B ist inversstabil.

Ein solches Splitting liegt z.B. dann vor, wenn B Diagonalmatrix mit positiven Diagonalelementen ist. Ist B außerdem die Diagonale von M, so ist w^{(k+1)} der Iterationsvektor des Jacobi-Verfahrens zur Lösung von q + Mz = 0. Man nennt die BS-Methode in diesem Fall projeziertes Jacobi-Verfahren, da w^{(k+1)} auf den nichtnegativen Orthanten des IR^{n} projeziert wird.

Dazu ein Beispiel:

Definiere

M := ( \begin{matrix}2 & 2 & 1 \\ 2 & 3 & 1 \\ 1 & 1 & 1 \end{matrix} ) , q := ( \begin{matrix}−1 \\ 2 \\ −2 \end{matrix} ) , B = ( \begin{matrix}2 \\ & 3 \\ & & 1 \end{matrix} ) , C = ( \begin{matrix}0 & 2 & 1 \\ 2 & 0 & 1 \\ 1 & 1 & 0 \end{matrix} )

Es soll eine Lösung des LKP (q, M) bestimmt werden, wobei bzgl. M das Splitting (B, C) durchgeführt wird. Das projezierte Jacobi-Verfahren lautet in diesem Fall:

z^{ν+1} = max\{0, z^{ν} − B^{−1}(q + Mz^{ν})\} mit B^{−1}(q + Mz^{ν}) = ( \begin{matrix}−\frac{1}{2} + z_{1}^{ν} + z_{2}^{ν} + \frac{1}{2}z_{3}^{ν} \\ \frac{2}{3} + \frac{2}{3}z_{1}^{ν} + z_{2}^{ν} + \frac{1}{3}z_{3}^{ν} \\ −2 + z_{1}^{ν} + z_{2}^{ν} + z_{3}^{ν} \end{matrix} ) ⇔

⇔  z_{1}^{ν+1} = max(0, \frac{1}{2} − z_{2}^{ν} − \frac{1}{2}z_{3}^{ν}) z_{2}^{ν+1} = max(0, −\frac{2}{3} − \frac{2}{3}z_{1}^{ν} − \frac{1}{3}z_{3}^{ν}) z_{3}^{ν+1} = max(0, 2 − z_{1}^{ν} − z_{2}^{ν})

Dieses Verfahren läßt sich leicht programmieren und führt auf den Lösungsvektor (0, 0, 2)^{t}.

Das Listing findet sich im Anhang: PROJAC.C.

3.2 Ein klassischer Ansatz

Das in diesem Abschnitt dargestellte Hildreth-Verfahren stammt aus dem Jahr 1957. Hildreth’s Ziel war es, eine quadratische Funktion unter linearen Restriktionen zu maximieren. Er gab einen Algorithmus zur Lösung des dualen quadratischen Programms an, das unter den gegebenen Voraussetzungen äquivalent zum LKP (q, M) ist. Das Hildreth-Verfahren ist eine BS-Methode und kann ähnlich wie das projezierte Jacobi-Verfahren als projeziertes Gauss-Seidel-Verfahren interpretiert werden.

3.2.1 Das LKP (q, M) als duales quadratisches Programm

Es seien C ∈ IR^{n×n} symmetrisch und negativ definit (d.h. −C ist positiv definit), b ∈ IR^{n}, d ∈ IR^{m}, G ∈ IR^{m×n}. Ausgangspunkt für dieses Kapitel ist folgende Optimierungsaufgabe:

maximiere g(x) = b^{t}x + \frac{1}{2}x^{t}Cx
unter den N.B. Gx ≥ d.

(PQP)

Das (PQP) ist äquivalent zu: minimiere −g(x) unter den N.B. Gx ≥ d.

Definiere M := −GC^{−1}G^{t} und q := −d − GC^{−1}b. Es können zwei Fälle eintreten:

1. G hat vollen Zeilenrang ⇒ M pos. definit (denn G^{t}x ≠ 0 ∀ x ≠ 0 ⇒ x^{t}Mx > 0 ∀ x ≠ 0).

2. G hat keinen vollen Zeilenrang ⇒ M pos. semidefinit mit m_{ii} > 0, denn sei g_{i}^{t} i-te Spalte von G^{t} ⇒ m_{ii} = −g_{i}^{t}C^{−1}g_{i} > 0.

Mit diesen Bezeichnungen betrachte man das Problem

minimiere f(z) = q^{t}z + \frac{1}{2}z^{t}Mz
unter den N.B. z ≥ 0.

(DQP)

(DQP) heißt das zu (P QP) “duale Programm”. Dementsprechend ist (P QP) das zu (DQP) “primäre Programm”.

Satz 3.7 :  1. x^{max} ist Lösung des (P QP) ⇒ Es gibt eine Lösung z^{min} des (DQP) mit x^{max} = −C^{−1}(G^{t}z^{min} + b). 2. z^{min} sei Lösung des (DQP), x^{max} := −C^{−1}(G^{t}z^{min} + b) ⇒ x^{max} löst das (P QP).

Beweis: \underline{Aussage 1.}

x^{max} sei Lösung des (P QP), dann existiert ein z^{min} ∈ IR^{m}, so daß die KKT-Bedingungen erfüllt sind (vergleiche Lemma 2.19):

v := Gx^{max} − d  (3.2)
v ≥ 0  (3.3)
−b − Cx^{max} − G^{t}z^{min} = 0  (3.4)
z^{min} ≥ 0  (3.5)
v^{t}z^{min} = 0  (3.6)

C ist negativ definit ⇒ C^{−1} existiert und (3.4) läßt sich nach x^{max} auflösen:

x^{max} = −C^{−1}(G^{t}z^{min} + b).  (3.7)

Diese Substitution führt zu den Bedingungen:

Mz^{min} − v = −q  (3.8)
v ≥ 0  (3.9)
z^{min} ≥ 0  (3.10)
v^{t}z^{min} = 0  (3.11)

Also: z^{min} löst das LKP (q, M) .

⟹^{Korollar 2.20} z^{min} ist Lösung des (DQP).

\underline{Aussage 2.}

z^{min} löst das (DQP) ⟹^{Korollar 2.20} z^{min} löst das LKP (q, M) .

Setze v := Mz^{min} + q, dann sind (3.8)-(3.11) erfüllt. Setzt man weiterhin x^{max} := −C^{−1}(G^{t}z^{min} + b), dann gelten (3.2)-(3.6) ⟹^{Lemma 2.19} x^{max} löst das (P QP). □

3.2.2 Das Hildreth-Verfahren

Das (DQP) soll iterativ gelöst werden. Mittels Satz 3.7 erhält man daraus eine Lösung des (PQP). Wähle also einen Startvektor z^{(0)} ≥ 0. Zweckmäßig ist z^{(0)} = 0, da im Optimalvektor stets einige Komponenten 0 sind. z^{(k+1)} wird aus z^{(k)} komponentenweise berechnet, jeder Iterationsschritt k besteht also aus n Einzelschritten.

Für i = 1, ..., n setze \overline{z}_{j} := \{ \begin{matrix}z_{j}^{(k+1)} & , für j < i \\ z_{j}^{(k)} & , für j > i \end{matrix} konstant.

z_{i}^{(k+1)} bestimmt sich als Wert von z_{i} in \frac{∂f}{∂z_{i}} = 0 ⇔ q_{i} + ∑_{{j=1 \atop j≠i}}^{m} m_{ij}\overline{z}_{j} + m_{ii}z_{i} = 0.

Ist z_{i} < 0, dann setze z_{i}^{(k+1)} = 0.

Das ganze lautet als Algorithmus formuliert:

k = 0, 1, 2, 3, ...

i = 1, 2, 3, ..., m

\begin{aligned}w_{i}^{(k+1)} & = & −\frac{1}{m_{ii}} (∑_{j=1}^{i−1} m_{ij}z_{j}^{(k+1)} + ∑_{j=i+1}^{m} m_{ij}z_{j}^{(k)} + q_{i}) \\ z_{i}^{(k+1)} & = & max(0, w_{i}^{(k+1)}). \end{aligned}  (3.12)

w_{i}^{(k+1)} ist wohldefiniert, da m_{ii} ≠ 0. Über (3.12) wird eine stetige Funktion F : IR^{m} ⟶ IR^{m} definiert mit z^{(k+1)} = F (z^{(k)}). Daß F stetig ist, folgt aus der einfachen Eigenschaft: H, G : IR^{n} ⟶ IR stetig ⇒ max{H, G} stetig.

Lemma 3.8 : Das (DQP) hat eine eindeutige Lösung z^{min}, wenn M positiv definit ist.

Beweis: f(z) ist nach unten beschränkt (nach Satz 2.16), und die Menge der zulässigen Werte ist ein nichtleeres Polyhedron. Nach Satz 2.22 nimmt f ihr Infimum an, es existiert also eine Lösung.

Zur Eindeutigkeit: Annahme: z^{min} und ẑ^{min} sind Lösungen mit z^{min} ≠ ẑ^{min}. Setze u := αz^{min} + (1 − α)ẑ^{min}, dann folgt f(u) < αf(z^{min}) + (1 − α)f(ẑ^{min}) = f(z^{min}).

Dies ist ein Widerspruch dazu, daß z^{min} Minimalwert von f ist. □ Im Originalartikel [Hildreth ] finden sich Konvergenzbeweise für den positiv (semi-)definiten Fall. Das Hildreth-Verfahrens ist ein Sonderfall des PSOR-Algorithmus, dessen Konvergenz wir im weiteren Verlauf dieses Kapitels beweisen.

3.2.3 Erweiterter Hildreth-Algorithmus

Die Aufgabe lautet weiterhin:

maximiere g(z) = b^{t}z + \frac{1}{2}z^{t}Cz
unter den N.B. Gz ≥ d.

(PQP)

wie in 3.2.1 definiert.

Neben dem dualen Programm (DQP) gibt es noch eine weitere äquivalente Optimierungsaufgabe: Von der Matrix −C bilde die Cholesky-Zerlegung −C = D^{t}D. Setze z := D^{−1}x − C^{−1}b. Mit den Bezeichnungen A := −GD^{−1} und c := −d − GC^{−1}b bringt man das (PQP) auf die Form:

minimiere \frac{1}{2}∥x∥^{2}
unter den N.B. Ax ≤ c.

(3.13)

Diese Transformation ist leicht einzusehen: Zum einen gilt: Ax ≤ c ⇔ −GD^{−1}x ≤ −d − GC^{−1}b ⇔ Gz = G(D^{−1}x − C^{−1}b) ≥ d. Und zum anderen: −g(z) = −\frac{1}{2}z^{t}Cz − b^{t}z = −\frac{1}{2}z^{t}C(D^{−1}x − C^{−1}b) − b^{t}(D^{−1}x − C^{−1}b) =

= −\frac{1}{2}z^{t}CD^{−1}x + \frac{1}{2}z^{t}b − b^{t}D^{−1}x + b^{t}C^{−1}b =

= −\frac{1}{2}x^{t}(D^{−1})^{t}CD^{−1}x + \frac{1}{2}b^{t}D^{−1}x + \frac{1}{2}b^{t}D^{−1}x − \frac{1}{2}b^{t}C^{−1}b − b^{t}D^{−1}x + b^{t}C^{−1}b = \frac{1}{2}x^{t}x + \frac{1}{2}b^{t}C^{−1}b.

Die Erweiterung des Hildreth-Verfahrens wurde 1980 von Lent und Censor vorgeschlagen. Aufgrund der besonderen Struktur dieser Optimierungsaufgabe ist eine geometrische Interpretation möglich. Der Algorithmus ist durch die Verwendung fast-zyklischer Folgen weniger restriktiv als das klassische Verfahren, und außerdem läßt sich die Konvergenz durch Relaxation verbessern. Wir verwenden den erweiterten Algorithmus zur Lösung eines “Quadratischen Netzwerkes” (Abschnitt 4.2.1). Durch passende Zuweisung bringen wir dort das LKP (q, M) auf die Form (3.13). Der erweiterte Algorithmus bezieht sich nur auf (3.13) und nur dafür wird nachfolgend die Konvergenz gezeigt. Im Hildreth-Algorithmus (3.12) werden pro Iterationsschritt (Index k) die Zeilen der Matrix M verarbeitet (Index i). Letztere Folge ist die Kontrollfolge.

Definition 3.9 : Sei I := {1, ..., m} endliche Indexmenge. Eine Folge {i_{k}}_{k} heißt fast-zyklisch in I, falls (i) i_{k} ∈ I ∀ k ≥ 0 und (ii) es ein N ∈ IN gibt (die “fast-zyklische Konstante”) , so daß für alle k ≥ 0 gilt: I ⊆ {i_{k+1}, ..., i_{k+N}}. Eine fast-zyklische Folge in I heißt zyklisch, wenn N = m.

Bei der Programmierung des Verfahrens (in Kapitel 6) benutzen wir nur zyklische Folgen; pro Iterationsschritt werden somit alle Komponenten ,,abgearbeitet“. Bei anderen Anwendungen kann es numerisch sinnvoll sein, auf fast-zyklische Folgen zurückzugreifen. Als Beispiel für eine solche Anwendung läßt sich das Problem des kleinsten Abstands ^{2} nennen.

Nun zum Verfahren selbst.

Das erweiterte Hildreth-Verfahren

Schritt 0:

Initialisierung

z^{(0)} ≥ 0 sei beliebig, x^{(0)} := −A^{t}z^{(0)}.

Schritt 1: r_{k} := min\{z_{i_{k}}^{(k)}, ω_{k} \frac{c_{i_{k}} − ∑_{j=1}^{m} a_{i_{k}j}x_{j}^{(k)}}{∥a_{i_{k}}∥^{2}} \}

x^{(k+1)} = x^{(k)} + r_{k}a_{i_{k}}, z^{(k+1)} = z^{(k)} − r_{k}e_{i_{k}}

e_{i} ist der i-te Einheitsvektor.

{ω_{k}}_{k} ist die positive Folge von Relaxationsparametern.

Schritt 2:

Abbruchbedingung

Ist die Folge {i_{k}}_{k} fast-zyklisch in I, dann bricht das Verfahren ab.

Für das erweiterte Hildreth-Verfahren gelten die Annahmen: a_{i} ≠ 0 ∀ i = 1, ..., n (a_{i} i-te Zeile von A) und S ≠ ∅, mit S = {x : Ax ≤ c}.

Satz 3.10 : Ist {i_{k}}_{k} zyklische Folge und ω_{k} = 1 für alle k, dann ist das erweiterte Hildreth- Verfahren der bereits bekannte Hildreth-Algorithmus (3.12).

Beweis: Nach Definition ist a_{i} = A^{t}e_{i}. Es gilt x^{(k)} = −A^{t}z^{(k)} ∀ k, was mittels Induktion aus x^{(k+1)} = x^{(k)} + r_{k}a_{i_{k}} = −A^{t}z^{(k)} + r_{k}A^{t}e_{i_{k}} = −A^{t}z^{(k+1)} folgt. k sei gegeben, wir setzen i = i_{k}:

\begin{matrix}z^{(k+1)} = z^{(k)} − r_{k}e_{i} & = & z^{(k)} − [min\{z_{i}^{(k)}, \frac{c_{i} − ⟨a_{i}, x^{(k)}⟩}{∥a_{i}∥^{2}}\}]e_{i} = \\ = & z^{(k)} & −[max\{0, z_{i}^{(k)} − \frac{c_{i} − ⟨a_{i}, x^{(k)}⟩}{∥a_{i}∥^{2}}\} − z_{i}^{(k)}]e_{i}. \end{matrix}  (3.14)

Nun ist

2Das Problem des kleinsten Abstand lautet:

minimiere \frac{1}{2} (x − b)^{t}(x − b)
unter den N.B. Mx ≥ 0

Seien b ∈ IR^{m} und M ∈ IR^{n×m}. Definiere den Kegel K := {y ∈ IR^{m} : My ≥ 0}. Gesucht ist ein x ∈ K, das bzgl. der euklidischen Norm den geringsten Abstand zu b hat. D.h.:

\begin{aligned}−⟨a_{i}, x^{(k)}⟩ = ⟨a_{i}, A^{t}z^{(k)}⟩ & = & ⟨a_{i}, ∑_{j=1}^{m} a_{j}z_{j}^{(k)}⟩ = \\ = ∑_{j=1}^{m}⟨a_{i}, a_{j}⟩z_{j}^{(k)} & = & ∑_{j=1}^{m} m_{ij}z_{j}^{(k)}. \end{aligned}  (3.15)

Setze M := AA^{t} ∈ IR^{m×m} = (m_{ij})_{i,j}, dabei ist m_{ij} = ⟨a_{i}, a_{j}⟩ und m_{ii} = ⟨a_{i}, a_{i}⟩ ≠ 0 ∀ i ∈ I. Setzt man (3.15) in (3.14) ein:

z_{j}^{(k+1)} = z_{j}^{(k)} für j ≠ i
z_{j}^{(k+1)} = max\{0, z_{i}^{(k)} − \frac{1}{m_{ii}} (c_{i} + ∑_{j=1}^{m} m_{ij}z_{j}^{(k)})\} für j = i (3.16)

Es ist M = AA^{t} = GD^{−1}(GD^{−1})^{t} = GD^{−1}(D^{−1})^{t}G^{t} = −GC^{−1}G^{t}. Also ist M die in (3.12) verarbeitete Matrix und (3.16) eben das klassische Hildreth-Verfahren. □

Das erweiterte Hildreth-Verfahren kann geometrisch so interpretiert werden: Setze wieder ω_{k} = 1. Nach Definition ist z^{(k)} ≥ 0 ∀ k. Betrachte den k-ten Schritt, in dem x^{(k+1)} aus x^{(k)} berechnet wird. Die Kontrollfolge sei zyklisch, also i = i_{k} ∀ k. H_{i} = {x ∈ IR^{m} : a_{i}^{t}x ≤ c_{i}} sei der Halbraum, der äquivalent zur i-ten Restriktion in (3.13) ist. Unterscheide 2 Möglichkeiten:

1. x^{(k)} ∉ H_{i} ⇒ c_{i} − a_{i}^{t}x^{(k)} < 0 ∧ r_{k} = \frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}}

⇒ x^{(k+1)} ist die orthogonale Projektion von x^{(k)} auf die Hyperebene ∂H_{i}:

Abb. 1 : x^{(k+1)} für den Fall x^{(k)} ∉ H_{i}

[figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure]

H_{i}

∂H_{i}

x^{(k+1)}

x^{(k)}

[figure]

Rechnerisch zeigt man das schnell:

x^{(k+1)} = x^{(k)} + \underbrace{\frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}}}_{r_{k}} a_{i}.

(i) \underline{x^{(k+1)} ∈ ∂H_{i}}

c_{i} − a_{i}^{t}x^{(k+1)} = c_{i} − (a_{i}^{t}x^{(k)} + a_{i}^{t}a_{i} \frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}} ) = 0.

(ii) x^{(k+1)} − x^{(k)} ist orthogonal zu ∂H_{i} Es sei ξ ∈ ∂H_{i}, d.h. c_{i} − a_{i}^{t}ξ = 0 ⇔ a_{i}^{t}ξ = c_{i}. Zu zeigen: ⟨ξ − x^{(k+1)}, x^{(k+1)} − x^{(k)}⟩ = 0. Beweis:

⟨ξ − x^{(k+1)}, x^{(k+1)} − x^{(k)}⟩ = ⟨ξ − x^{(k)} − \frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}} a_{i}, \frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}} a_{i}⟩ =

⟨ξ − x^{(k)}, \frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}} a_{i}⟩ − \frac{(c_{i} − a_{i}^{t}x^{(k)})^{2}}{∥a_{i}∥^{2}} = \frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}} (c_{i} − a_{i}^{t}x^{(k)}) − \frac{(c_{i} − a_{i}^{t}x^{(k)})^{2}}{∥a_{i}∥^{2}} = 0.

2. x^{(k)} ∈ H_{i}. Auch hier nähert die Approximation sich orthogonal ∂H_{i} und zwar um den Betrag z^{(k)}a_{i}, allerdings nicht über ∂H_{i} hinaus:

Abb. 2 : x^{(k+1)} für den Fall x^{(k)} ∈ H_{i}
[figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure]

H_{i}

∂H_{i}

x^{(k+1)}
x^{(k)}

Beide Aussagen sind trivial: ∥x^{(k+1)} − x^{(k)}∥ = ∥z^{(k)}a_{i}∥, sowie für ξ, μ ∈ ∂H_{i} :

⟨ξ − μ, x^{(k+1)} − x^{(k)}⟩ = z_{i}^{(k)} \underbrace{⟨a_{i}, ξ − μ⟩}_{0} = 0.

Im Hildreth-Algorithmus gilt also immer x^{(k+1)} ∈ H_{i}, gleichzeitig wird z_{i}^{(k+1)} so klein wie möglich gemacht. Das zu (3.13) duale Programm lautet:

maximiere ϕ(z) = −\frac{1}{2}z^{t}AA^{t}z − c^{t}z
unter den Nebenbedingungen: z ≥ 0

(3.17)

Dies ist klar nach Abschnitt 3.2.1. \frac{1}{2}x^{t}x ist durch 0 nach unten beschränkt für alle x ∈ S. Nach Satz 2.22 nimmt \frac{1}{2}x^{t}x ihr Minimum auf S an, und aus Satz 3.7 folgert man, daß es ein z^{max} gibt mit ϕ(z^{max}) = max_{z≥0} ϕ(z).

Lemma 3.11 : Gilt 0 < ω_{k} ≤ 2 − ε, ε > 0, so erfüllen die beim erweiterten Hildreth-Verfahren berechneten Folgen \{z^{(k)}\}_{k}, \{x^{(k)}\}_{k} und \{r_{k}\}_{k} die Eigenschaften: (i) z^{(k)} ≥ 0 ∀ k. (ii) ϕ(z^{(k+1)}) ≥ ϕ(z^{(k)}). (iii) lim_{k⟶∞}[ϕ(z^{(k+1)}) − ϕ(z^{(k)})] = 0. (iv) r_{k}_{\overrightarrow{k⟶∞}}0. (v) (x^{(k+1)} − x^{(k)})_{\overrightarrow{k⟶∞}}0 ∧ (z^{(k+1)} − z^{(k)})_{\overrightarrow{k⟶∞}}0.

Beweis: Wir setzen wieder i = i_{k}, die Kontrollfolge sei also zyklisch. \underline{zu (i):} Per Induktion: z^{(0)} ≥ 0 gilt nach Voraussetzung, r_{k} ≤ z_{i}^{(k)} (Schritt 1) ⇒ z_{i}^{(k+1)} = z_{i}^{(k)} − r_{k} ≥ 0, die restlichen Komponenten bleiben konstant. \underline{zu (ii):} Man setzt x^{(k)} = −A^{t}z^{(k)} ∀ k ein:

ϕ(z^{(k+1)}) − ϕ(z^{(k)}) = −\frac{1}{2}(z^{(k+1)})^{t}AA^{t}z^{(k+1)} − c^{t}z^{(k+1)} + \frac{1}{2}(z^{(k+1)})^{t}AA^{t}z^{(k)} + c^{t}z^{(k)} =

= −r_{k}a_{i}^{t}x^{(k)} − \frac{1}{2}r_{k}^{2}∥a_{i}∥^{2} + c^{t}r_{k}e_{i} = ∥a_{i}∥^{2}r_{k} [\frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}} − \frac{1}{2}r_{k}].

Nach Schritt 1 gilt: r_{k} ≤ ω_{k} \frac{c_{i} − a_{i}^{t}x^{(k)}}{∥a_{i}∥^{2}} .

Die Ungleichung ist genau dann strikt erfüllt, wenn r_{k} = z_{i}^{(k)} ⇒

ϕ(z^{(k+1)}) − ϕ(z^{(k)}) ≥ ∥a_{i}∥^{2}(r_{k})^{2}(\frac{1}{ω_{k}} − \frac{1}{2}) ≥ 0,  (3.18)

denn 0 < ω_{k} < 2. \underline{zu (iii):} nach (ii) ist ϕ monoton steigend und wegen (3.17) durch max_{z≥0} ϕ(z) nach oben beschränkt. Folglich konvergiert ϕ(z^{(k)}).

\underline{zu (iv):} Sei ∥a_{i}∥^{2} ≥ a > 0 ∀ i. Aus (3.18) folgt dann:

ϕ(z^{(k+1)}) − ϕ(z^{(k)}) ≥ ar_{k}^{2} (\frac{1}{ω_{k}} − \frac{1}{2}) ≥ ar_{k}^{2}(\frac{1}{2 − ε} − \frac{1}{2}) = ar_{k}^{2} \frac{ε}{4 − 2ε} > \frac{1}{4}r_{k}^{2}ε.

mit (iii) ergibt sich r_{k}_{\overrightarrow{k⟶∞}}0. \underline{zu (v):} folgt sofort aus Schritt 1 und (iv). □ Es soll eine Folge von Mengen {S^{(k)}}_{k} konstruiert werden. Dazu zunächst

Definition 3.12 : Eine Folge von Schlupfvektoren \{q^{(k)}\}_{k} sei folgendermaßen gegeben: q^{(0)} = 0 q_{j}^{(k+1)} = \{ \begin{matrix}q_{j}^{(k)} & , wenn j ≠ i_{k} \\ −\frac{r_{k}∥a_{i}∥^{2}}{ω_{k}} + c_{i} − a_{i}^{t}x^{(k)} & , wenn j = i_{k} = i \end{matrix}

i_{k} ist der Index der im k-ten Schritt berechneten Variablen. Weiter sei:

c^{(k)} = q^{(k)} + Ax^{(k)} ∀ k.  (3.19)

Definition 3.13 : Die Folge der Zulässigkeitsmengen \{S^{(k)}\}_{k} ist definiert als: S^{(k)} = \{x : Ax ≤ c^{(k)}\}. \{S^{(k)}\}_{k} konvergiert gegen S, wenn c^{(k)}_{\overrightarrow{k⟶∞}}c.

Nach Konstruktion ist S^{(k)} ≠ ∅ für alle k.

Lemma 3.14 : \{i_{k}\}_{k} sei fast-zyklisch in I = \{1, 2, ..., m\}. Es existieren ε_{1}, ε_{2} > 0 mit ε_{1} ≤ ω_{k} ≤ 2 − ε_{2} ∀ k ⇒ c^{(k)}_{\overrightarrow{k⟶∞}}c.

Beweis: t ∈ I sei fest. Für k gegeben bezeichne l = l(k) den größten Index mit i_{l} = t und l ≤ k. {i_{k}}_{k} ist fast zyklisch ⇒ k − N ≤ l, N ist die fast-zyklische Konstante. (Vorausgesetzt ist natürlich k > N).

Abb. 3 : Zur Erläuterung des Index i_{l}

t

i_{l} i_{k} Wegen Definition 3.12 gilt: q_{t}^{(k)} = q_{t}^{(k−1)} = ... = q_{t}^{(l+1)}. Setze dies ein in (3.19) und benutze Definition 3.19:

c_{t}^{(k)} = q_{t}^{(k)} + a_{i}^{t}x^{(k)} = q_{t}^{(l+1)} + a_{i}^{t}x^{(k)} = −\frac{r_{k}∥a_{t}∥^{2}}{ω_{k}} + c_{t} + a_{i}^{t}(x^{(k)} − x^{(l)}).

Mit k ⟶ ∞ gilt auch l ⟶ ∞ ⟹^{Lemma 3.11} −\frac{r_{l}∥a_{t}∥^{2}}{ω_{k}} ⟶ 0.

Ebenfalls wegen Lemma 3.11: a_{i}^{t}(x^{(k)} − x^{(l)}) − ∑_{j=1}^{k−1} a_{i}^{t}(x^{(j+1)} − x^{(j)})_{\overrightarrow{j⟶∞}}0.

Somit bleibt: c_{t}^{(k)}_{\overrightarrow{k⟶∞}}c_{t} ⇒ c^{(k)} ⟶ c. □

Der sich anschließende Satz 3.15 ist nicht nur ein Hilfsmittel für den Beweis von Satz 3.16; er liefert darüberhinaus eine Aussage über jeden Iterationsvektor, was z.B. dann von Nutzen sein kann, wenn der Algorithmus nicht konvergiert.

Satz 3.15 : Ist x^{(k)} k-ter Iterationsvektor des erweiterten Hildreth-Verfahrens, dann minimiert x^{(k)} die Funktion \frac{1}{2}∥x∥^{2} unter der Nebenbedingung x ∈ S^{(k)}.

Beweis: x^{(k)} minimiert die verallgemeinerte Lagrange-Funktion L(x, z^{(k)}) (siehe zur Definition den Beweis zu Satz 2.18, nähere Erläuterungen finden sich in [Everett ]), also: f(x^{(k)}) + ⟨Ax^{(k)} − c, z^{(k)}⟩ ≤ f(x) + ⟨Ax − c, z^{(k)}⟩, mit f(x) = \frac{1}{2}∥x∥^{2} ∀ x ∈ IR^{n}. ⇒ f(x^{(k)}) ≤ f(x) + ⟨Ax, z^{(k)}⟩ − ⟨Ax^{(k)}, z^{(k)}⟩. Für k = 0, 1, 2, 3, ... und j ∈ I := {1, ..., m} zeigt man per Induktion: q_{j}z_{j}^{(k)} = 0. q_{j} sind die Schlupfvektoren aus Definition 3.12. Daraus folgt: f(x^{(k)}) ≤ f(x) + ⟨Ax, z^{(k)}⟩ − ⟨Ax^{(k)} + q^{(k)}, z^{(k)}⟩ = f(x) + ⟨Ax − c^{(k)}, z^{(k)}⟩ ∀ x ∈ IR^{n}. Ist x ∈ S^{(k)} ⇒ Ax − c^{(k)} ≤ 0, da z^{(k)} ≥ 0 ∀ k ⇒ f(x^{(k)}) ≤ f(x) ∀ x ∈ S^{(k)}. □

Nach Satz 2.3 ist bekannt, daß für T = {y : Ay ≤ d} und T̃ = {y : Ay ≤ d̃} eine Konstante α existiert, so daß es zu jedem y ∈ T ein ỹ ∈ T̃ gibt mit: ∥y − ỹ∥ ≤ α∥ [d − d̃]^{+}∥. Auf diese Eigenschaft greifen wir im Beweis zu Satz 3.16 zurück.

Es seien x^{min} = min_{x∈S} ∥x∥, x̂^{(k)} ∈ S^{(k)} sei ein Element mit ∥x̂^{(k)} − x^{(k)}∥ = min_{x∈S^{(k)}} ∥x − x^{(k)}∥

und x_{min}^{(k)} ∈ S^{(k)} sei ein Element mit ∥x_{min}^{(k)} − x^{min}∥ = min_{x∈S^{(k)}} ∥x − x^{min}∥.

Satz 3.16 : Es gelten die Voraussetzungen: (a) S = \{x : Ax ≤ c\} ≠ ∅ a_{i} ≠ 0 ∀ i ∈ I = \{1, ..., m\}. (b) \{i_{k}\}_{k} ist fast-zyklisch in I. (c) Es existieren ε_{1}, ε_{2} > 0 mit ε_{1} ≤ ω_{k} ≤ 2 − ε_{2} ∀ k. Daraus können wir die Schlußfolgerungen ziehen: 1. ∥x^{(k)}∥_{\overrightarrow{k⟶∞}}∥x^{min}∥. 2. ∥x̂^{(k)} − x^{(k)}∥_{\overrightarrow{k⟶∞}}0. 3. ∥x̂^{(k)}∥_{\overrightarrow{k⟶∞}}∥x^{min}∥.

Beweis: \underline{zu 1.:} Schreibe

∥x^{min}∥ ≤ ∥x̂^{(k)}∥ ≤ ∥x̂^{(k)} − x^{(k)}∥ + ∥x^{(k)}∥  (3.20)

Als Folgerung aus Satz 2.3, erhält man mit T = S^{(k)}, T̃ = S, x^{(k)} = y und α_{1} passend:

∥x^{(k)} − x̂^{(k)}∥ ≤ α_{1}∥ [c^{(k)} − c]^{+}∥.  (3.21)

Zu ε > 0 existiert nach Lemma 3.14 eine Konstante K_{1}, so daß ∥x^{min}∥ ≤ ε + ∥x^{(k)}∥ ∀ k ≥ K_{1}. Umgekehrt gilt: ∥x^{(k)}∥ ≤ ∥x_{min}^{(k)}∥ ≤ ∥x_{min}^{(k)} − x^{min}∥ + ∥x^{min}∥. Die erste Ungleichung ist eine Folgerung aus Satz 3.15. Wie oben ergibt sich mit T = S, T̃ = S^{(k)}, x^{min} = y und α_{2} konstant passend:

∥x_{min}^{(k)} − x^{min}∥ ≤ α_{2}∥ [c − c^{(k)}]^{+}∥ ⟹^{Lemma 3.14} ∃ K_{2} : ∥x^{(k)}∥ ≤ ε + ∥x^{min}∥ ∀ k ≥ K_{2}.

Damit ist 1. bewiesen. \underline{zu 2.:} Wende Lemma 3.14 auf (3.21) an. \underline{zu 3.:} Folgt aus 1., 2. und (3.20). □

Nach diesen aufwendigen Vorbereitungen liefert uns Satz 3.17 eine Aussage über die Konvergenz des erweiterten Hildreth-Verfahrens.

Satz 3.17 : Sind (a) − (c) aus Satz 3.16 gegeben, dann konvergiert die Folge {x^{(k)}}_{k} des erwei- terten Hildreth-Verfahrens gegen eine Lösung von (3.13).

Beweis: Zu zeigen ist lim_{k⟶∞} x^{(k)} = x^{min}.

Schreibe ∥x^{min} − x^{(k)}∥ ≤ ∥x^{min} − x̂^{(k)}∥ + ∥x^{(k)} − x̂^{(k)}∥. Der zweite Summand konvergiert nach Satz 3.16 2. mit k ⟶ ∞ gegen 0. Aus Satz 3.16 3. ergibt sich für die Folge {x̂^{(k)}}_{k} ⊂ S:

∥x̂^{(k)}∥ ⟶ ∥x^{min}∥ mit k ⟶ ∞.

S ist abgeschlossen und konvex.

x^{min} ist eindeutig bestimmt, also x̂^{(k)} ⟶ x^{min}, w.z.b.w. □

3.3 Das projezierte SOR-Verfahren

Im LKP (q, M) sei M ∈ IR^{n×n} positiv definit und symmetrisch, q ∈ IR^{n}. Als Lösungsverfahren benutzen wir die BS-Methode:

{z^{(0)} ≥ 0 \atop z^{(k+1)} = max\{0; w^{(k+1)}\}}

In diesem Abschnittes ist w^{(k+1)} der Iterationsvektor des SOR-Verfahrens. Aus diesem Grund bezeichnet man die Methode als projeziertes SOR-Verfahren. Zur Formulierung des Verfahrens und für den Beweis von Satz 3.18 geben wir einleitend eine zum LKP (q, M) äquivalente Problemstellung an:

Suche z, y ∈ IR^{n} mit :
Mz − y = −q
z^{t}y = 0
y ≥ 0
z ≥ 0.

(3.22)

Klar ist: {z löst das LKP (q, M) ⇔ z und y := Mz + q löst (3.22)}. Im PSOR-Algorithmus werden die Vektorenfolgen {z^{(k)}}_{k} und {y^{(k)}}_{k} gebildet, die gegen ein Vektorenpaar z, y konvergieren, welches (3.22) löst.

PSOR-Algorithmus

Schritt 0:

Initialisierung

Wähle einen Startvektor z^{(0)} ≥ 0

und den Relaxationsparameter ω : 0 < ω < 2.

Schritt 1: Bilde die Folgen von Vektoren z^{(k)}, r^{(k)}, y^{(k)} : k = 1, 2, ...

i = 1, ..., n
r_{i}^{(k+1)} = −q_{i} − ∑_{j=1}^{i−1} m_{ij}z_{j}^{(k+1)} − ∑_{j=i}^{n} m_{ij}z_{j}^{(k)}  (3.23)

z_{i}^{(k+1)} = max\{0; z_{i}^{(k)} + ω \frac{r_{i}^{(k+1)}}{m_{ii}} \}  (3.24)
y_{i}^{(k+1)} = −r_{i}^{(k+1)} + m_{ii}(z_{i}^{(k+1)} − z_{i}^{(k)}).  (3.25)

Satz 3.18 : Für die Folgen \{z^{(k)}\}_{k} und \{y^{(k)}\}_{k}, die durch den PSOR-Algorithmus berechnet werden, gilt: z^{(k)}_{\overrightarrow{k⟶∞}}z^{∗} und y^{(k)}_{\overrightarrow{k⟶∞}}q + Mz^{∗}, z^{∗} ist eine Lösung des LKP (q, M) .

Beweis: Definiere die Vektoren z^{(k+1,l)} ∈ IR^{n} über

z_{i}^{(k+1,l)} = \{ \begin{matrix}z_{i}^{(k+1)}, & für 1 ≤ i ≤ l \\ z_{i}^{(k)}, & für l < i ≤ n. \end{matrix}

Zur Vereinfachung setze z^{l} := z^{(k+1,l)} ∀ l ∈ {1, ..., n}. Des weiteren definiere

ω_{k+1,i} := \{ \begin{matrix}\frac{[z_{i}^{(k+1)} − z_{i}^{(k)}]m_{ii}}{r_{i}^{(k+1)}} & , für r_{i}^{(k+1)} ≠ 0 \\ ω & , sonst. \end{matrix}

Mit diesen Bezeichnungen ergibt sich im PSOR-Algorithmus:

r_{i}^{(k+1)} = −q − Mz^{i−1} ∧ z_{i}^{(k+1)} = z_{i}^{(k)} + ω_{k+1,i} \frac{r_{i}^{(k+1)}}{m_{ii}}.

Unter den gegebenen Voraussetzungen ist f(z) = q^{t}z + \frac{1}{2}z^{t}Mz streng konvex und nach Satz 2.17 nach unten beschränkt. f(z^{i}) − f(z^{i−1}) = (z^{i} − z^{i−1})^{t}(q + Mz^{i−1}) + \frac{1}{2}(z^{i} − z^{i−1})^{t}M(z^{i} − z^{i−1}) = = \frac{1}{2}m_{ii}(z_{i}^{i} − z_{i}^{i−1})^{2} − (z_{i}^{i} − z_{i}^{i−1})r_{i}^{(k+1)} = \frac{1}{2}m_{ii}(z_{i}^{(k+1)} − z_{i}^{(k)})^{2} − (z_{i}^{(k+1)} − z_{i}^{(k)})r_{i}^{(k+1)} =

= \frac{1}{2m_{ii}} ω_{k+1,i}^{2}[r_{i}^{(k+1)}]^{2} − \frac{1}{m_{ii}} ω_{k+1,i}[r_{i}^{(k+1)}]^{2} = \frac{[r_{i}^{(k+1)}]^{2}}{m_{ii}} ω_{k+1,i}(\frac{ω_{k+1,i}}{2} − 1) ≤ 0.

Die letzte Ungleichung ist erfüllt wegen 0 ≤ ω_{k+1,i} ≤ ω < 2. Die Folge {f(z^{(k)})}_{k} ist monoton fallend, und da sie nach unten beschränkt ist, konvergiert sie: f(z^{(k)}) ⟶ f_{min}. Wir beweisen als nächstes die Abschätzung:

|z_{i}^{(k+1)} − z_{i}^{(k)}| ≤ \sqrt{\frac{f(z^{i−1}) − f(z^{i})}{m(−\frac{1}{2} + \frac{1}{ω})}},  (3.26)

m := min_{i} m_{ii}. Für z_{i}^{(k+1)} = z_{i}^{(k)} ist die Ungleichung trivialerweise richtig, setze also Œ z_{i}^{(k+1)} ≠ z_{i}^{(k)}:

f(z^{i−1}) − f(z^{i}) =^{s.o.} \frac{[r_{i}^{(k+1)}]^{2}}{m_{ii}} (ω_{k+1,i} − \frac{ω_{k+1,i}^{2}}{2} ) = \frac{(z_{i}^{(k+1)} − z_{i}^{(k)})^{2}m_{ii}}{ω_{k+1,i}^{2}} (ω_{k+1,i} − \frac{ω_{k+1,i}^{2}}{2} ) =

= (z_{i}^{(k+1)} − z_{i}^{(k)})^{2}m_{ii}(\frac{1}{ω_{k+1,i}} − \frac{1}{2}) ≥ (z_{i}^{(k+1)} − z_{i}^{(k)})^{2}m(\frac{1}{ω} − \frac{1}{2}) ⇒ (3.26).

Aus f(z^{(k)})_{\overrightarrow{k⟶∞}}f_{min} ergibt sich deshalb (z_{i}^{(k+1)} − z_{i}^{(k)})_{\overrightarrow{k⟶∞}}0.

Die Folge \{z^{(k)}\}_{k} hat also einen Häufungspunkt z^{HP}. \{z^{(k_{i})}\}_{i} sei eine Teilfolge mit z^{(k_{i})}_{\overrightarrow{i⟶∞}}z^{HP}. Aus dem PSOR-Algorithmus folgt: r^{(k_{i})}_{\overrightarrow{i⟶∞}}r = −q − Mz^{HP} und y^{(k_{i})}_{\overrightarrow{i⟶∞}} − r. Außerdem gilt:

• z^{HP} ≥ 0.

• r ≤ 0.

Beweis: Annahme: r > 0. Dann gibt es ein ε > 0 sowie l, m ∈ IN mit r_{l}^{(k_{i})} ≥ ε ∀ k_{i} ≥ m. ⇒ z_{l}^{(k_{i})} = z_{l}^{(k_{i}−1)} + ω \frac{r_{l}^{(k_{i})}}{m_{ll}} ≥ z_{l}^{(k_{i}−1)} + \frac{εω}{m_{ll}}. Dies ist ein Widerspruch zu z_{i}^{(k+1)} − z_{i}^{(k)}_{\overrightarrow{k⟶∞}}0 ∀ i.

• r^{t}z^{HP} = 0.

Beweis: Annahme: r^{t}z^{HP} ≠ 0. ⇒ ∃ ε > 0 und l, m ∈ IN : r_{l}^{(k_{i})} ≤ −ε und z_{l}^{(k_{i})} ≥ ε ∀ k_{i} ≥ m. ⇒ 0 < ε ≤ z_{l}^{(k_{i})} = z_{l}^{(k_{i}−1)} + ω \frac{r_{l}^{(k_{i})}}{m_{ll}} . ⇒ |z_{l}^{(k_{i})} − z_{l}^{(k_{i}−1)}| ≥ ω \frac{ωε}{m_{ll}} ∀ k_{i} ≥ m. Dies ist ein Widerspruch zu z_{i}^{(k+1)} − z_{i}^{(k)}_{\overrightarrow{k⟶∞}}0 ∀ i.

Nun ist z^{(k)} ∈ X = {z ∈ IR^{n} : f(z) ≤ f(z^{(0)})} ∀ k. X ist kompakt und f ist streng konvex. Demnach ist z^{HP} der einzige Häufungspunkt von {z^{(k)}}_{k} in X und eindeutige Lösung des LKP (q, M) :

y^{(k)}_{\overrightarrow{k⟶∞}}q + Mz^{HP} ≥ 0

z^{(k)}_{\overrightarrow{k⟶∞}}z^{HP} ≥ 0

(q + Mz^{HP})^{t}z^{HP} = 0. □

Setzt man ω = 1 im PSOR-Algorithmus, dann ergibt sich das Hildreth-Verfahren. Für den Fall, daß M positiv definit ist, ist damit also auch die Konvergenz des Hildreth-Verfahrens bewiesen.

3.3.1 Konvergenzgeschwindigkeit des Verfahrens

Das PSOR-Verfahren konvergiert gegen eine Lösung des LKP (q, M) . Die Konvergenzgeschwindigkeit R des Verfahrens definieren wir in Abhängigkeit von M, q und ω:

R(M, q, ω) = − log\{sup_{e^{(0)}} lim sup_{k→∞} ||e^{(k)}||^{\frac{1}{k}} \}.  (3.27)

Formal stimmt diese Definition der Konvergenzgeschwindigkeit (2.22) des SOR-Verfahrens überein. In Satz 3.20 wird ein inhaltlicher Zusammenhang beider Größen hergestellt. e^{(k)} = z^{(k)} − z ist der Fehler im k-ten Iterationsschritt. Dabei sind z und y Lösungsvektoren von (3.22). B_{ω}(M) sei die Iterationsmatrix des SOR-Verfahrens bzgl. Mw = −q mit Relaxationsparameter ω. Für I := {1, ..., n} ist I(z) := {i ∈ I : z_{i} ≠ 0} der Träger von z. e^{(k)} hängt ab von M, q, ω und e^{(0)}. Für beliebiges e^{(0)} gilt: e^{(k)}_{\overrightarrow{k⟶∞}}0.

Lemma 3.19 : z^{(0)} sei gegeben. Dann existiert ein k_{0} ∈ IN, so daß für alle k ≥ k_{0} gilt: z_{i}^{(k)} > 0, wenn i ∈ Z := I(z) z_{i}^{(k)} = 0 und y_{i}^{(k)} > 0, wenn i ∈ Y := I(y). (3.28)

Beweis: Sei z^{(0)} gegeben. r^{(k)} wird im PSOR-Algorithmus berechnet. Aus z^{(k)}_{\overrightarrow{k⟶∞}}z und r^{(k)}_{\overrightarrow{k⟶∞}} − y ergibt sich die Existenz von ε > 0 und m ∈ IN,

so daß für alle k ≥ m gilt:

z_{i}^{(k)} > 0, wenn i ∈ Z
y_{i}^{(k)} > 0 und r_{i}^{(k)} < −ε, wenn i ∈ Y.

(3.28) erhält man nun mit Hilfe von (3.25):

Ist z_{i}^{(m)} ≤ \frac{ωε}{m_{ii}} ⇒ z_{i}^{(k)} = 0 ∀ k ≥ m.

Sei darum z_{i}^{(m)} > \frac{ωε}{m_{ii}} ⇒ z_{i}^{(m+1)} = z_{i}^{(m)} + ω \frac{r_{i}^{(k+1)}}{m_{ii}} < z_{i}^{(m)} − \frac{ωε}{m_{ii}} .

l_{i} ∈ IN sei Index mit z_{i}^{(m)} < l_{i} \frac{ωε}{m_{ii}} und z_{i}^{(m+l)} > 0 ∀ l < l_{i}.

Für alle i ∈ Y gilt z_{i}^{(m)} < max_{i∈Y} \{\frac{l_{i}ωε}{m_{ii}}\} .

Setze \overline{l} := max_{i∈Y} l_{i}, dann folgt \overline{l} > max_{i∈Y} \{\frac{m_{ii}z_{i}^{(m)}}{εω} \}.

Für k ≥ m + \overline{l} + 1 ist (3.28) somit erfüllt.

Mit k_{0} := m + 1 + \overline{l} folgt die Behauptung. □

Satz 3.20 : M und q seien so gewählt, daß |z_{i}| + |y_{i}| > 0, (3.29) 1 ≤ i ≤ n. z, y ∈ IR^{n} sind Lösungsvektoren von (3.22). Weiter sei Z := I(z), dann gilt: R(M, q, ω) = \{ \begin{matrix}∞ & , falls Z = ∅ \\ R_{∞}(B_{ω}[M_{ZZ}]) & , sonst \end{matrix} R_{∞}(B_{ω}[M_{ZZ}]) ist die Konvergenzgeschwindigkeit des SOR-Verfahrens mit der SOR- Matrix B_{ω}(M_{ZZ}).

Beweis: z^{(0)} sei gegeben. k_{0} sei aus Lemma 3.19. Für alle k ≥ k_{0} ergibt sich mit Lemma 3.19, (3.23), (3.24) und (3.29):

e_{Z}^{(k+1)} = B_{ω}[M_{ZZ}]) e_{Z}^{(k)}

e_{Y}^{(k+1)} = 0.

Aus der ersten Gleichung folgt nach Abschnitt 2.4.3 und (3.27), daß R(M, q, ω) = R_{∞}(M_{ZZ}). Aus der zweiten Gleichung folgt, daß man unmittelbar ein Lösung hat, somit R(M, q, ω) = ∞. □ Eigenschaft (3.29) erfüllen fast alle Matrizen M und Vektoren q. Nach Satz 3.20 können wir also einen optimalen Relaxationsparameter bestimmen wie beim SOR-Verfahren (vergleiche Kapitel 2).

Marktgleichgewicht (Economic Equilibrium)

4.1 Was ist Marktgleichgewicht ?

Netzwerke, die wir mit LKP-Verfahren lösen wollen, gehören zu den Marktgleichgewichtsmodellen. In Abschnitt 4.1.1 soll verdeutlicht werden, was ein ökonomisches Modell ist.

4.1.1 Ein kurzer Exkurs in die Makroökonomik

Grundsätzliches

Ökonomisches Handeln bedeutet Anbieten und Nachfragen. Beides findet auf Märkten statt. Durch die interdependente Untersuchung von Angebot und Nachfrage auf verschiedenen Märkten ist es möglich, gesamtwirtschaftliche Größen zu analysieren. Beispiele für gesamtwirtschaftliche Größen sind Arbeitslosigkeit, Inflation oder Preisniveau.

Ein Markt ist charakterisiert durch die wirtschaftlich handelnden ( Wirtschaftssubjekte ) und durch bestimmte Umweltfaktoren.

Da es unmöglich ist, ein System zu analysieren, in dem jeder einzelne Markt Berücksichtigung findet, werden Elementarmärkte zu Aggregaten zusammengefaßt. Üblich ist folgende Einteilung:

1. Gütermarkt: die Zusammenfassung aller produzierten Güter und Dienstleistungen

2. Arbeitsmarkt

3. Geldmarkt

4. Wertpapiermarkt

Wirtschaftssubjekte, also die wirtschaftlich handelnden, können sein:

KAPITEL 4. MARKTGLEICHGEWICHT (ECONOMIC EQUILIBRIUM)

1. Private Haushalte 2. Unternehmungen

3. Staaten

Ökonomische Modelle unterscheiden sich in den berücksichtigten Märkten, den wirtschaftlich handelnden und den unterstellten Umweltfaktoren. Sind alle Märkte im Modell berücksichtigt, so spricht man von einem Totalmodell. Untersucht man nur einen einzelnen Markt, dann macht man eine Partialanalyse, wobei man Abhängigkeiten zwischen den Märkten ausschließt.

Bestandteile eines ökonomischen Modells

In einem ökonomischen Modell sind Angebot und Nachfrage definiert durch Verhaltensgleichungen. Eine solche Verhaltensgleichung ist z.B. die Hypothese, daß die Nachfrage x auf einem Gütermarkt durch die Höhe des Güterpreises π bestimmt ist: x = x(π). Verhaltensgleichungen können frei gewählt werden, solange sie plausibel oder empirisch belegt sind; sie unterliegen jedoch der Budgetrestriktion. Auf einem Gütermarkt bedeutet dies, daß die Summe der Preise aller gekauften Güter der aller verkauften Güter entsprechen muß.

Zu jedem ökonomischen Modell gehören neben den Verhaltensgleichungen noch Identitäten und Relationen. Im einzelnen:

• \underline{Identitäten} sind Gleichungen, die wegen den verwendeten Definitionen stets erfüllt sind. So ist die Gesamtnachfrage x nach einem Gut identisch gleich der Summe der Nachfragen aller Wirtschaftssubjekte: x ≡ ∑_{i=1}^{n} x_{i}.

• \underline{institutionelle Relationen} spiegeln gesetzliche Zusammenhänge wider. So hängt etwa die Steuerzahlung eines Wirtschaftssubjektes von seinem Einkommen ab.

• \underline{technische Relationen} spiegeln technische Zusammenhänge wider; klassisches Beispiel hierfür sind die Produktionsfunktionen, nach denen die Menge der produzierten Güter von der Menge der verarbeiteten Rohstoffe abhängig ist.

Um das Modell vollständig zu machen, muß schließlich noch das Verhältnis festgelegt werden, in dem Angebot und Nachfrage stehen sollen.

Marktgleichgewicht

Stimmen Angebot und Nachfrage überein, so ist der Markt im Gleichgewicht. Dies liegt z.B. dann vor, wenn man unterstellt, daß die Wirtschaftssubjekte mit den angebotenen/nachgefragten Mengen auf die Preise reagieren. Genauer gesagt heißt das: herrscht Nachfrageüberschuß, so konkurrieren die Nachfrager den Preis hinauf, herrscht Angebotsüberschuß, konkurrieren die Anbieter den Preis hinunter. Preisvariationen führen schließlich zu einem Gleichgewicht. Man nennt diese Reaktion den walrasianischen Anpassungsmechanismus.

Netzwerkoptimierung ist ein typisches Marktgleichgewichtsmodell. Bei der Optimierung einer Netzwerkstruktur betrachtet man einen Gütermarkt, auf dem gleichrangige Wirtschaftssubjekte (i.f. bezeichnet als Regionen) mit einem einzigen Gut handeln. Produktion findet nicht statt, jede Region verfügt über einen Vorrat an Handelsware. Bestimmt werden sollen die benötigten Transportmengen zwischen den Regionen und insbesondere die regionalen Güterpreise, so daß der Markt im Gleichgewicht ist. ^{1}

4.1.2 Netzwerkoptimierung

Wir führen eine Partialanalyse durch und betrachten einen reinen Gütermarkt. Die Regionen 1 bis n handeln mit einer einzigen Ware. Die Nachfrage d_{i} und das Angebot s_{i} in Region i hängen vom dortigen Preis π_{i} des Gutes ab: d_{i} = d_{i}(π_{i}) und s_{i} = s_{i}(π_{i}). Der Nachfrageüberschuß, oder auch die effektive Nachfrage x_{i} ergibt sich dann als:

x_{i}(π_{i}) = d_{i}(π_{i}) − s_{i}(π_{i}).  (4.1)

x : IR^{n} ⟶ IR^{n}, x(π) := (x_{1}(π_{1}), ..., x_{n}(π_{n})) ist eine Funktion, die abhängt von den Preisen in allen Regionen. Die x_{i}(π_{i}) sind Verhaltensgleichungen, für die gilt:

i) Die x_{i} sind konvex, stetig differenzierbar und monoton fallend (x(π) ≤ x(π̂) ∀ π̂ ≤ π).

ii) Die effektive Nachfrage bei einem Preis von 0 ist nicht-negativ: x_{i}(0) ≥ 0.

iii) Es gibt einen Preisvektor π̂ ∈ IR^{n}, π̂ ≥ 0 mit ∑_{i=1}^{n} x_{i}(π̂) < 0.

Läßt man Im/Exporte nicht zu, dann ist x_{i}(π_{i}) = 0 (Nachfrage hat keinen Sinn), und die Gleichgewichtspreise ergeben sich als Schnittpunkte der Angebots/Nachfragefunktionen: d_{i}(π_{i}) = s_{i}(π_{i}) ( siehe dazu Abb. 5 ).

Um eine Einheit des Gutes von Region i nach Region j zu befördern, seien Transportkosten in Höhe von t_{ij}, i ≠ j aufzuwenden. Der direkte Transport soll kostengünstiger sein, als der Transfer über eine dritte Region. Das heißt, es gilt die △-Ungleichung

t_{ij} + t_{jk} ≥ t_{ik} ∀ i, j, k  (4.2)

und es sind die Transportbedingungen erfüllt:

π_{i} + t_{ij} ≥ π_{j}  (4.3)
(π_{i} − π_{j} + t_{ij})z_{ij} = 0 ∀ i, j.  (4.4)

z_{ij} bezeichnet dabei die Stückzahl des Gutes, die von Region i nach Region j transportiert wird. Nach (4.3) ist der Preis π_{j}, der in j erzielt wird, nicht höher als die Ausgaben, das Gut in i zu kaufen und nach j zu transportieren. Ein Export von i nach j lohnt also nicht mehr. Ist das Netzwerk im Gleichgewicht und findet (noch) Export statt (z_{ij} > 0), dann gilt (4.3) als Gleichung. Gilt in (4.3) keine Gleichheit, dann wird auch nicht exportiert (z_{ij} = 0). (4.3) soll in Matrizenschreibweise formuliert werden. Definiere dazu die Beschränkungsmatrix A:

1Wie diese Werte volkswirtschaftlich interpretiert und ausgewertet werden, soll hier nicht erläutert werden (zu diesem Zweck vergleiche man [Fuhrmann ]).

A = ( \begin{matrix}1 & −1 \\ ⋮ & & ... \\ ⋮ & & & ... \\ 1 & & & & −1 \end{matrix} \} n − 1 ; \begin{matrix}−1 & 1 \\ & ⋮ & −1 \\ & ⋮ & & ... \\ & 1 & & & −1 \end{matrix} \} n − 1 ; ⋮  ⋮ ⋮  ⋮ ; \begin{matrix}−1 & & & & 1 \\ & ... & & & ⋮ \\ & & ... & & ⋮ \\ & & & −1 & 1 \end{matrix} \} n − 1 ) ∈ IR^{n(n−1)×n}.

Ordnet man die Vektoren des IR^{n(n−1)} so an: (t_{12}, ..., t_{1n}, t_{21}, t_{23}, ..., t_{2n}, ..., t_{n1}, ..., t_{n,n−1})^{t},dann ist (4.3) äquivalent zu

Aπ ≥ −t.  (4.5)

Wir halten fest, was in diesem Kapitel als “Netzwerkproblem” bezeichnet wird: Bestimme π^{∗} und z^{∗}, so daß

x(π^{∗}) + A^{t}z^{∗} = 0  (4.6)
Aπ^{∗} ≥ −t, π^{∗} ≥ 0  (4.7)
π^{∗}^{t}A^{t}z^{∗} + t^{t}z^{∗} = 0, z^{∗} ≥ 0.  (4.8)

−A^{t}z ist der Nettoimport jeder Region, d.h.

−(A^{t}z)_{k} = − ∑_{{i=1 \atop i≠k}}^{n} z_{ki} + ∑_{{i=1 \atop i≠k}}^{n} z_{ik}.  (4.9)

−(A^{t}z)_{k} ist die innerhalb des gesamten Netzwerkes nach k transportierte Warenmenge abzüglich der Menge, die von Region k in alle anderen Regionen exportiert wird. (4.6) ist die Marktgleichgewichtsbedingung, was bedeutet, daß die effektive Nachfrage x_{k}(π_{k}) von Region k ihrem Nettoimport entsprechen soll. (4.7) und (4.8) sind die Transportbedingungen.

4.1.3 Der Baum-Wald-Algorithmus

Der Baum-Wald-Algorithmus ist ein Verfahren zur Lösung des Netzwerkproblems. Er hat nichts mit linearer Komplementarität zu tun, und soll an dieser Stelle zur näheren Definition des Problems dienen. Insofern wird auf Beweise verzichtet und der Artikel von [Glassey ] zur Vertiefung empfohlen. Dem Verfahren liegt folgende Anschauung zugrunde: Jede Region k im Netzwerk stellen wir als einen Punkt dar. Findet ein Export von j nach k statt (z_{jk} > 0), dann sind die Punkte j und k über den Ast (j, k) miteinander verbunden.

Definition 4.1 : Ein Baum besteht aus k Punkten 1 bis k, die mit k − 1 Ästen miteinander verbunden sind: Jeder Ast hat eine Richtung. Den Ast, der i und j in Richtung von i nach j verbindet, bezeichnen wir mit (i, j). Ein Baum enthält keine geschlossenen Wege, d.h. es existiert kein N ∈ IN, so daß der Baum die Äste (i_{l−1}, i_{l}), l = 1, ..., N enthält mit {i_{0}, ..., i_{N}} ⊆ {1, ..., k} und i_{0} = i_{N}. Ein Wald ist eine Menge von Bäumen.

Lemma 4.2 : Es existiert eine Lösung z^{∗} des Netzwerkproblems mit folgender Eigenschaft: lauten die Gleichgewichtswerte z_{ij}^{∗}, i, j ∈ {1, ..., n} und trägt man alle Äste (i, j) mit z_{ij}^{∗} > 0 in die Graphik ein, dann ergibt sich ein Wald.

Beweis: □

Hat man also z.B. 7 Regionen, dann hat z^{∗} möglicherweise die Darstellung

[figure] [figure] [figure] [figure] [figure] [figure]

1

2 3

4

5 6

7

Abb. 4 : 7 Regionen-Beispiel

Die zu einem Baum verbundenen Regionen bilden jeweils eine Handelsgruppe. An dieser Darstellung erkennt man, daß eine Region entweder exportiert oder importiert – jedoch nicht beides. Die Programmierung des Baum-Wald-Algorithmus erfolgt nach dem Prinzip, daß solange Handelsgruppen aufgebaut werden, bis eine Gleichgewichtslösung bestimmt ist.

Der Baum-Wald-Algorithmus

Schritt 0:

Initialisierung

Setze π_{i} = a_{i}, y_{i} = 0, i = 1, ..., n; z_{ij} = 0 für alle Äste (i, j). D.h. jeder Knoten bildet eine Handelsgruppe von einem Teilnehmer.

Schritt 1: Suche zwei Knoten j und k, so daß π_{j} + t_{jk} < π_{k}.

Existiert kein Paar mit dieser Eigenschaft ⇒ die aktuellen Werte sind eine Gleichgewichtslösung.

Schritt 2: Die Handelsgruppe E enthalte j, die Handelsgruppe I enthalte k. Suche die Exportregion e ∈ E und die Importregion i ∈ I, so daß

\overline{t}_{ei} = π_{e} + t_{ei} − π_{i} = min_{j∈E,k∈I}\{π_{j} + t_{jk} − π_{k}\}

(a) z_{ei} steigt bei Einhaltung der Gleichgewichtsbedingungen an, bis entweder (b) oder (c) erfüllt ist:

(b) π_{e}(z_{ei}) + t_{ei} = π_{i}(z_{ei}). π_{k}(z_{ei}) ist dabei der Gleichgewichtspreis in Region k als Funktion in z_{ei}. E und I werden mit dem Ast (e, i) verbunden und bilden eine neue Handelsgruppe. Gehe zu Schritt 1.

(c) z_{jk} = 0 für einen Ast (j, k) in E oder I. Trenne diese Handelsgruppe. Wenn dieser Ast sich in I befand, nenne den Teilbaum, der i enthält, I, befand sich dieser Ast in E, dann bezeichne den Teilbaum, der e enthält mit E. Gehe zu Schritt 2(a).

Bemerkungen zum Algorithmus und Hinweise für die Programmierung finden sich in Kapitel 6.

Satz 4.3 : In einer endlichen Zahl von Iterationsschritten führt der Baum-Wald-Algorithmus zur Lösung des Netzwerkproblems.

Beweis: Der Beweis wird geführt, indem man das Netzwerkproblem als Optimierungsaufgabe formuliert mit einer Zielfunktion Φ. Der Wert von Φ steigt bei Durchlauf der Schritte 1 und 2 strikt an. Verschiedene Werte von Φ ergeben eine neue Zusammensetzung der Handelsgruppen. Da es nur endlich viele Zusammensetzungen gibt, bricht das Verfahren ab. □

4.2 Verfahren zur Netzwerkoptimierung

4.2.1 “Quadratisches Netzwerk”

Wir wollen das Netzwerkproblem (4.6)-(4.8) in der Weise umformen, daß sich ein LKP (q, M) ergibt. Die Verhaltensgleichungen d_{i}(π_{i}) und s_{i}(π_{i}) aus Abschnitt 4.1.2 seien linear.

Es folgt, daß auch die effektive Nachfrage x(π) eine Funktion der Form x(π) = Dπ + c ist. D ∈ IR^{m×m} ist dabei eine Diagonalmatrix. Ist y_{k} = −(A^{t}z)_{k}, k = 1, ..., n die Nettoimportmenge aus (4.9), dann schreibt sich (4.6):

π_{i} = a_{i} − b_{i}y_{i},  (4.10)

mit passenden Konstanten a_{i} und b_{i}. Im Modell ist a_{i} der (positive) Gleichgewichtspreis ohne Im/Exporte und b_{i} die (positive) Elastizität von Angebot und Nachfrage. a_{i} und b_{i} lassen sich aus d_{i}(π_{i}) und s_{i}(π_{i}) herleiten; a_{i} ist der Preis im Schnittpunkt von d_{i}(π_{i}) und s_{i}(π_{i}):

[figure] [figure]

Abb. 5 : a_{i} als Schnittpunkt von d_{i}(π_{i}) und s_{i}(π_{i}).

d_{i}/s_{i}

π_{i}

d_{i}(π_{i})

s_{i}(π_{i})

a_{i}

Ist y_{i} < 0, so wird y_{i} nicht als Import, sondern als Export interpretiert. t und z seien definiert wie in 4.1.2. Es gelten (4.6)-(4.8). Mit den eingeführten Bezeichnungen stellt sich das Netzwerkproblem nun folgendermaßen: Bestimme aus den gegebenen Daten (a_{i}, b_{i}, t_{ij}, i = 1, ..., n; j = 1, ..., n) die Gleichgewichtspreise und -mengen. Also: Gesucht ist ein z^{∗} ∈ IR^{n(n−1)} und ein π^{∗} ∈ IR^{n} mit:

\begin{aligned}z_{ij}^{∗} & ≥ & 0 \\ π_{i}^{∗} + t_{ij} − π_{j}^{∗} & ≥ & 0 \\ z_{ij}^{∗}(π_{i}^{∗} + t_{ij} − π_{j}^{∗}) & = & 0 ∀ i, j. \end{aligned}  (4.11)

(4.9) und (4.10) führen auf die Zuweisung:

π_{i} = a_{i} − b_{i}(∑_{j=1}^{n} z_{ji} − ∑_{j=1}^{n} z_{ij}) i = 1, ..., n.  (4.12)

In Abschnitt 4.1.2 wurde die Beschränkungsmatrix A ∈ IR^{n(n−1)×n} eingeführt. Setzen wir die vorgegebenen Werte ein in die Diagonalmatrix B = ( \begin{matrix}b_{1} & & & & 0 \\ & . & . & . \\ 0 & & & & b_{n} \end{matrix} ) ∈ IR^{n×n}

und den Vektor (a_{1}, ..., a_{n}) ∈ IR^{n}, dann schreibt sich (4.12) in der Form

π = BA^{t}z + a.

In Matrizenform lautet (4.11):

z ≥ 0
Aπ ≥ −t
z^{t}(Aπ + t) = 0.

Setze π = BA^{t}z + a ein ⇒

\begin{aligned}z & ≥ & 0 \\ ABA^{t}z + Aa + t & ≥ & 0 \\ z^{t}(ABA^{t}z + Aa + t) & = & 0. \end{aligned}  (4.13)

Also erhält man eine Lösung von (4.11), indem man das LKP (q, M) löst mit q = Aa + t und M = ABA^{t}. Die Vektoren des IR^{n(n−1)} sind nach 4.1.2 so angeordnet: (ξ_{12}, ..., ξ_{1n}, ξ_{21}, ξ_{23}, ..., ξ_{n1}, ..., ξ_{n,n−1})^{t}.

M ∈ IR^{n(n−1)×n(n−1)} ist positiv-semidefinit, da die b_{i} > 0 sind.

4.2.2 Nicht-lineares Netzwerk

Im letzten Abschnitt ging es um lineares Netzwerkoptimieren mit der effektiven Nachfrage x(π) = Dπ + c. Mit passenden Zuweisungen ergab sich aus (4.6)-(4.8) ein LKP (q, M) . Jetzt betrachte man (4.6)-(4.8) und x(π) sei eine nicht-lineare Funktion in dem Sinn, daß jede Komponente eine nicht-lineare Funktion ist: x(π) = (x_{1}(π_{1}), ..., x_{n}(π_{n}))^{t}. Nach wie vor gelten die Eigenschaften i) - iii) aus Abschnitt 4.1.2. Es seien:

y := (π, z)^{t} und G(y) := ( \begin{matrix}−x(π) − A^{t}z \\ Aπ + t \end{matrix} ).

Ist y^{∗} = (π^{∗}, z^{∗}) eine Lösung des nichtlinearen Komplementaritätsproblems

\begin{aligned}G(y) & ≥ & 0 \\ y & ≥ & 0 \\ y^{t}G(y) & = & 0, \end{aligned} (NKP)

so erfüllt y^{∗} die Gleichgewichtsbedingungen (4.6)-(4.8) unter der Voraussetzung π^{∗} > 0. Diese Bedingung ist bereits bei recht schwachen Voraussetzungen erfüllt, wie Lemma 4.4 zeigt:

Lemma 4.4 : Für alle i sei x_{i}(0) > 0. Löst y^{∗} = (π^{∗}, z^{∗}) das (NKP) ⇒ π^{∗} > 0.

Beweis: Es gilt x_{i}(π_{i}) + A_{i}^{t}z ≤ 0 ( nach der Definition des (NKP) ). Hinreichend ist zu zeigen: π_{i} = 0 impliziert A_{i}^{t}z ≥ 0 (Was bedeutet, daß Region i bei einem Preis von 0 nicht importieren kann). Annahme: A_{i}^{t}z < 0. Dann existiert ein j mit z_{ji} > 0. Aus der Komplementärbedingung in (4.11) folgt: π_{i} − π_{j} = t_{ij}. π_{i} = 0 wurde angenommen, also π_{j} < 0. Dies ist ein Widerspruch. □

Definition 4.5 : X sei Teilmenge eines normierten linearen Raumes, X′ sei sein Dualraum. C ⊆ X bezeichne eine nicht-leere, abgeschlossene konvexe Teilmenge von X; bzgl. C definiert man N_{C}(·) als den Normalen-Kegel-Operator : N_{C}(x) = \{ \begin{matrix}\{y ∈ X′ | y^{t}(c − x) ≤ 0 ∀ c ∈ C\} , x ∈ C \\ ∅ , x ∉ C. \end{matrix} Eine verallgemeinerte Gleichung ist eine Inklusion der Form 0 ∈ g(x) + N_{C}(x), mit g : X ⟶ X′.

Mit den Bezeichnungen aus Definition 4.5 setze: X = IR^{n}, X′ = IR^{n}, C = IR_{+}^{n}. Bestimmt werden soll eine Lösung von

0 ∈ G(y) + NIR_{+}^{n}(y).  (4.14)

Lemma 4.6 : y^{∗} ist Lösung des (NKP) ⇔ y^{∗} löst (4.14).

Beweis: ′′ ⟹′′

y^{∗} sei eine Lösung des (NKP). ⇒ −G(y^{∗}) ∈ NIR_{+}^{n}(y^{∗}), denn für c ∈ IR_{+}^{n} beliebig ist: −G(y^{∗})^{t}(c − y^{∗}) = −G(y^{∗})^{t}c ≤ 0.

Setze u := −G(y^{∗}), dann gilt 0 = G(y^{∗}) + u und 0 ∈ G(y^{∗}) + NIR_{+}^{n}(y^{∗}) ist bewiesen.

′′ ⟸′′

y^{∗} löse (4.14) ⇒ ∃ u ∈ NIR_{+}^{n}(y^{∗}) mit 0 = G(y^{∗}) + u. Insbesondere ist NIR_{+}^{n}(y^{∗}) nicht leer, also y^{∗} ∈ IR_{+}^{n}.

u^{t}(c − y^{∗}) ≤ 0 ∀ c ∈ IR_{+}^{n} ⇒ (c − y^{∗})^{t}G(y^{∗}) ≥ 0 ∀ c ∈ IR_{+}^{n}.

Sei i ∈ {1, ..., n} und definiere c̃_{j} = y_{j}^{∗} ∀ j ≠ i und c̃_{i} = 1 + y_{i}^{∗} also c ∈ IR_{+}^{n}.

Setzt man dieses c̃ oben für c ein, so folgt G_{i}(y^{∗}) ≥ 0 ⇒ G(y^{∗}) ≥ 0.

Mit c = 0 erhält man −(y^{∗})^{t}G(y^{∗}) ≥ 0. Aus y^{∗} ≥ 0 und G(y^{∗}) ≥ 0 folgt (y^{∗})^{t}G(y^{∗}) = 0. □

Um die Bedeutung von Lemma 4.6 für die Lösung des (NKP) zu erklären, sei zunächst kurz an das Newton-Verfahren erinnert: Man will eine Nullstelle der nicht-linearen Funktion g : IR^{n} ⟶ IR^{n} berechnen. Ist z^{(k)} eine Näherung an diese Nullstelle, so erhält man eine bessere Näherung z^{(k+1)} als Lösung des Systems 0 = g(z^{(k)}) + ∇g(z^{(k)})(z − z^{(k)}).

Dieses Linearisierungsprinzip soll auf (4.14) angewendet werden. Ist also y^{(k)} eine Approximation an die Lösung von (4.14), so ergibt sich eine bessere Approximation y^{(k+1)} als Lösung der linearisierten verallgemeinerten Gleichung

0 ∈ L_{k}(y) + NIR_{+}^{n}(y)  (4.15)

mit L_{k}(y) = G(y^{(k)}) + ∇G(y^{(k)})(y − y^{(k)}). Konvergiert {y^{(k)}}_{k} gegen y^{∗}, so ist y^{∗} Lösung von (4.14). y^{(k+1)} wiederum ist nach Lemma 4.6 eine Lösung von (4.15) genau dann, wenn y^{(k+1)}

y ≥ 0, L_{k}(y) ≥ 0, y^{t}L_{k}(y) = 0  (4.16)

k = 0, 1, 2, 3, ... löst.

{y^{(k)}}_{k} ist also eine Folge von LKP-Lösungen.

Zu zeigen bleibt ihre Konvergenz.

Mit den Begriffen der Netzwerkstruktur kann der Konvergenzbeweis recht anschaulich durchgeführt werden. y = (π, z)^{t} ist i.f. gegeben, y^{+} = (π^{+}, z^{+})^{t} wird berechnet. y^{+} ist Lösung des LKP (q, M) mit

q = ( \begin{matrix}−x(π) + π^{t}∇x(π) \\ t \end{matrix} ) und M = ( \begin{matrix}−∇x(π) & −A^{t} \\ A & 0 \end{matrix} ).

Setzt man y = y^{(k)} und y^{+} = y^{(k+1)}, dann ist y^{+} eine Lösung von (4.16).

Da x(π) nicht-linear ist, entsteht bei der Iteration eine nicht-konstante Folge von Preisvektoren. Aus den folgenden Lemmata geht hervor, daß diese Folge quadratisch konvergiert. Im weiteren bezeichnen wir A_{i}^{t}z := (A^{t}z)_{i} ∀ z ∈ IR^{n(n−1)}.

Lemma 4.7 : Nach jedem Iterationsschritt überschreitet oder entspricht der Nachfrageüber- schuß in Region i die dorthin erfolgenden Nettoimporte, d.h. x_{i}(π_{i}^{+}) ≥ −A_{i}^{t}z^{+} ∀ i.

Beweis: Ist π_{i}^{+} = 0, dann argumentiert man wie in Lemma 4.4, um zu zeigen, daß Region i nicht importieren kann, somit A_{i}^{t}z^{+} ≥ 0 erfüllt ist. Vorausgesetzt wurde ( Bedingung ii) ) x_{i}(0) ≥ 0. ⇒ x_{i}(0) + A_{i}^{t}z^{+} ≥ 0. Ist π_{i}^{+} > 0, dann ergibt sich mit der Komplementaritätsbedingung:

(q + My^{+})^{t}y^{+} = 0. ⇒ π_{i}^{+}(−x_{i}(π_{i}) + x′_{i}(π_{i})(π_{i} − π_{i}^{+}) − A_{i}^{t}z^{+}) = 0.

⟹^{π_{i}^{+}>0} −x_{i}(π_{i}) − A_{i}^{t}z^{+} + x′_{i}(π_{i})(π_{i} − π_{i}^{+}) = 0.

Da die x_{i} konvex sind, folgt mit dem Satz über konvexe Funktionen x_{i}(π_{i}^{+}) ≥ x_{i}(π_{i})+x′_{i}(π_{i})(π_{i}^{+}−π_{i}). D.h. x_{i}(π_{i}^{+}) + A_{i}^{t}z^{+} ≥ 0. □

Lemma 4.8 : Aus π_{i}^{+} < π_{i} folgt A_{i}^{t}z > A_{i}^{t}z^{+}. D.h. wenn in einer Region der Preis sinkt, dann steigt der dorthin erfolgende Nettoimport.

Beweis: Annahme: A_{i}^{t}z ≤ A_{i}^{t}z^{+}.

⟹^{Lemma 4.7} x_{i}(π_{i}) + A_{i}^{t}z^{+} ≥ x_{i}(π_{i}) + A_{i}^{t}z ≥ 0

(π und z wurden im vorangegangenem Iterationsschritt berechnet, Lemma 4.7 gilt also auch für diese Werte). Aus der Bestimmungsgleichung für y^{+} erhalten wir:

x_{i}(π_{i}) + x′_{i}(π_{i})(π_{i}^{+} − π_{i}) + A_{i}^{t}z^{+} ≤ 0 ⇒ x′_{i}(π_{i})(π_{i}^{+} − π_{i}) ≤ 0.

Aus x′_{i}(π_{i}) < 0 ( denn x_{i}(π_{i}) ist für alle i monoton fallend ) folgt nun π_{i}^{+} ≥ π_{i}. □

Lemma 4.9 : Gilt z_{kj}^{+} > 0 und π_{j} > π_{j}^{+}, so folgt π_{k} > π_{k}^{+}. Gilt z_{jk} > 0 und π_{j} > π_{j}^{+}, so folgt π_{k} > π_{k}^{+}.

Beweis: Im ersten Fall gilt: π_{k} + t_{kj} ≥ π_{j} und π_{k}^{+} + t_{kj} = π_{j}^{+}. ⇒ π_{k} − π_{k}^{+} ≥ π_{j} − π_{j}^{+} > 0. Im zweiten Fall haben wir: π_{j} + t_{jk} = π_{k} und π_{j}^{+} + t_{jk} ≥ π_{k}^{+}. ⇒ 0 > π_{j}^{+} − π_{j} ≥ π_{k}^{+} − π_{k}. □

(4.15) hat für alle k eine Lösung y = (π, z)^{t} mit der Eigenschaft

{∃ j : z_{ij} > 0 ⇒ z_{ki} = 0 ∀ k}, bzw. {∃ k : z_{ki} > 0 ⇒ z_{ij} = 0 ∀ j}.

Im Klartext heißt das, daß eine Region entweder importiert oder exportiert, aber nicht beides.

Satz 4.10 : Für k = 0, 1, 2, 3, ... bilden die die ersten n Komponenten der Lösungen von (4.16) eine konvergente Folge von Preisvektoren.

Beweis: Ist die Funktion der effektiven Nachfrage x konstant oder linear, dann bricht das Verfahren (4.16) nach dem ersten Schritt ab. Für den Fall, daß x nicht-linear ist, zeigen wir, daß die Folge von Preisvektoren nach oben beschränkt und monoton steigend ist:

1.Beschränktheit:

Annahme: in der j-ten Komponente sei die Folge der Preisvektoren unbeschränkt. Aus (4.7) folgt dann, daß die Folge in jeder Komponente unbeschränkt ist.

⇒ In einem Iterationsschritt gilt ∑_{i=1}^{n} x_{i}(π_{i}) < 0 und ∑_{i=1}^{n} A_{i}^{t}z = 0.

⟹^{Lemma 4.7} ∑_{i=1}^{n} x_{i}(π_{i}) ≥ 0 in jedem Schritt der Iteration.

Das ist ein Widerspruch. 2.Monotonie: Annahme: ∃ j : π_{j}^{+} < π_{j} ⇒ ∃ Region i mit A_{i}^{t}z^{+} < 0 und π_{i}^{+} < π_{i}, denn: Gilt A_{j}^{+}z^{+} < 0, dann erfüllt j diese Eigenschaft, sonst: A_{j}^{+}z^{+} ≥ 0. Nach Lemma 4.8 gilt: A_{j}^{t}z > A_{j}^{t}z^{+} ⇒ A_{j}^{t}z > 0. Also existiert eine Region i, i ≠ j mit z_{ji} > 0 ⟹^{Lemma 4.9} π_{i} > π_{i}^{+}. i ist Importeur, fazit kein Exporteur, d.h. A_{i}^{t}z < 0 ⟹^{Lemma 4.8} A_{i}^{t}z^{+} < 0. Zum weiteren Beweis definiere folgende Indexmengen: I := {i | A_{i}^{t}z^{+} < 0 ∧ π_{i}^{+} < π_{i}} ist die Menge der Importeure mit fallenden Preisen. J := {j | z_{ji} > 0 für ein i ∈ I} sind die nach I exportierenden Regionen. K := {k | z_{jk} > 0 für ein j ∈ J} sind die Regionen, in die die Regionen aus J importieren.

⟹^{Lemma 4.9} π_{j}^{+} < π_{j} ∀ j ∈ J. ⇒ π_{k}^{+} < π_{k} ∀ k ∈ K. ⟹^{Lemma 4.8} A_{k}^{t}z > A_{k}^{t}z^{+} ∀ k ∈ K.

Alle k ∈ K sind Importeure, also keine Exporteure. ⇒ A_{k}^{t}z^{+} < 0 ∀ k ∈ K ⇒ K ⊆ I. Da die Regionen, in die die Regionen aus J exportieren, in K enthalten sind, ergibt sich die Ungleichung

∑_{k∈K} A_{k}^{t}z + ∑_{j∈J} A_{j}^{t}z ≤ 0. Mit Lemma 4.8 folgt ∑_{k∈K} A_{k}^{t}z^{+} + ∑_{j∈J} A_{j}^{t}z^{+} < 0.

Da K in I enthalten ist, ergibt sich ∑_{i∈I} A_{i}^{t}z^{+} + ∑_{j∈J} A_{j}^{t}z^{+} < 0.

Andererseits gilt aber ∑_{i∈I} A_{i}^{t}z^{+} + ∑_{j∈J} A_{j}^{t}z^{+} ≥ 0, und das ist ein Widerspruch.

Daher π_{j}^{+} ≥ π_{j} ∀ j. □

Satz 4.11 : Die Folge der Preisvektoren {π^{(k)}}_{k} konvergiere gegen π^{∗}. Ist ∇x(π^{∗}) nicht singulär und ist ∇x(π) lipschitzstetig in einer offenen konvexen Umgebung von π^{∗}, dann konvergiert das Verfahren quadratisch.

Beweis: U(π^{∗}) sei die Umgebung von π^{∗}, in der ∇x(π) nicht-singulär ist. ⇒ ∃ a_{2} > 0 : a_{2}|π^{+} − π^{∗}| ≤ |∇x(π)(π^{+} − π^{∗})| ∀ π^{+}, π ∈ U(π^{∗}). x(π) ist nicht ansteigend, also gilt ∇x(π) ≤ 0. {π^{(k)}}_{k} ist monoton steigend (nach Beweis von Satz 4.10)

⇒ (π^{+} − π^{∗}) ≤ 0 ⇒ a_{2}|π^{+} − π^{∗}| ≤ ∑_{i=1}^{n} x′_{i}(π_{i})(π_{i}^{+} − π_{i}^{∗}).

KAPITEL 4. MARKTGLEICHGEWICHT (ECONOMIC EQUILIBRIUM)

Aus den Komplementärbedingungen folgt: x′_{i}(π_{i}) + x′_{i}(π_{i})(π_{i} − π_{i}^{+}) − A_{i}^{t}z^{+} ≥ 0 ⇒ x′_{i}(π_{i})π_{i}^{+} ≤ −x_{i}(π_{i}) + x′_{i}(π_{i})π_{i} − A_{i}^{t}z^{+}.

⇒ ∑_{i=1}^{n} x′_{i}(π_{i})π_{i}^{+} ≤ ∑_{i=1}^{n} x′_{i}(π_{i})π_{i} − ∑_{i=1}^{n} x_{i}(π_{i}).

⇒ a_{2}|π^{+} − π^{∗}| ≤ ∑_{i=1}^{n} x′_{i}(π_{i})π_{i} − ∑_{i=1}^{n} x_{i}(π_{i}) − ∑_{i=1}^{n} x′_{i}(π_{i})π_{i}^{∗}.

Mit ∑_{i} x(π^{∗}) = 0 erhält man:

\begin{aligned}a_{2}|π^{+} − π^{∗}| & ≤ & ∑_{i} x_{i}(π_{i}^{∗}) − ∑_{i} x′_{i}(π_{i})(π_{i}^{∗} − π_{i}) − ∑_{i} x_{i}(π_{i}) ≤ \\ & ≤ & |x(π^{∗}) − x(π) − ∇x(π)(π^{∗} − π)|. \end{aligned}  (4.17)

∇x(π) ist lipschitzstetig ⇒ ∃ a_{1} > 0 mit

|x(π^{∗}) − x(π) − ∇x(π)(π^{∗} − π)| ≤ a_{1}|π − π^{∗}|^{2},  (4.18)

wenn π hinreichend nahe bei π^{+} ist. (4.17) und (4.18) ergeben zusammen die Behauptung. □

Der Konvergenzbeweis von Luo und Tseng (1991)

Die Konvergenzbeweise aus Kapitel 3 werden nun weniger restriktiv geführt. Die Bedingung, daß die Matrix M im LKP (q, M) positiv definit sein muß, wird abgeschwächt zu M positiv semidefinit. Für den Beweis setze die Prämisse

Die Funktion f(z) = q^{t}z + \frac{1}{2}z^{t}Mz (d.h. die Zielfunktion des (DQP)) ist auf dem Zulässigkeitsbereich X = {z ∈ IR^{n} : z ≥ 0} nach unten beschränkt.

(P)

( Es sei bemerkt, daß sich (P) nach Satz 2.17 erübrigt, wenn M positiv definit ist. ) Wir greifen zurück auf die Bezeichnungen in Abschnitt 3.1.1. Führt man bezüglich M das reguläre Splitting (B, C) durch, dann ist

z^{(k+1)} = [z^{(k+1)} − (Bz^{(k+1)} + Cz^{(k)} + q)]^{+}  (5.1)

eine Basis-Splitting Methode. Luo und Tseng bewiesen 1991, daß die BS-Methode (5.1) konvergiert, wenn (B, C) ein reguläres Splitting ist. Sie analysierten dazu das Verhalten der Iterationswerte am Rand des Zulässigkeitsbereiches X.

5.1 Eigenschaften regulärer Splittings

M ∈ IR^{n×n} ist in diesem Abschnitt stets symmetrisch und positiv semidefinit. (B, C) sei reguläres Splitting von M. Definiere die Funktion F_{B} : X ⟶ X

F_{B}(z) = {y ∈ IR^{n} : y = [y − (By + Cz + q)]^{+}} ∀ z ∈ X.  (5.2)

L sei die Lösungsmenge des LKP (q, M) .

KAPITEL 5. DER KONVERGENZBEWEIS VON LUO UND TSENG (1991)

Lemma 5.1 : Es existiert ein d^{∗} ∈ IR^{n} : Mz^{∗} + q = d^{∗} ∀ z^{∗} ∈ L.

Beweis: Es seien z_{1}, z_{2} ∈ L. Setze w_{i} = Mz_{i} + q, i = 1, 2. ⇒ w_{1} − w_{2} = M(z_{1} − z_{2}). Da M positiv semidefinit und z_{1}, z_{2} ∈ L : 0 ≤ (z_{1} − z_{2})^{t}M(z_{1} − z_{2}) = −z_{1}^{t}w_{2} − z_{2}^{t}w_{1} ≤ 0 ⇒ (z_{1} − z_{2})^{t}M(z_{1} − z_{2}) = 0.

Man betrachte die Funktion f(z) := z^{t}Mz mit M wie oben.

Gilt für z_{min} ∈ IR^{n} : f(z_{min}) = 0, so folgt ∇f(z_{min}) = 0, da z_{min} auf IR^{n} globales Minimum ist.

Aus ∇f(z) = 2Mz = (M + M^{t})z erhält man: (M + M^{t})(z_{1} − z_{2}) = 0. M ist symmetrisch, folglich Mz_{1} = Mz_{2}. ⇒ Behauptung. □

Lemma 5.2 : Es gelten folgende Aussagen: (a) F_{B} : X ⟶ X ist wohldefiniert und bildet jeden Punkt aus X auf einen Punkt aus X ab. (b) Es sei z ∈ X. Dann gilt: f(y) − f(z) ≤ \frac{1}{2}(y − z)^{t}(C − B)(y − z) mit y = F_{B}(z).

Beweis: \underline{zu (a):} Sei z ∈ X, dann ist zu beweisen, daß es genau ein y ∈ IR^{n} gibt mit y = [y − (By + Cz + q)]^{+}:

y = [y − (By + Cz + q)]^{+}

⟺
y ≥ 0 ∧ y ≥ y − (By + Cz + q) ∧ {(By + Cz + q)_{i} > 0 ⇒ y_{i} = 0}
⟺

\begin{aligned}y & ≥ & 0 \\ By + Cz + q & ≥ & 0 \\ y^{t}(By + Cz + q) & = & 0 \end{aligned}

y ist also Lösung des LKP (Cz + q, B).

Aus M positiv-semidefinit und (B − C) positiv definit folgt, daß B positiv definit ist. Nach Satz 2.30 existiert dann genau eine Lösung y.

\underline{zu (b):} M = B + C ⇒ f(y) − f(z) = (y − z)^{t}(By + Cz + q) + \frac{1}{2}(y − z)^{t}(C − B)(y − z) ∀ y, z ∈ IR^{n} (vergleiche Beweis zu Lemma 3.3).

Darüberhinaus gilt nach (5.2): y = F_{B}(z) ⇔ y ∈ X ∧ (By)_{i} + (Cz)_{i} + q_{i} ≥ 0 ∧ {(By)_{i} + (Cz)_{i} + q_{i} > 0 ⇒ y_{i} = 0} ∀ i = 1, ..., n. Demzufolge ist: −z^{t}(By + Cz + q) ≤ 0 ∧ y^{t}(By + Cz + q) = 0. Daraus ergibt sich die Behauptung. □

Für Q ∈ IR^{n×n} symmetrisch setze Kern(Q) = {z ∈ IR^{n} : Qz = 0} und span(Q) = {z ∈ IR^{n} : ∃ y ∈ IR^{n} : Qy = z}. Es ist klar, daß IR^{n} = Kern(Q) ⊕ span(Q) .

Was eine konvergente Matrix ist, war bereits Gegenstand von Definition 2.38. Es wurde gezeigt, daß das iterative Verfahren x^{(k+1)} = B(M) x^{(k)} zur Lösung von Mx = b konvergiert, wenn die Iterationsmatrix B(M) konvergent ist. Grundvoraussetzung war die Regularität von M. Ist M singulär, so ist die Konvergenz des Verfahrens nicht mehr äquivalent zur Konvergenz von B(M). Keller zeigte für diesen Fall, daß das Verfahren Nx^{(k+1)} = (N −M)x^{(k)}+b mit einer nicht-singulären Matrix N ∈ IR^{n×n} genau dann gegen eine Lösung von Mx = b konvergiert, wenn

• Kern(M) und ( Kern(M) )^{C} invariante Unterräume von M sind mit IR^{n} = Kern(M) ⊕ ( Kern(M) )^{C}.

• Für T = I − N^{−1}M gilt: T z = z ∀ z ∈ Kern(M) und ∥T^{k}z∥_{\overrightarrow{k⟶∞}}0 ∀ z ∉ Kern(M) .

(siehe [Keller ]). Wir formulieren diesen Satz als Definition:

Definition 5.3 : Sei Q ∈ IR^{n×n} symmetrisch. T ∈ IR^{n×n} heißt konvergent in Q, wenn 1. T z = z  ∀ z ∈ Kern(Q) . 2. ∥T^{k}z∥_{\overrightarrow{k⟶∞}}0  ∀ z ∈ span(Q) .

Im Beweis zu Lemma 5.5 verwenden wir die Hauptvektoren einer Matrix T ∈ IR^{n×n}.

Definition 5.4 : Es sei k ∈ IN. Gilt für v ∈ IR^{n} : (T −λI)^{k−1}v ≠ 0 und (T −λI)^{k}v = 0, dann heißt v Hauptvektor der Stufe k von T zum Eigenwert λ. Der von allen Hauptvektoren zu einem Eigenwert λ von T aufgespannte Teilraum heißt invarianter Unterraum von T zum Eigenwert λ.

Jede komplexe Matrix T ∈ C^{n×n} hat n linear unabhängige Hauptvektoren.

Lemma 5.5 : (Satz von Keller) Ist Q ∈ IR^{n×n} symmetrisch und N ∈ IR^{n×n} nicht singuläre Matrix, für die N + N^{t} − Q positiv definit ist. Dann gilt: T = I − N^{−1}Q ist konvergent in Q ⇔ Q ist positiv semi-definit.

Beweis: λ sei Eigenwert von T und v der entsprechende Eigenvektor: T v = λv. N ist nicht-singulär, also gilt für jede beliebige Matrix Q ∈ IR^{n×n} : Qv = (1 − λ)Nv. Klar ist: Qv = 0 ⇔ λ = 1. ′′ ⟹′′ T sei konvergent in Q. Wir nehmen an, Q ist nicht positiv-semidefinit; dann sei z^{(0)} ∈ IR^{n} ein Vektor mit z^{(0)}^{t}Qz^{(0)} < 0.

Definiert man z^{(k+1)} = T z^{(k)}, so läßt sich zeigen, daß {z^{(k)}^{t}Qz^{(k)}}_{k} monoton fallend ist. Es folgt z^{(k)}^{t}Qz^{(k)} ≤ z^{(0)}^{t}Qz^{(0)} < 0 ∀ k = 0, 1, 2, 3, ... z^{(0)} = u^{(0)} + w^{(0)} ist mit u^{(0)} ∈ Kern(Q) und w^{(0)} ∈ span(Q) eine eindeutige Zerlegung. T ist in Q konvergent, also z^{(k)} = T^{k}z^{(0)} = u^{(0)} + w^{(k)}, mit w^{(k)} = T^{k}w^{(0)}_{\overrightarrow{k⟶∞}}0.

⇒ z^{(k)}^{t}Qz^{(k)} = (u^{(0)} + w^{(k)})^{t}Q(u^{(0)} + w^{(k)}) = (u^{(0)} + w^{(k)})^{t}Qw^{(k)}_{\overrightarrow{k⟶∞}}0.

⇒ ∃ K ∈ IN : z^{(k)}^{t}Qz^{(k)} > z^{(0)}^{t}Qz^{(0)} ∀ k ≥ K. Dies ist ein Widerspruch !

′′ ⟸′′ Q sei positiv semi-definit. v sei ein Eigenvektor von T mit v ∈ span(Q) . ⇒ v ∉ Kern(Q) ⇒ Qv ≠ 0 ⇒ λ ≠ 1 und v^{t}Qv > 0.

Aus Qv = (1 − λ)Nv folgt \frac{1}{1 − λ} = \frac{v^{t}Nv}{v^{t}Qv}.

Addiere auf beiden Seiten das Komplexkonjugierte (N ist reell, d.h. v^{t}Nv = v^{t}N^{t}v):

2 Re \frac{1}{1 − λ} = \frac{v^{t}(N + N^{t})v}{v^{t}Qv}.

Mit den Bezeichnungen λ = α + iβund P = N + N^{t} − Q kommt man auf:

\frac{2(1 − α)}{(1 − α)^{2} + β^{2}} = 1 + \frac{v^{t}P v}{v^{t}Qv} > 1 ⇒ |λ|^{2} = α^{2} + β^{2} < 1.

span(Q) ist der invariante Unterraum von T, den die Hauptvektoren von T mit Eigenwerten ungleich 1 aufspannen.

Für diese Eigenwerte gilt also |λ| < 1, folglich ∥T^{k}z∥_{\overrightarrow{k⟶∞}}0 ∀ z ∈ span(Q) .

Es bleibt zu zeigen, daß es zum Eigenwert λ = 1 keine Hauptvektoren der Stufe größer als 1 gibt.

Dazu zeigt man, daß für die Folge z^{(k+1)} = T z^{(k)}, z^{(0)} beliebig, lim_{k⟶∞} Qz^{(k)} = 0 gilt. Nimmt man an, daß es zum Eigenwert λ = 1 Hauptvektoren einer Stufe ungleich 1 gibt, dann heißt das: Es existiert ein Hauptvektor w ∈ IR^{n} und ein u ∈ IR^{n} mit Qw ≠ 0, T w = w + u und T u = u. Setze z^{(0)} = w ⇒ z^{(k)} = w + ku. u ∈ Kern(Q) , also Qz^{(k)} = Qw. Dies ist ein Widerspruch zu lim_{k⟶∞} Qz^{(k)} = 0. □

Im Beweis von Lemma 5.6 und in Abschnitt 5.2 wird der symmetrische Teil einer Matrix verwendet; S ∈ IR^{n×n} ist symmetrischer Teil von M ∈ IR^{n×n}, wenn S symmetrisch ist und wenn gilt: x^{t}Sx = x^{t}Mx ∀ x ∈ IR^{n}.

Lemma 5.6 : Sei Q ∈ IR^{m×m} symmetrisch und positiv semi-definit. (N, H) sei reguläres Split- ting von Q. Dann ist N positiv definit und es gilt: (a) Es existiert ein ρ ∈ (0, 1) und ein τ > 0 mit ∥ ∏_{j=1}^{k}(I − θ^{j}QN^{−1})z∥ ≤ τ(1 − \overline{θ}(1 − ρ))^{k}∥z∥ ∀ k ≥ 1 ∀ z ∈ span(Q) dabei ist \overline{θ} ∈ (0, 1], und \{θ^{j}\}_{j} ist eine Folge von Skalaren aus [\overline{θ}, 1]. (b) Es gibt ein Δ ≥ 1 mit ∥ ∏_{j=1}^{k}(I − θ^{j}N^{−1}Q)z∥ ≤ Δ∥z∥ ∀ k ≥ 1 ∀ z ∈ IR^{m} \overline{θ} und \{θ^{j}\}_{j} wie oben.

Beweis: Aus Q = N + H folgt N − H = 2N − Q, d.h. der symmetrische Teil von N − H ist N + N^{t} − Q; dieser ist positiv definit, da N − H positiv definit ist. D.h. N + N^{t} ist positiv definit, und damit ist N positiv definit. Insbesondere existiert N^{−1}. Mit dem Satz von Keller folgt I − N^{−1}Q ist in Q konvergent. ⇒

(I − N^{−1}Q)z ∈ span (Q)  ∀ z ∈ span (Q),  (5.3)

und es gibt ein ρ ∈ (0, 1), das nur von N und Q abhängt mit

∥(I − N^{−1}Q)z∥ ≤ ρ∥z∥  ∀ z ∈ span (Q).  (5.4)

\underline{zu (a)}

Betrachte ein z ∈ span(Q) , ein k ≥ 1, ein \overline{θ} ∈ (0, 1] und eine Folge \{θ^{k}\}_{k} aus [\overline{θ}, 1]. Q ist symmetrisch positiv semidefinit, darum existiert ein w mit z = Qw und

∥w∥ ≤ \frac{1}{λ_{min}} ∥z∥.

λ_{min} ≥ 0 ist der kleinste Eigenwert von Q. w^{0} sei die orthogonale Projektion von w auf span(Q) . Die Vektoren w^{1}, ..., w^{k} werden rekursiv definiert:

w^{j} = (I − θ^{j}N^{−1}Q)w^{j−1}, j = 1, ..., k.  (5.5)

⇒ z = Qw^{0}, w^{0} ∈ span(Q) .

Mittels Induktion folgt aus (5.3) und (5.5) w^{j} ∈ span(Q) ∀ j = 0, 1, ..., k.

Für alle j ∈ {1, ..., k} erhält man über w^{j−1} ∈ span(Q) , (5.3) und (5.5) ∥w^{j}∥ = ∥(I − θ^{j}N^{−1}Q)w^{j−1}∥ = ∥(1 − θ^{j})w^{j−1} + θ^{j}(I − N^{−1}Q)w^{j−1}∥ ≤

≤ (1 − θ^{j})∥w^{j−1}∥ + θ^{j}∥(I − N^{−1}Q)w^{j−1}∥ ≤ (1 − θ^{j})∥w^{j−1}∥ + θ^{j}ρ∥w^{j−1}∥ ≤^{θ≥\overline{θ}} (1 − \overline{θ}(1 − ρ))∥w^{j−1}∥

⇒ ∥w^{k}∥ ≤ (1 − \overline{θ}(1 − ρ))^{k}∥w^{0}∥ ≤ (1 − \overline{θ}(1 − ρ))^{k}∥w∥.

Da ∏_{j=1}^{k}(I − θ^{j}QN^{−1})z = Qw^{k} und ∥w∥ ≤ \frac{1}{λ_{min}} ∥z∥

erhält man

∥ ∏_{j=1}^{k}(I − θ^{j}QN^{−1})z∥ = ∥Qw^{k}∥ ≤ ∥Q∥(1 − \overline{θ}(1 − ρ))^{k}∥w∥ ≤ \frac{∥Q∥}{λ_{min}} (1 − \overline{θ}(1 − ρ))^{k}∥z∥.

Setze τ = \frac{∥Q∥}{λ_{min}}, dann folgt Aussage (a).

\underline{zu (b)}

Sei z ∈ IR^{n}, \overline{θ} und \{θ^{j}\}_{j} wie oben. Man zerlege z in z = x + y, x ∈ span(Q) und y ∈ Kern(Q) .

I − N^{−1}Q ist konvergent in Q ⇒ (I − N^{−1}Q)y = y,

also (I − θ^{j}N^{−1}Q)y = (1 − θ^{j})y + θ^{j}(I − N^{−1}Q)y = y, .j = 1, ..., k. Man erhält:

∏_{j=1}^{k}(I − θ^{j}N^{−1}Q)y = y.  (5.6)

Setze z^{0} = x und definiere z^{1}, ..., z^{k} über z^{j} = (I − θ^{j}N^{−1}Q)z^{j−1}, j = 1, ..., k. Induktiv folgt mit (5.3) z^{j} ∈ span(Q) ∀ j. Mittels (5.4) kommt man auf

∥z^{j}∥ = ∥(I − θ^{j}N^{−1}Q)z^{j−1}∥ ≤ (1 − θ^{j})∥z^{j−1}∥ + θ^{j}∥(I − N^{−1}Q)z^{j−1}∥ ≤ (1 − \overline{θ}(1 − ρ))∥z^{j−1}∥, j = 1, ..., k. ⇒

∥ ∏_{j=1}^{k}(I − θ^{j}N^{−1}Q)x∥ = ∥z^{k}∥ ≤ (1 − \overline{θ}(1 − ρ))^{j}∥z^{0}∥ = (1 − \overline{θ}(1 − ρ))^{j}x^{(k)}∥.

Unter Ausnutzung dieser Ungleichung, (5.6) und den Trivialitäten ∥x∥ ≤ ∥z∥ und ∥y∥ ≤ ∥z∥ ergibt sich:

∥ ∏_{j=1}^{k}(I − θ^{j}N^{−1}Q)z∥ ≤ ∥ ∏_{j=1}^{k}(I − θ^{j}N^{−1}Q)x∥ + ∥y∥ ≤ (1 − \overline{θ}(1 − ρ))^{j}∥z∥ + ∥z∥ ≤ 2∥z∥.

Damit ist (b) bewiesen. □

Korollar 5.7 : Sei b ∈ span(Q) y^{(0)} ∈ IR^{m}. Dann konvergiert die Folge y^{(k+1)} = (I − N^{−1}Q)y^{(k)} − N^{−1}b, k = 0, 1, 2, 3, ... geometrisch.

Beweis: Aus Lemma 5.6 (a) folgt:

y^{(k)} = (I − N^{−1}Q)^{k}y^{(0)} − (I − N^{−1}Q)^{k−1}N^{−1}b ⇒

0 ≤ ∥y^{(k)}∥ ≤ τρ^{k}∥y^{(0)}∥ − τρ^{k−1}∥N^{−1}b∥ = ρ^{k}(τ∥y^{(0)}∥ − \frac{τ}{ρ}∥N^{−1}b∥)_{\overrightarrow{k⟶∞}}0, da ρ < 1. □

Im Verlauf dieses Kapitels wird oft mit Untermatrizen gerechnet. Sind I, J ⊆ {1, ..., n} Indexmengen, dann bezeichnet M_{IJ} ∈ IR^{|I|×|J|} die Untermatrix von M, die die Zeilen entsprechend I und die Spalten entsprechend J enthält. Ist J = {1, ..., n}, so schreibt man M_{I·}, bzw. M_{·J} für I = {1, ..., n}.

Korollar 5.8 : (B, C) sei reguläres Splitting der Matrix M. (a) Für jede nicht-leere Indexmenge J ⊆ {1, ..., n} existiert ein ρ_{J} ∈ (0, 1) und ein τ_{J} > 0 mit: ∥(I − M_{JJ}(B_{JJ})^{−1})^{k}z∥ ≤ τ_{J}(ρ_{J})^{k}∥z∥ ∀ k ≥ 1 ∀ z ∈ span(M_{JJ}) . (b) Es existiert ein Δ ≥ 1, so daß für jede nicht-leere Indexmenge J ⊆ {1, ..., n} gilt: ∥(I − (B_{JJ})^{−1}M_{JJ})^{k}z∥ ≤ Δ∥z∥ ∀ k ≥ 1 ∀ z.

Beweis: (B − C) ist positiv definit ⇒ (B_{JJ} − C_{JJ}) ist positiv definit, (B_{JJ}, C_{JJ}) ist folglich reguläres Splitting von M_{JJ}. Setzt man \overline{θ} = 1 in Lemma 5.6, dann ergibt sich die Behauptung. □

Es sei I_{L} = {i ∈ {1, ..., n} : d_{i}^{∗} = 0} mit d^{∗} aus Lemma 5.1. Damit bekommt man: M_{I_{L}}z^{∗} + q_{I_{L}} = 0 ∀ z^{∗} ∈ L. Für alle z^{∗} ∈ L gilt darüberhinaus z^{∗} = [z^{∗} − d^{∗}]^{+}. Für i ∉ I_{L} folgt daraus d_{i}^{∗} > 0 und z_{i}^{∗} = 0 ∀ z^{∗} ∈ L.

Lemma 5.9 : Ist J ⊆ I_{L} nicht-leere Indexmenge, dann gilt span(M_{JJ^{C}}) ⊆ span(M_{JJ}) und q_{J} ∈ span(M_{JJ}) . J^{C} ist das Komplement von J bzgl. {1, ..., n}.

Beweis: Nehme an, daß es ein i ∈ J^{C} gibt, so daß M_{Ji} ∉ span(M_{JJ}) . v_{1}, ..., v_{|J|} ∈ IR^{|J|} sei eine Orthonormalbasis von IR^{|J|} (diese existiert, da M_{JJ} symmetrisch ist). Es seien α_{1}, ..., α_{|J|} ∈ IR so gewählt, das M_{Ji} = α_{1}v_{1} + ... + α_{|J|}v_{|J|}. Wegen M_{Ji} ∉ span(M_{JJ}) existiert ein j ∈ {1, ..., |J|} mit α_{j} ≠ 0 ∧ v_{j} ∈ Kern(M_{JJ}) , v_{j} ≠ 0 ⇒ v_{j}^{t}M_{Ji} = α_{j} ≠ 0 und M_{JJ}v_{j} = 0. Es gibt also Œ ein u ∈ IR^{|J|} mit u^{t}M_{Ji} < 0 ∧ M_{JJ}u = 0.

Definiere x ∈ IR^{n} über x_{j} := \{ \begin{matrix}u_{j} & , j ∈ J \\ 1 & , j = i \\ 0 & , sonst \end{matrix}

⇒ x^{t}Mx = u^{t}M_{JJ}u + 2u^{t}M_{Ji} <^{s.o.} 0.

Damit hat man einen Widerspruch zu M positiv semidefinit. Für alle z^{∗} ∈ L ist q_{I_{L}} = −M_{I_{L}}z^{∗}, folglich q_{J} ∈ span(M_{I_{L}}) . ⇒ q_{J} ∈ span(M_{JJ}) . □

5.2 Konvergenz des Verfahrens z^{(k+1)} = F_{B}(z^{(k)})

(B, C) sei ein reguläres Splitting und F_{B} wie in (5.2) definiert. Es soll gezeigt werden, daß die durch das Verfahren

z^{(k+1)} = F_{B}(z^{(k)}), k = 0, 1, 2, 3, ...  (5.7)

entstehende Vektorenfolge {z^{(k)}}_{k} gegen ein z^{∗} ∈ L konvergiert. Die Wohldefiniertheit der Vektorfolge wurde in Lemma 5.2 (a) gezeigt.

5.2.1 Motivation

Eine Folgerung aus Lemma 5.2 (b) ist

f(z^{(k+1)}) ≤ f(z^{(k)}) − \frac{1}{2}(z^{(k+1)} − z^{(k)})^{t}(B − C)(z^{(k+1)} − z^{(k)}) ≤ f(z^{(k)}) − λ\frac{1}{2}∥z^{(k+1)} − z^{(k)}∥^{2}.

λ > 0 ist der kleinste Eigenwert des symmetrischen Teils von (B − C). Aufgrund der Prämisse (P) ist f für alle k nach unten beschränkt ist. Durch Aufsummieren erhält man:

∑_{k=0}^{∞} ∥z^{(k+1)} − z^{(k)}∥^{2} < ∞.  (5.8)

Lemma 5.10 : {(a) z^{(k+1)} − z^{(k)}_{\overrightarrow{k⟶∞}}0. \atop (b) z^{(k)} − [z^{(k)} − Mz^{(k)} − q]^{+}_{\overrightarrow{k⟶∞}}0.}

Beweis:

(a) ergibt sich aus (5.8).

\underline{zu (b):} Für k ⟶ ∞ ist z^{(k)} − [z^{(k)} − Mz^{(k)} − q]^{+} = (−1)(−z^{(k)} + [z^{(k)} − Mz^{(k)} − q − B(z^{(k+1)} − z^{(k)}) + (z^{(k+1)} − z^{(k)})]^{+}, da [ · ]^{+} lipschitzstetig ist.

Aus z^{(k+1)} = F_{B}(z^{(k)}) folgt d.h. z^{(k+1)} − z^{(k)} = z^{(k)} − [z^{(k)} − Mz^{(k)} − q]^{+}, und mit (a) erhält man die Behauptung. □

Jeder Grenzwert z^{∗} von {z^{(k)}}_{k} erfüllt z^{∗} = [z^{∗} − Mz^{∗} − q]^{+} ⇒ z^{∗} ∈ L. Die Eigenschaft gemäß Lemma 5.10 (a) reicht nicht aus, um zu zeigen, daß {z^{(k)}}_{k} konvergiert, denn es kann z.B. der Fall auftreten, daß z^{(k+1)} − z^{(k)} sich wie \frac{1}{k} verhält, was zur Folge hätte, daß ∥z^{(k)}∥_{\overrightarrow{k⟶∞}}∞.

Φ(z) = min_{z^{∗}∈L} ∥z − z^{∗}∥ sei der Abstand, den z von der Menge L hat ∀ z ∈ X.

Die Folge {Mz^{(k)}}_{k} konvergiert und die Folge {z^{(k)}}_{k} kommt beliebig nahe an L. Diese beiden Aussagen werden bewiesen im

Lemma 5.11 : {(a) Mz^{(k)} + q_{\overrightarrow{k⟶∞}}d^{∗}. \atop (b) Φ(z^{(k)})_{\overrightarrow{k⟶∞}}0.}

Beweis: Beweis von Teil (a): Wir brauchen 3 Vorbemerkungen, die nicht schwer zu beweisen sind (siehe [Pang ]) :

(1) Ist P ⊆ IR^{n} ein Polyhedron, A ∈ IR^{m×n}, dann hat jeder Häufungspunkt der Folge Az^{(k)} mit {z^{(k)}}_{k} ⊂ P die Darstellung Aẑ mit ẑ ∈ P.

(2) Ist C ∈ IR^{n×n} symmetrisch und positiv semidefinit, A ∈ IR^{m×n}.

Dann gilt: \{ACA^{t}z^{(k)}\}_{\overrightarrow{k⟶∞}}ACA^{t}ẑ ⇔ A^{t}z^{(k)}_{\overrightarrow{k⟶∞}}A^{t}ẑ.

(3) M ∈ IR^{n×n} sei symmetrisch und positiv semi-definit.

I ⊆ {1, ..., n} sei eine Indexmenge.

Dann gilt M_{II}z_{I}^{(k)}_{\overrightarrow{k⟶∞}}M_{II} ẑ_{I} ⇔ M_{·I}z_{I}^{(k)}_{\overrightarrow{k⟶∞}}M_{·I} ẑ_{I} .

Mit diesen vorbereitenden Bemerkungen beweisen wir, daß aus der Prämisse (P) (a) folgt. Nach Lemma 5.2 (a) ist {z^{(k)}}_{k} eindeutig definiert.

{Mz^{(k)}}_{k} hat einen Häufungspunkt Ist {z^{(k)}}_{k} beschränkt, dann ist diese Behauptung klar. Nehme also an, {z^{(k)}}_{k} ist unbeschränkt ⇒ ∃ nicht-leere Indexmenge I ⊆ {1, ..., n}

und eine Teilfolge \{z^{(k_{j})}\}_{j} mit z_{i}^{(k_{j})}_{\overrightarrow{j⟶∞}}∞ ∀ i ∈ I und \{z_{i}^{(k_{j})}\}_{j} ist beschränkt ∀ i ∉ I.

Sei J = I^{C} ⇒ M_{II}z_{I}^{(k_{j})} = − [(q + C(z^{(k_{j}−1)} − z^{(k_{j})}))_{I} + M_{IJ}z_{J}^{(k_{j})}] ∀ k_{j} groß genug. Dies schließt man aus der Komplementaritätsbedingung. Also ist {M_{II}z_{I}^{(k_{j})}}_{j} beschränkt. Mit Vorbemerkung (3) folgt: {M_{·I}z_{I}^{(k_{j})}}_{j} hat einen Häufungspunkt. Da {Mz^{(k_{j})} = M_{·I}z^{(k_{j})} + M_{·J}z^{(k_{j})}}_{j} und {z_{J}^{(k_{j})}}_{j} beschränkt sind, ergibt sich: {Mz^{(k_{j})}}_{j} hat einen Häufungspunkt ⇒ {Mz^{(k)}}_{k} hat einen Häufungspunkt z^{HP}.

\{Mz^{(k_{i})}\}_{i} sei eine Teilfolge mit Mz^{(k_{i})}_{\overrightarrow{i⟶∞}}z^{HP} ⇒ ∃ eine (evtl. leere) Indexmenge I ⊆ \{1, ..., n\}

und eine Teilfolge {z^{(l_{i})}}_{i} von {z^{(k_{i})}}_{i} mit z_{j}^{(l_{i})} ⟶ ∞ ∀ j ∈ I und {z_{j}^{(l_{i})}}_{i} ist für alle j ∈ I beschränkt.

Mit J = I^{C} erhält man wieder

M_{II}z_{I}^{(l_{j})} = − [(q + C(z^{(l_{j}−1)} − z^{(l_{j})}))_{I} + M_{IJ}z_{J}^{(l_{j})}] ,  (5.9)

wenn l_{j} groß genug ist.

Nach Vorbemerkung (3) hat {M_{·I}z_{I}^{(l_{i})}}_{i} einen Häufungspunkt der Gestalt M_{·I} ẑ_{I} mit ẑ_{I} ≥ 0. Œ ist M_{·I}z_{I}^{(l_{i})}_{\overrightarrow{i⟶∞}}M_{·I} ẑ_{I}.

{z_{J}^{(l_{i})}}_{i} ist beschränkt, also konvergiert {z_{j}^{(l_{i})}}_{i} gegen ein ẑ_{j} ≥ 0 ∀ j ∈ J ⇒ z^{HP} = Mẑ.

\underline{ẑ löst das LKP (q, M)} ẑ ≥ 0 ist klar. Für l_{i} ⟶ ∞ ergibt sich über (5.9): (q + Mẑ)_{I} = 0, da nach Lemma 5.10 (a) z^{(l_{i}−1)} − z^{(l_{i})} ⟶ 0. Außerdem gilt für alle l_{i}: 0 ≤ (q + Cz^{(l_{i}−1)} + Bz^{(l_{i})})_{J} = (q + C(z^{(l_{i}−1)} − z^{(l_{i})}) + Mz^{(l_{i})})_{J}. Für l_{i} ⟶ ∞ hat man also: z_{J}^{(l_{i})} ⟶ ẑ_{J} und Mz^{(l_{i})} ⟶ Mẑ ⇒ (q + Mẑ)_{J} ≥ 0 und ẑ_{J}^{t}(q + Mẑ)_{J} = 0. Also löst ẑ das LKP (q, M) . Zusammenfassend stellt man fest: Ist z^{HP} ein Häufungspunkt der Folge {Mz^{(k)}}_{k}, so existiert eine Lösung ẑ des LKP (q, M) mit z^{HP} = Mẑ. Nach Lemma 5.1 ist z^{HP} eindeutig bestimmt ⇒ Mz^{(k)}_{\overrightarrow{k⟶∞}}Mẑ, mit ẑ ∈ L.

Zum Beweis von Teil (b) benutzt man (a) und die obere Lipschitzeigenschaft von L. Näheres dazu in [Luo ]. □

Korollar 5.12 : Zu jedem δ > 0 gibt es ein k_{δ} ∈ IN mit: Φ(z^{(k)}) ≤ δ ∀ k ≥ k_{δ}  (5.10) ∥z^{(k+1)} − z^{(k)}∥ ≤ δ ∀ k ≥ k_{δ}  (5.11) z_{i}^{(k)} ≤ δ ∀ k ≥ k_{δ} ∀ i ∉ I_{L}.  (5.12)

Beweis: (5.10) und (5.11) sind eine Folgerung aus den Lemmata 5.10 (a) und 5.11 (b). (5.12) ergibt sich aus (5.10). □ Es lassen sich also folgende Aussagen machen:

• {z^{(k)}}_{k} nähert sich mit beliebig kleinem Abstand der Lösungsmenge L; da jedoch nicht bekannt ist, ob {z^{(k)}}_{k} beschränkt ist, können wir daraus nicht auf die Konvergenz der Folge schließen.

• Für i ∉ I_{L} folgt d_{i}^{∗} > 0 und z_{i}^{∗} = 0. Mit Lemma 5.10 (b) und 5.11 (a) ergibt sich z_{i}^{(k)}_{\overrightarrow{k⟶∞}}0.

{z_{i}^{(k)}}_{k} ist monoton fallend in {z_{i} ≥ 0}, also existiert ein N ∈ IN : z_{i}^{(k)} = 0 ∀ k ≥ N.

• Gilt die Implikation: {i ∈ I_{L} ⇒ z_{i}^{∗} > 0} ( man sagt auch “ z ist nicht ausgeartet” ), dann läuft die LKP-Lösung hinaus auf das unrestriktive Problem: Suche ein z ∈ IR^{n} mit: M_{I_{L}}z + q_{I_{L}} ≥ 0 und z^{t}(M_{I_{L}}z + q_{I_{L}}) = 0.

z^{(k+1)} = F_{B}(z^{(k)}) reduziert sich in diesem Fall auf

Bz^{(k+1)} + Cz^{(k)} + q = 0 ⇔ z^{(k+1)} = B^{−1}(Cz^{(k)} + q) = (I − B^{−1}M)z^{(k)} − B^{−1}q.

Nach Korollar 5.7 konvergiert dieses Verfahren geometrisch.

Problematisch sind also die Koordinaten z_{i}^{(k)}, i ∈ I_{L}, die um die untere Schranke 0 herum beschränkt sind, die 0 jedoch nicht annehmen.

Ein sinnvoller Konvergenzbeweis muß folglich die Koordinaten der Iterationswerte getrennt untersuchen. Dies führt auf den Ansatz der

5.2.2 Separation der Koordinaten

Definiere die Konstante

β = 1 + max_{J⊆I_{L}} \sqrt{|J^{C}|} \{(\frac{τ_{J}∥(B_{JJ})^{−1}∥ ∥M_{JJ}∥}{1 − ρ_{J}} + Δ + 1) ∥(B_{JJ})^{−1}B_{JJ^{C}}∥ + \frac{τ_{J}∥(B_{JJ})^{−1}∥ ∥M_{JJ}∥}{1 − ρ_{J}}\} .

Die verwendeten Konstanten sind bekannt aus Korollar 5.8.

Lemma 5.13 : (Separation der Koordinaten) Sei J ⊆ I_{L}. Gibt es s, t ∈ IN, s ≥ t ≥ 0, so daß z_{i}^{(k)} > 0 ∀ i ∈ J ∀ k = t, t + 1, ..., s, so gilt für jedes z^{∗} ∈ L: ∥z_{J}^{(s)} − z_{J}^{∗}∥ ≤ Δ∥z_{J}^{(t)} − z_{J}^{∗}∥ + β max_{k∈{t,...,s}} ∥z_{J}^{(k)} − z_{J}^{∗}∥_{∞}.

Beweis: Dieser Beweis wird sehr genau von [Luo ] geführt. Theoretisch ist die Beweisführung nicht schwierig, allerdings sind in geschickter Weise Gleichungen umzuformen und Hilfsvektoren zu definieren. Wir beschränken uns an dieser Stelle auf die Darstellung der wichtigsten Schritte. Für die genauen Details ziehe man [Luo ] zu Rate.

Für s = t ist die Aussage klar, da Δ ≥ 1. Man untersuche also den Fall s > t.

\{ \begin{matrix}z_{i}^{(k)} > 0 & ∀ i ∈ J, k = t + 1, ..., s \\ z^{(k+1)} = [z^{(k+1)} − (Bz^{(k+1)} + Cz^{(k)} + q)]^{+} & ∀ k \end{matrix} \}

⇒ B_{J·}z^{(k+1)} + C_{J·}z^{(k)} + q_{J} = 0, k = t, ..., s − 1.

Da M_{J·} = B_{J·} + C_{J·} gilt, ist dies äquivalent zu B_{J·}(z^{(k+1)} − z^{(k)}) + M_{J·}z^{(k)} + q_{J} = 0, k = t, ..., s − 1. Benutzt man Lemma 5.1, so folgt wegen J ⊆ I_{L} : M_{J·}z^{∗} + q_{J} = 0. Subtrahiere beide Gleichungen und multipliziere sie mit (B_{JJ})^{−1}:

(B_{JJ})^{−1}B_{J·}(z^{(k+1)} − z^{(k)}) + (B_{JJ})^{−1}M_{J·}(z^{(k)} − z^{∗}) = 0.

Über die Definition G := I − (B_{JJ})^{−1}M_{JJ} bekommt man daraus:

z_{J}^{(s)} − z_{J}^{∗} = G^{h}(z_{J}^{(t)} − z_{J}^{∗}) − ∑_{k=0}^{h−1} G^{h−k−1}(B_{JJ})^{−1}M_{JJ^{C}}(z_{J^{C}}^{(t+k)} − z_{J^{C}}^{∗}) −

− ∑_{k=0}^{h−1} G^{h−k−1}(B_{JJ})^{−1}B_{JJ^{C}}(z_{J^{C}}^{(t+k+1)} − z_{J^{C}}^{(t+k)}), mit h = s − t.

Um diese Formel abzukürzen, definiert man

y^{(k)} = (B_{JJ})^{−1}B_{JJ^{C}}(z_{J^{C}}^{(t+k)} − z_{J^{C}}^{∗}), ∀ k = 0, ..., h

und H = I − M_{JJ}(B_{JJ})^{−1} ⇒ G = (B_{JJ})^{−1}HB_{JJ} ⇒ G^{h−k−1} = (B_{JJ})^{−1}H^{h−k−1}B_{JJ} ∀ k.

Die obige Gleichung verkürzt sich damit auf z_{J}^{(s)} − z_{J}^{∗} = G^{h}(z^{(t)} − z_{J}^{∗}) −

− ∑_{k=0}^{h−1}(B_{JJ})^{−1}H^{h−k−1}M_{JJ^{C}}(z_{J^{C}}^{(t+k)} − z_{J^{C}}^{∗}) + ∑_{k=0}^{h−1}(B_{JJ})^{−1}H^{k−h−1}M_{JJ}y^{(k)} − y^{(h)} + G^{h−1}y^{(0)}.

Es ist J ⊆ I_{L}, daher ergibt sich aus Korollar 5.7 und mit den dort definierten Konstanten: ∥H^{h−k−1}z∥ ≤ τ_{J}(ρ_{J})^{h−k−1}∥z∥ ∀ z ∈ span(M_{JJ}) und ∥G^{h}z∥ ≤ Δ∥z∥ ∀ z. Mit Hilfe der Aussage span(M_{JJ^{C}}) ⊆ span(M_{JJ}) aus Lemma 5.9 errechnet man folgende Abschätzung

∥z_{J}^{(s)} − z_{J}^{∗}∥ ≤ Δ∥z_{J}^{(t)} − z_{J}^{∗}∥ + τ_{J}∥(B_{JJ})^{−1}∥ ∥M_{JJ^{C}}∥(1 − ρ_{J})^{−1} max_{r∈{t,...,s−1}} ∥z_{J^{C}}^{(r)} − z_{J^{C}}^{∗}∥ +

+τ_{J}∥(B_{JJ})^{−1}∥ ∥M_{JJ}∥(1 − ρ_{J})^{−1} max_{r∈{1,...,h−1}} ∥y^{(k)}∥ + ∥y^{(h)}∥ + Δ∥y^{(0)}∥.

Laut Vorgabe gilt: ∥y^{(k)}∥ ≤ ∥(B_{JJ})^{−1}B_{JJ^{C}}∥ ∥z_{J}^{(t+k)} − z_{J^{C}}^{∗}∥ ∀ k. Also: ∥z_{J}^{(s)} − z_{J}^{∗}∥ ≤ Δ∥z_{J}^{(t)} − z_{J}^{∗}∥ +

+ max_{r∈\{1,..,s\}} ∥z_{J}^{(r)}−z_{J}^{∗}∥_{∞} \{\frac{τ_{J}∥(B_{JJ})^{−1}∥ ∥M_{JJ^{C}}∥}{1 − ρ_{J}} + (\frac{τ_{J}∥(B_{JJ})^{−1}∥ ∥M_{JJ}∥}{1 − ρ_{J}} + Δ + 1) ∥(B_{JJ})^{−1}∥ ∥B_{JJ^{C}}∥\} =

= Δ∥z_{J}^{(t)} − z_{J}^{∗}∥ + max_{k∈\{1,..,s\}} ∥z_{J}^{(k)} − z_{J}^{∗}∥_{∞} \frac{β − 1}{max_{J⊆I_{L}} \sqrt{|J^{C}|}} ≤ Δ∥z_{J}^{(t)} − z_{J}^{∗}∥ + max_{k∈\{1,..,s\}} ∥z_{J}^{(k)} − z_{J}^{∗}∥_{∞}β. □

Wir benötigen für die sich anschließenden Beweise noch folgende weitere Konstanten: σ_{0} = 1 und σ_{r} = Δ + 3 + β + (β + 1)σ_{r−1}, r = 1, ..., n.

Lemma 5.14 : Setze δ > 0. k_{δ} ∈ IN sei ein Skalar, der (5.10) und (5.11) erfüllt. Existieren r ∈ {1, ..., n}, eine nicht-leere Indexmenge J ⊆ I_{L}, sowie k, k′ ∈ IN mit k′ > k ≥ k_{δ}, die z_{i}^{(k)} > σ_{r}δ ∀ i ∈ J  (5.13) z_{i}^{(t)} ≤ σ_{r−1}δ ∀ i ∉ J ∀ t = k, k + 1, ..., k′ − 1  (5.14) erfüllen, dann gilt: (a) z_{i}^{(k′)} > σ_{r−1}δ ∀ i ∈ J. (b) Es gibt ein z^{∗} ∈ L : ∥z^{(t)} − z^{∗}∥_{∞} ≤ σ_{r}δ ∀ t = k, k + 1, ..., k′ − 1.

Beweis: z^{∗} ∈ L sei so gewählt, daß Φ(z^{(k)}) = ∥z^{(k)} − z^{∗}∥. Aus (5.10) folgt dann

∥z^{(k)} − z^{∗}∥ < δ.  (5.15)

Mit (5.14) ergibt sich für alle i ∉ J :

z_{i}^{∗} ≤ z_{i}^{(k)} + ∥z^{(k)} − z^{∗}∥ ≤ σ_{r−1}δ + ∥z^{(k)} − z^{∗}∥ ⟹^{(5.15)} 0 ≤ z_{i}^{∗} ≤ σ_{r−1}δ + δ.

Für t = k, ..., k′ − 1 ist (5.14) vorausgesetzt: 0 ≤ z^{(t)} ≤ σ_{r−1}δ ⇒

|z_{i}^{(t)} − z_{i}^{∗}| ≤ σ_{r−1}δ + δ ∀ i ∉ J ∀ t = k, ..., k′ − 1.  (5.16)

Induktiv soll als nächstes gezeigt werden, daß für t = k, k + 1, ..., k′ − 1 gilt:

z_{i}^{(t)} > σ_{r−1}δ + δ ∀ i ∈ J.  (5.17)

\underline{t = k} (5.17) erhält man für diesen Fall unmittelbar wegen (5.13) und σ_{r} ≥ σ_{r−1} + 1. Induktionsschritt: Vorausgesetzt sei, (5.17) ist für t = k, k + 1, ..., s erfüllt, wobei s ∈ {t, t + 1, ..., t′ − 2}. Zu zeigen: (5.17) gilt auch für t = s + 1. Nach (5.17) ist z_{i}^{(t)} > 0 ∀ i ∈ J ∀ t = k, ..., s

⟹^{Lemma 5.13} ∥z_{J}^{(s)} − z_{J}^{∗}∥ ≤ Δ∥z_{J}^{(k)} − z_{J}^{∗}∥ + β max_{t∈{k,...,s}} ∥z_{J}^{(t)} − z_{J}^{∗}∥_{∞}

⟹^{(5.15), (5.16)}

∥z_{J}^{(s)} − z_{J}^{∗}∥ ≤ Δδ + β(σ_{r−1}δ + δ).  (5.18)

Ist i ∈ J, so hat man: z_{i}^{(s+1)} ≥ z_{i}^{(k)} − ∥z_{J}^{(k)} − z_{J}^{(s+1)}∥ ≥ z_{i}^{(k)} − (∥z_{J}^{(k)} − z_{J}^{∗}∥ + ∥z_{J}^{∗} − z_{J}^{(s)}∥ + ∥z_{J}^{(s)} − z_{J}^{(s+1)}∥). Mit (5.11), (5.13) und (5.15) folgert man:

z_{i}^{(s+1)} > σ_{r}δ − (δ + ∥z_{J}^{∗} − z_{J}^{(s)}∥ + δ) ≥^{(5.18)} σ_{r}δ − (δ + (Δδ + βσ_{r−1}δ + βδ) = σ_{r−1}δ + δ.

Soweit der Induktionsbeweis. (5.17) gilt also ∀ t = k, ..., k′ − 1. ⇒ (5.18) gilt für alle s = k, ..., k′ − 1. Dies kombiniert mit (5.16) führt auf ∥z^{(t)} − z^{∗}∥_{∞} ≤ (Δ + βσ_{r−1} + β)δ. Ausgenutzt wurde hier: β > 1 und ∥z∥_{∞} ≤ ∥z∥. Mit Δ + βσ_{r−1} + β ≤ σ_{r} folgt Behauptung (b). Aus (5.18) erhält man für alle i ∈ J mit t = k′ − 1 :

z_{i}^{(k′)} ≥ z_{i}^{(k′−1)} − ∥z^{(k′−1)} − z^{(k′)}∥ > σ_{r−1}δ + δ − ∥z^{(k′−1)} − z^{(k′)}∥ ≥^{(5.11)} σ_{r−1}δ + 2δ > σ_{r−1}δ.

Dies ist Behauptung (a). □

Lemma 5.15 : Setze δ > 0. Der Skalar k_{δ} ∈ IN erfülle (5.10) und (5.11). Existiert ein r ∈ \{1, ..., n\}, so daß für die nicht-leere Indexmenge J ⊆ I_{L} |J| ≥ |I_{L}| − r + 1 gilt und außerdem ein t ≥ k_{δ} existiert mit z_{i}^{(t)} > σ_{r}δ ∀ i ∈ J und z_{i}^{(t)} ≤ σ_{r−1}δ ∀ i ∉ J, dann gibt es ein z^{∗} ∈ L und ein \overline{t} ≥ t mit ∥z^{(k)} − z^{∗}∥_{∞} ≤ σ_{r}δ ∀ k ≥ \overline{t}. (5.19)

Beweis: Durch Induktion über r: \underline{r = 1} J = I_{L}. Bedingung (5.14) ist wegen (5.12) erfüllt, wenn t′ > t, also

∃ z^{∗} ∈ L : ∥z^{(k)} − z^{∗}∥ ≤ σ_{r}δ ∀ k ≥ t. Nehme nun an, die Behauptung gilt für r = 1, 2, ..., h − 1 für ein h ≥ 2.

J ⊆ I_{L} und t ∈ IN seien beliebig, aber fest mit |J| ≥ |I_{L}| − h + 1

und z_{i}^{(t)} > σ_{h}δ ∀ i ∈ J ∧ z_{i}^{(t)} ≤ σ_{h−1}δ ∀ i ∉ J. Unterscheide 2 Fälle: \underline{1.Fall:} z_{i}^{(k)} ≤ σ_{h−1}δ ∀ i ∉ J ∀ k ≥ t

z_{i}^{(t)} > σ_{h}δ ∀ i ∈ J ist vorausgesetzt. ⟹^{Lemma 5.11(b)} ∃ z^{∗} ∈ L : ∥z^{(k)} − z^{∗}∥ ≤ σ_{h}δ ∀ k ≥ t.

Folglich ist (5.18) erfüllt für k = h. \underline{2.Fall:} Es existiert ein k > t und ein i ∉ J mit z_{i}^{(k)} > σ_{h−1}δ.

Sei t′ = min{k : k > t ∧ z_{i}^{(k)} > σ_{h−1}δ, i ∉ J}

⇒ z_{i}^{(k)} ≤ σ_{h−1}δ ∀ i ∉ J ∀ k = t, t + 1, ..., t′ − 1.

nach Voraussetzung ist z_{i}^{(t)} > σ_{h}δ ∀ i ∈ J.. Mit Lemma 5.11 (a) folgert man:

z_{i}^{(t′)} > σ_{h−1}δ ∀ i ∈ J.  (5.20)

Man betrachte die h + 1 Intervalle

[0, σ_{0}δ), (σ_{0}δ, σ_{1}δ], (σ_{1}δ, σ_{2}δ], ..., (σ_{h−2}δ, σ_{h−1}δ], (σ_{h−1}δ, ∞).

Aus (5.20) und z_{i}^{(t′)} > σ_{h−1}δ für ein i ∉ J ergibt sich, daß das Intervall (σ_{h−1}δ, ∞) |J| + 1 Elemente der Menge {z_{1}^{(t′)}, ..., z_{n}^{(t′)}} enthält.

Aus (5.12) und σ_{0} = 1 ergibt sich, daß das Intervall [0, σ_{0}δ) n − |I_{L}| Elemente derselben Menge enthält. Da |J| ≥ |I_{L}| − h + 1 vorausgesetzt wurde, enthalten die restlichen h − 1 Intervalle höchstens h − 2 Elemente der Menge. Es gibt folglich ein j ∈ {1, 2, ..., h − 1} mit z_{i}^{(t′)} ∉ (σ_{j−1}δ, σ_{j}δ] ∀ i. Man spricht hier auch vom Ablagefachprinzip. h′ sei größter Index j mit dieser Eigenschaft, d.h. (σ_{h′}δ, ∞) enthält mindestens |J| + h − h′ Elemente der Menge {z_{1}^{(t′)}, ..., z_{n}^{(t′)}}. J′ sei die Indexmenge dieser Elemente: J′ = {i : z_{i}^{(t′)} > σ_{h′−1}δ} ⇒ z_{i}^{(t′)} > σ_{h′}δ ∀ i ∈ J′ und z_{i}^{(t′)} ≤ σ_{h′−1}δ ∀ i ∉ J′ und |J′| ≥ |J| + h − h′ ≥ |I_{L}| + 1 − h′.

⟹^{(5.12)} J′ ⊆ I_{L}.

Es ist h′ < h. Wendet man daher die Induktionsvoraussetzung an auf h′, t′ und J′, so folgt, daß es ein z^{∗} ∈ L und ein \overline{t} ≥ t′ gibt mit ∥z^{(k)} − z^{∗}∥_{∞} ≤ σ_{h′}δ ∀ k ≥ \overline{t}. Damit ist die Behauptung auch für k = h gezeigt. □

Die Lemmata 5.13-5.15 bilden die Grundlage für den Beweis der Schlüsselaussage in diesem Abschnitt.

Lemma 5.16 : Zu jedem δ > 0 existiert ein z^{∗} ∈ L und ein k̂ ∈ IN mit ∥z^{(k)} − z^{∗}∥ < σ_{n}δ + δ ∀ k ≥ k̂.

Beweis: k_{δ} ∈ IN sei ein Skalar, der (5.10) und (5.11) erfüllt. Für \overline{k} ∈ IN mit \overline{k} ≥ k_{δ} können 2 Fälle vorkommen: \underline{1.Fall:} z_{i}^{(k)} ≤ σ_{n}δ ∀ i ∀ k ≥ \overline{k}. \underline{2.Fall:} Es existiert ein t ≥ \overline{k} und ein i mit z_{i}^{(t)} > σ_{n}δ. \underline{zum 1.Fall:}

Sei z^{∗} ∈ L mit Φ(z^{(\overline{k})}) = ∥z^{(\overline{k})} − z^{∗}∥ ⟹^{(5.14)} 0 ≤ z_{i}^{∗} ≤ z_{i}^{(\overline{k})} + ∥z^{(\overline{k})} − z^{∗}∥ ≤ σ_{n}δ + δ.

Mit 0 ≤ z_{i}^{(k)} ≤ σ_{n}δ ∀ i ∀ k ≥ \overline{k} erhält man ∥z^{(k)} − z^{∗}∥_{∞} ≤ σ_{n}δ + δ ∀ k ≥ \overline{k}. \underline{zum 2.Fall:} Eines der n Intervalle (σ_{0}δ, σ_{1}δ], (σ_{1}δ, σ_{2}δ], ..., (σ_{n−1}δ, σ_{n}δ] enthält kein Element aus der Menge {z_{1}^{(t)}, ..., z_{n}^{(t)}}. j ∈ {1, ..., n} sei der Index mit z_{i}^{(t)} ∉ (σ_{j−1}δ, σ_{j}δ] ∀ i. r sei maximaler Index dieser j. Weiter sei J = {i : z_{i}^{(t)} > σ_{r}δ}. Es folgt |J| ≥ n − r + 1 sowie z_{i}^{(t)} > σ_{r}δ ∀ i ∈ J ∧ z_{i}^{(t)} ≤ σ_{r−1}δ ∀ i ∉ J. Aus (5.16) und σ_{r} ≥ 1 folgt J ⊆ I_{L}. Also erfüllen r, J und t die Bedingungen in Lemma 5.15, folglich existiert ein z^{∗} ∈ L und ein \overline{t} ≥ t mit ∥z^{(k)} − z^{∗}∥_{∞} ≤ σ_{r}δ ∀ k ≥ \overline{t}. Aus σ_{r} ≤ σ_{n} folgt die Behauptung, wenn man k̂ = \overline{t} setzt. □

Als Conclusio dieses Abschnitts notieren wir:

Satz 5.17 : Das Verfahren (5.7) ist wohldefiniert und erzeugt eine Folge von Vektoren {z^{(k)}}_{k}, die gegen ein Element aus L konvergiert.

Beweis: Die Wohldefiniertheit ist eine Konsequenz aus Lemma 5.2(a).

Sei ε > 0 ⟹^{Lemma 5.16} ∃ z^{∗} ∈ L ∧ k̂ ∈ IN : ∥z^{(k)} − z^{∗}∥ < \frac{ε}{2} ∀ k ≥ k̂.

⇒ Für alle k_{1}, k_{2} > k̂ :

∥z^{(k_{1})} − z^{(k_{2})}∥_{∞} ≤ ∥z^{(k_{1})} − z^{∗}∥ + ∥z^{∗} − z^{(k_{2})}∥ ≤ \frac{ε}{2} + \frac{ε}{2} = ε.

Daraus folgt, daß {z^{(k)}}_{k} eine Cauchy-Folge und damit konvergent ist. □

5.3 Die SOR-Erweiterung des Konvergenzbeweises

Um die in Kapitel 4 erläuterten Marktgleichgewichtsprobleme zu lösen, verwenden wir das projezierte SOR-Verfahren und das erweiterte Hildreth-Verfahren. Beide Methoden lassen sich wie (5.7) formulieren. Der oben geführte Konvergenzbeweis ist also übertragbar. (B, C) sei ein Splitting von M und die Relaxationsparameter ω_{1} und ω_{2} seien so gewählt, daß:

\begin{aligned}ω_{2} & > & 0 \\ x^{t}[B − C + (1 − ω_{2})M]x & > & 0 ∀ x ∈ IR^{n} \\ 0 < & ω_{1} & ≤ min\{1, ω_{2}\}. \end{aligned}  (5.21)

D ∈ IR^{n×n} sei eine Diagonalmatrix mit positiven Diagonalelementen. Definiere die Vektorenfolge {z^{(k)}}_{k} über

z^{(0)} ∈ X und z^{(k+1)} = (1 − ω^{(k)})z^{(k)} + ω^{(k)}x^{(k)}.  (5.22)

Dabei ist y = x^{(k)} eine Lösung von

y = [y − D(By + Cz^{(k)} + q)]^{+}  (5.23)

und ω^{(k)} ∈ [ω_{1}, ω_{2}] sei so gewählt, daß z^{(k+1)} ∈ X. Ist (B, C) ein reguläres Splitting, und erfüllt ω_{2} die Abschätzung 0 < ω_{2} < 1 + \frac{1}{∥Q^{− \frac{1}{2}} MQ^{− \frac{1}{2}} ∥},

dann lassen sich ω_{1} und ω_{2} wählen wie in (5.21). Q ist dabei der symmetrische Teil von M. (5.7) erhält man aus (5.22)/(5.23), indem man ω_{1} = ω_{2} = 1 und D = I setzt. Die Multiplikation mit der Matrix D ändert nichts an den Iterationsvektoren im Verfahren (5.7). b_{i} bzw. c_{i} bezeichne die i-te Zeile von B bzw. C. Definiere τ_{i} := d_{ii}(b_{i}^{t}y + c_{i}^{t}z^{(k)} + q_{i}). Damit ergibt sich die Äquivalenzkette:

y löst (5.23) ⟺

y ∈ X ∧ {y_{i} = 0 ⇒ τ_{i} > 0 ∀ i = 1, ..., n}

⟺^{d_{ii}>0}

y ∈ X ∧ \{y_{i} = 0 ⇒ \frac{τ_{i}}{d_{ii}} > 0 ∀ i = 1, ..., n\}

⟺

y = [y − (By + Cz^{(k)} + q)]^{+}, bzw. y = F_{B}(z^{(k)}).

Da der Iterationswert also unabhängig von D ist, kann D beliebig verändert werden. Durch Modifikation der in Abschnitt 5.2.2 verwendeten Argumente beweist man die Konvergenz des Verfahrens (5.22)/(5.23).

Satz 5.18 : (B, C) sei ein Splitting von M. Für die Skalare ω_{1} und ω_{2} gelte (5.21). D sei Diago- nalmatrix mit positiven Diagonalelementen, z^{(0)} ∈ X. Dann konvergiert die Vekto- renfolge {z^{(k)}}_{k}, die durch (5.22)/(5.23) entsteht.

Beweis: Der Beweis dieses Satzes verläuft mit Ausnahme von geringfügigen Änderungen genauso wie der Konvergenzbeweis des Verfahrens (5.7). Aus diesem Grund zitieren wir Abschnitt 5.2.

Es ist 2B = (B − C + (1 − ω_{2})M) + ω_{2}M und ω_{2}M ist positiv semi-definit. ⟹^{(5.21)} B ist positiv definit. Nach Lemma 5.2 (a) ist F_{B} eine wohldefinierte Punkt-zu-Punktabbildung. Für k = 0, 1, 2, 3, ... sei x^{(k)} Lösung von (5.23). Dann läßt sich schreiben: x^{(k)} = F_{B}(z^{(k)}), k = 0, 1, 2, 3, ... also ist {x^{(k)}}_{k} wohldefiniert ⟹^{(5.22)} {z^{(k)}}_{k} ist wohldefiniert.

Der Beweis, daß {z^{(k)}}_{k} konvergiert, kann genauso geführt werden, wie der Beweis zu Satz 5.17 : Setzt man in Lemma 5.6 Q = M_{JJ}, N = \frac{1}{ω} B_{JJ}, θ = \frac{ω_{1}}{ω_{2}} ,

dann gilt folgende allgemeinere Version von Korollar 5.8:

(a) Sei J ⊆ {1, ..., n} eine nicht-leere Indexmenge, dann existieren ρ_{J} ∈ (0, 1) und τ_{J} > 0 mit

∥ ∏_{h=s}^{s+r−1}(I − ω^{h}M_{JJ}(B_{JJ})^{−1})z∥ ≤ τ_{J}(ρ_{J})^{r}∥z∥ ∀ r ≥ 1 ∀ s ≥ 0 ∀ z ∈ span(M_{JJ}) .

(b) ∃ Δ ≥ 1, so daß für jede nicht-leere Indexmenge J ⊆ {1, ..., n} gilt:

∥ ∏_{h=s}^{s+r−1}(I − ω^{h}(B_{JJ})^{−1})M_{JJ}z∥ ≤ Δ∥z∥ ∀ r ≥ 1 ∀ s ≥ 0 ∀ z.

Mit (a) und (b) beweist man eine modifizierte Version von Satz 5.13, in der anstelle von “ z_{i}^{(k)} > 0 ∀ i ∈ J ∀ k = t, ..., s ” steht “ x_{i}^{(k)} > 0 ∀ i ∈ J ∀ k = t, ..., s ”, und mit einer neuen Konstanten β, die nicht nur von M und B, sondern auch von ω_{1} und ω_{2} abhängt. Aus (5.21), (5.22), (5.23) und ω^{(k)} ∈ [ω_{1}, ω_{2}] läßt sich f(z^{(k+1)}) ≤ f(z^{(k)}) − λ\frac{1}{2}∥z^{(k+1)} − z^{(k)}∥ herleiten; λ ist hier der kleinste Eigenwert des symmetrischen Teils von \frac{1}{ω_{1}} [(1 − ω_{2})M + B − C]. Lemma 5.10 und 5.11 gelten also ebenfalls.

⟹^{(5.22)} z^{(k+1)} − x^{(k)} = (1 − \frac{1}{ω^{(k)}} ) (z^{(k+1)} − z^{(k)}) ∀ k

⇒ ∥z^{(k+1)} − x^{(k)}∥ ≤ max \{1, \frac{1}{ω_{1}} \} ∥z^{(k+1)} − z^{(k)}∥.

Definiere σ_{0} = max\{1, \frac{1}{ω^{(k)}} \} und die übrigen σ_{i} wie gehabt, dann beweist man Lemma 5.14 ebenfalls analog (verwendet werden beim Beweis die Lemmata 5.10, 5.11 und 5.13). Lemma 5.15 und 5.16 werden bewiesen mit Lemma 5.10, 5.11 und 5.14. □

Mit Satz 5.18 läßt sich nun die Konvergenz der Verfahren aus Kapitel 3 beweisen.

M ∈ IR^{n×n} habe positive Diagonalelemente, M := (m_{ij})_{i,j}. Man betrachte das Verfahren

z_{i}^{(k+1)} = [z_{i}^{(k)} − \frac{ω}{m_{ii}} (∑_{j<i} m_{ij}z_{j}^{(k+1)} + ∑_{j≥i} m_{ij}z_{j}^{(k)} + q_{i})]^{+}, i = 1, ..., n.  (5.24)

ω ∈ (0, 2) ist der Relaxationsparameter. (5.24) ist eine spezielle Form von (5.22)/(5.23), was deutlich wird, wenn man folgende Zuweisungen trifft: ω_{1} = ω_{2} = 1, D_{M} ∈ IR^{n×n} sei die Diagonale von M, L sei linke untere △-Matrix von M.

(B, C) und D seien:

B = \frac{1}{ω} (D_{M} + L), C = (1 − \frac{1}{ω})D_{M} + L^{t}, D = ωD_{M}^{−1}.

B − C = (2\frac{1}{ω} − 1)D_{M} + L − L^{t} ist positiv definit ∀ ω ∈ (0, 2). Nach Satz 5.18 ist (5.24) konvergent.

(5.24) ist nun nichts anderes, als der PSOR-Algorithmus und für ω = 1 erhält man das Hildreth-Verfahren.

Numerik

Es geht in diesem Kapitel darum, die in Kapitel 3 dargestellten Verfahren zu benutzen, um die in Kapitel 4 vorgestellten Marktgleichgewichtsprobleme zu lösen. Das bedeutet, wir verwenden das erweiterte Hildreth-Verfahren, um ein Quadratisches Netzwerkproblem zu lösen und den PSOR-Algorithmus, um ein nichtlineares Netzwerkproblem zu lösen. Darüberhinaus wird der Baum-Wald-Algorithmus erläutert, allerdings ohne Programmlisting.

6.1 Lösung eines quadratischen Netzwerkes mit dem erweiterten Hildreth-Verfahren

6.1.1 Der numerische Ansatz

Gegeben ist ein Datensatz von Preiselastiziäten b_{i}, Preisen a_{i} und Transportkosten t_{ij}. Bestimmt werden sollen die Gleichgewichtspreise π_{i}^{∗}, i = 1, ..., n der Regionen 1 bis n. Bekannt ist der Zusammenhang:

π_{i}^{∗} = a_{i} − b_{i}y_{i}^{∗}.

y_{i}^{∗}, i = 1, ..., n sind die Nettoimportmengen im Marktgleichgewicht. Um ein LKP (q, M) zu formulieren, benötigt man die Bruttoimportmengen z_{ij}^{∗}, da für die y_{i}^{∗} weder eine Komplementaritätsbedingung, noch y_{j} ≥ 0 angenommen werden kann. Die Netzwerkgleichgewichtsbedingungen sind:

z_{ij}^{∗} ≥ 0
π_{i}^{∗} + t_{ij} − π_{j}^{∗} ≥ 0
z_{ij}^{∗}(π_{i}^{∗} + t_{ij} − π_{j}^{∗}) = 0.

Man definiert über die vorgegebenen Daten:

B := ( \begin{matrix}b_{1} & & & & 0 \\ & . & . & . \\ 0 & & & & b_{n} \end{matrix} ) ∈ IR^{n×n}, a := (a_{1}, ..., a_{n})^{t} ∈ IR^{n}.

KAPITEL 6. NUMERIK

A ist die aus Abschnitt 4.1.2 bekannte Beschränkungsmatrix:

A = ( \begin{matrix}1 & −1 \\ ⋮ & & ... \\ ⋮ & & & ... \\ 1 & & & & −1 \end{matrix} \} n − 1 ; \begin{matrix}−1 & 1 \\ & ⋮ & −1 \\ & ⋮ & & ... \\ & 1 & & & −1 \end{matrix} \} n − 1 ; ⋮  ⋮ ⋮  ⋮ ; \begin{matrix}−1 & & & & 1 \\ & ... & & & ⋮ \\ & & ... & & ⋮ \\ & & & −1 & 1 \end{matrix} \} n − 1 ) ∈ IR^{n(n−1)×n}.

Mit q := t + Aa und M = ABA^{t} läuft die Bestimmung der z_{ij}^{∗} hinaus auf die Lösung des LKP (q, M) (vergleiche Abschnitt 4.2.1). Um die π_{i}^{∗} zu berechnen, sind die Nettoimportmengen y_{i}^{∗} = −(A^{t}z^{∗})_{i}, i = 1, ..., n zu bestimmen. Dazu benutzen wir das zum LKP (q, M) primäre Programm (PQP*) und formen es äquivalent so um, daß das erweiterte Hildreth-Verfahren zur numerischen Lösung benutzt werden kann.

Wir setzen b := B^{−1}a und C = −B^{−1} und betrachten das Optimierungsproblem:

maximiere g(x) = b^{t}x + \frac{1}{2}x^{t}Cx
unter den Nebenbedingungen Ax ≥ −t.

(PQP*)

C = −B^{−1} ist negativ definit. Mit diesen Bezeichnungen gilt: M = −AC^{−1}A^{t} = ABA^{t} und q = t − AC^{−1}b = t + Aa. Das LKP (q, M) ist folglich das zum (PQP*) duale Programm. M ist positiv semidefinit.

Ist x^{∗} eine Lösung des (PQP*), dann gibt es nach Satz 3.7 eine Lösung z^{∗} des LKP (q, M) , so daß:

x^{∗} = −C^{−1}(A^{t}z^{∗} + b) = B(A^{t}z^{∗} + B^{−1}a) = BA^{t}z^{∗} + a ⇔ A^{t}z^{∗} = B^{−1}(x^{∗} − a).

Von −C = B^{−1} bilde man die Cholesky-Zerlegung B^{−1} = D^{t}D. Da B Diagonalmatrix ist mit positiven Diagonalelementen ⇒

B^{−1} = ( \begin{matrix}\frac{1}{b_{1}} & & & & 0 \\ & . & . & . \\ 0 & & & & \frac{1}{b_{n}} \end{matrix} ) ⇒ D = ( \begin{matrix}\frac{1}{\sqrt{b_{1}}} & & & & 0 \\ & . & . & . \\ 0 & & & & \frac{1}{\sqrt{b_{n}}} \end{matrix} ).

Man setzt weiter (vergleiche dazu Abschnitt 3.2.3):

x := D^{−1}ξ − C^{−1}b = D^{−1}ξ + Bb = D^{−1}ξ + a,

Q := −AD^{−1} und c := −d − AC^{−1}b = t + ABb = t + Aa.

Löst ξ^{∗} die Aufgabe

Minimiere \frac{1}{2}∥ξ∥^{2}
unter den Nebenbedingungen Qξ ≤ c,

(EHP)

dann löst x^{∗} = D^{−1}ξ^{∗} + a das (PQP*). Folglich gibt es eine Lösung z^{∗} des LKP (q, M) mit

A^{t}z^{∗} = B^{−1}(D^{−1}ξ^{∗} + a − a) = ( \frac{\sqrt{b_{1}}}{b_{1}} ξ_{1}^{∗} ; ⋮ ; \frac{\sqrt{b_{n}}}{b_{n}} ξ_{n}^{∗} ) ⇔ (A^{t}z^{∗})_{k} = \frac{\sqrt{b_{k}}}{b_{k}} ξ_{k}^{∗} ∀ k = 1, ..., n.

Ziel des Programms ist es, die Gleichgewichtspreise π_{k}^{∗} zu berechnen. Diese ergeben sich wegen y_{k}^{∗} = −(A^{t}z^{∗})_{k} als π_{k}^{∗} = a_{k} − b_{k}y_{k}^{∗} = a_{k} + \sqrt{b_{k}}ξ_{k}^{∗}.

Der numerische Ansatz lautet folglich:

• Löse das ( EHP ) mit dem erweiterten Hildreth-Verfahren, Bezeichne die Lösung als ξ^{∗} ∈ IR^{n}.

• ∀ k = 1, ..., n erhält man die Gleichgewichtspreise als π_{k}^{∗} = a_{k} − b_{k}y_{k}^{∗} = a_{k} + \sqrt{b_{k}}ξ_{k}^{∗}.

6.1.2 Numerische Ergebnisse und Bemerkungen

Damit das erweiterte Hildreth-Verfahren auf das (EHP) angewendet werden kann, müssen zwei Voraussetzungen erfüllt sein. Erstens darf keine Zeile der Matrix Q identisch Null sein; diese Bedingung ist aufgrund der Definition Q = −AD^{−1} erfüllt. Zweitens muß es einen zulässigen Vektor ξ ∈ IR^{n} geben mit Qξ ≤ c. Aufgrund der Problemstellung gibt es stets den zulässigen Vektor ξ_{0} := −Da, denn t ≥ 0 (Transportkosten sind nie negativ) ⇒

t ≥ 0 ⇔ Aa + t ≥ Aa ⇔ Aa + t ≥ −AD^{−1}(−Da) ⇔ Qξ_{0} ≤ c.

Die Datenstruktur im quadratischen Netzwerk erlaubt also immer die Anwendung des erweiterten Hildreth-Verfahrens. Anhand eines übersichtlichen Beispiels prüfen wir die durchgeführten Transformationen:

Für \underline{n=3} ergibt sich mit den Daten

B = ( \begin{matrix}0.04 & 0 & 0 \\ 0 & 0.16 & 0 \\ 0 & 0 & 0.25 \end{matrix} ) , a = (2, 9, 1)^{t}, t = (\frac{1}{2}, \frac{1}{2}, 1, 1, 1, 1)^{t}

der Lösungsvektor ξ^{∗} = (5.4788, −13.5107, 3.1914)^{t}.

Die Gleichgewichtspreise π_{1}^{∗}, π_{2}^{∗}, π_{3}^{∗} lauten also: π^{∗} = ( \begin{matrix}3.0957 \\ 3.5957 \\ 2.5957 \end{matrix} ).

Diese Werte wurden mit dem Programm EHV1.C berechnet (siehe Listing im Anhang). Die niedrige Dimension dieser Aufgabe erlaubt es, die in 6.1.1 durchgeführten Transformationen nachzurechnen und auf diese Weise zu interpretieren. Nach den Vorbemerkungen ist

x^{∗} = D^{−1}ξ^{∗} + a = ( \begin{matrix}0.2 & 0 & 0 \\ 0 & 0.4 & 0 \\ 0 & 0 & 0.5 \end{matrix} ) ξ^{∗} + a = ( \begin{matrix}3.096 \\ 3.596 \\ 2.596 \end{matrix} )

eine Lösung des (PQP*). Ax^{∗} = (−\frac{1}{2}, \frac{1}{2}, \frac{1}{2}, 1, −\frac{1}{2}, −1)^{t} bestätigt, daß Ax^{∗} ≥ −t gilt. Es gibt eine Lösung z^{∗} des LKP (q, M) mit A^{t}z^{∗} = B^{−1}(x^{∗} − a). Für dieses z^{∗} gilt darüberhinaus: Ax^{∗} + t = Mz^{∗} + q. Ax^{∗} + t = (0 , 1, \frac{3}{2}, 2, \frac{1}{2}, 0)^{t} ⇒ z_{i}^{∗} = 0 für i = 2, 3, 4, 5 ⇒ A^{t}z^{∗} = (z_{1}^{∗}, −z_{1}^{∗} − z_{6}^{∗}, z_{6}^{∗})^{t}.

Aus B^{−1}(x^{∗} − a) = ( \begin{matrix}27.4 \\ −33.78 \\ 6.38 \end{matrix} ) folgert man schließlich, daß

z^{∗} = (27.4, 0, 0, 0, 0, 6.38)^{t} eine Lösung des LKP (q, M) ist. Nach diesem übersichtlichen und schnell nachprüfbaren Beispiel geben wir nun ein umfangreicheres der Dimension \underline{n = 15}. Der Datensatz (a_{i}, b_{i}, t_{ij}, i, j = 1, ..., n; i ≠ j) sei gegeben wie folgt:

b_{i} = \frac{1}{i^{2}}
a_{i} = \frac{i mod 10}{i mod 10 + 1}
t_{ij} = \frac{|i − j|}{5}.

Das erweiterte Hildreth-Verfahren liefert den Lösungsvektor ξ_{9}^{∗} = 0.7426, ξ_{10}^{∗} = −8.1749, ξ_{11}^{∗} = 5.3076, ξ_{12}^{∗} = 1.3901 und ξ_{i}^{∗} = 0 für alle restlichen Komponenten. Die Gleichgewichtspreise lauten demzufolge: π^{∗} = (0.5, 0.6667, 0.75, 0.8, 0.83, 0.857, 0.875, 0.889, 0.982, 1.182, 0.982, 0.782, 0.75, 0.8, 0.83)^{t}. Wie man leicht nachrechnet ist für beide Beispiele die Bedingung (4.2) erfüllt. Der Relaxationsparameter ω wurde bei der Berechnung konstant ω = 1.1 gewählt. Wie aus nachfolgender Tabelle hervorgeht ist dies der effizienteste Relaxationsparameter:

ω  0.1  0.2  0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ; # Iterationen 32700 16800 10920 7980 6090 5040 4200 3570 2940 2310

ω 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 ; # Iterationen 1890 2100 2310 2730 3570 4830 6300 9660 18270

Das Listing dieses Beispiels ist im Anhang aufgeführt als EHV2.C.

Das erweiterte Hildreth-Verfahren (siehe Abschnitt 3.2.3) arbeitet nicht mit einem konstanten ω, sondern mit einer Folge von Relaxationsparametern {ω_{k}}_{k}. Interessant wäre die Überlegung, wie sich der Algorithmus effizienter gestalten läßt durch eine “geschickte” Wahl der Folge {ω_{k}}_{k}. Der Ursprungsartikel [Lent ] schweigt sich darüber allerdings aus.

6.2 Der PSOR-Algorithmus als Lösungsverfahren eines nichtlinearen Netzwerkproblems

6.2.1 Der numerische Ansatz

Vorgegeben ist der Transportkostenvektor t ∈ IR^{n(n−1)}, sowie die n nicht-linearen Nachfragefunktionen x_{i}(π_{i}) : IR ⟶ IR. Eine Lösung des nicht-linearen Netzwerkproblems (4.6)-(4.8) bestimmt man nach den Bemerkungen in Abschnitt 4.2.2 mittels einer Folge {y^{(k)}}_{k} von LKP-Lösungen:

y^{(0)} ≥ 0

y^{(k+1)} löst das LKP (q_{k}, M_{k}) mit

q_{k} = ( \begin{matrix}−x(π^{(k)}) + ∇x(π^{(k)})π^{(k)} \\ t \end{matrix} ) und M_{k} = ( \begin{matrix}−∇x(π^{(k)}) & −A^{t} \\ A & 0 \end{matrix} ) ∈ IR^{n^{2}×n^{2}}.

Dabei ist y^{(k)} = (π^{(k)}, z^{(k)})^{t}. A ∈ IR^{n(n−1)×n} ist die Beschränkungsmatrix.

Unser Ziel ist es, Lösungen des LKP (q_{k}, M_{k}) für k = 0, 1, 2, ... mit dem PSOR-Algorithmus zu bestimmen und dies möglichst so, daß die Folge der Lösungsvektoren {y^{(k)}}_{k} konvergiert.

Eine passende Nachfragefunktion x(π) muß gemäß Abschnitt 4.1.2

i) komponentenweise konvex, stetig differenzierbar und monoton fallend sein,
ii) x_{i}(0) ≥ 0 für alle i erfüllen und
iii) es muß ein π̂ ≥ 0 existieren mit ∑_{i=1}^{n} x_{i}(π̂) < 0.

x(π) ist so zu definieren, daß −\frac{∂x(π_{i})}{∂π_{i}} > 0 ∀ π_{i} ≥ 0 gilt,

und damit −∇x(π) positiv definit ist für alle π ≥ 0. Diese Eigenschaft ist schon durch die Bedingung i) gegeben, wenn man für x_{i}(π_{i}) strenge Monotonie fordert. Das Verfahren konvergiert nach dem Beweis von Luo und Tseng, weil die M_{k} für alle k positiv semidefinit sind.

Schwierigkeiten ergeben sich noch aus den Diagonalelementen von M_{k}, die 0 sind. Um dem zu begegnen, verwenden wir anstelle von M_{k} ein M̃_{k} ∈ IR^{n^{2}×n^{2}} mit

M̃_{k} = ( \begin{matrix}−∇x(π^{(k)}) & −A^{t} \\ A & εI \end{matrix} )

I ist die (n(n − 1), n(n − 1))-Einheitsmatrix, und ε > 0 wird möglichst klein gewählt ( so daß der Rechner M̃_{k} noch als positiv definit verarbeitet).

6.2.2 Numerische Ergebnisse und Bemerkungen

Betrachte die Funktion

x_{i}(π_{i}) = \frac{1}{π_{i} + 1} − i.  (6.1)

Es ist −\frac{∂x_{i}(π_{i})}{∂π_{i}} = \frac{1}{(π_{i} + 1)^{2}} , daher ist −∇x(π) positiv definit, und die Bedingungen i)-iii) gelten.

Der Vektor q berechnet sich nach

q_{i} = −x_{i}(π_{i}) + \frac{∂x_{i}(π_{i})}{∂π_{i}}π_{i} = 1 − \frac{2π_{i} + 1}{(π_{i} + 1)^{2}} , für i = 1, ..., n.

Die übrigen n(n − 1) Komponenten von q sind der Vektor t.

Programmiert man das Verfahren mit der Nachfragefunktion (6.1), so hat man schon nach dem ersten Iterationsschritt den Lösungsvektor 0. Dies folgt aus q_{k} ≥ 0 ∀ k = 0, 1, 2, 3, ... (siehe Listing PSOR1.C). Ein solches Ergebnis ist ebenso richtig wie uninteressant.

Ein Lösungsvektor ungleich 0 ergibt sich dann, wenn neben den Eigenschaften i)-iii) gilt:

q_{i} = −x_{i}(π_{i}) + \frac{∂x_{i}(π_{i})}{∂π_{i}}π_{i} < 0

für alle π_{i} ≥ 0 und für mindestens ein i.

Eine in diesem Sinne geschickter als (6.1) gewählte Nachfragefunktion ist die folgende:

x_{i}(π_{i}) = exp (−π_{i}) − \frac{1}{2i}.  (6.2)

i) und ii) sind erfüllt für die Funktion (6.2).

Wählt man π̂_{i} ≥ 0 so groß, daß exp(π̂_{i}) < \frac{1}{2}, so gilt auch iii).

Für \underline{n = 3} und den Transportkostenvektor t = (1, 2, 3, 4, 5, 6)^{t} berechnet man den Lösungsvektor

y^{∗} = (0.6931, 1.3863, 1.7918, 0, 0, 0, 0, 0, 0)^{t}.

Die Gleichgewichtspreise π_{1}^{∗}, π_{2}^{∗}, π_{3}^{∗} lauten also π^{∗} = ( \begin{matrix}0.6931 \\ 1.3863 \\ 1.7918 \end{matrix} ).

Vergleiche dazu das Listing PSOR2.C im Anhang.

6.3 Programmierhinweise zum Baum-Wald-Algorithmus

In Kapitel 4 wurde der Baum-Wald-Algorithmus angegeben. Zu den einzelnen Programmierschritten notieren wir hier einige Erläuterungen.

Bemerkungen zum Algorithmus: \underline{Zu Schritt 1:} Bringt die vorliegende Lösung das Netzwerk ins Gleichgewicht, dann finden sich keine zwei Knoten mehr mit π_{j} + t_{jk} < π_{k}. \underline{Zu Schritt 2:} Bei dem hier bestimmten Indize e handelt es sich tatsächlich um den Knoten einer exportierenden Region und bei i um den einer importierenden Region. Denn: Annahme: j ist kein exportierender Knoten (y_{j} > 0) ⇒ ∃ e ∈ E : z_{ej} > 0. Aus (4.7) und (4.9) folgt dann: \overline{t}_{ek} = π_{e} +t_{ek} −π_{k} ≤ π_{e} +t_{ej} +t_{jk} −π_{k} = π_{e} −π_{e} +π_{j} +t_{jk} −π_{k} = \overline{t}_{jk}. Somit ergibt sich der Widerspruch \overline{t}_{ek} ≤ \overline{t}_{jk} ∀ k ∈ I. Daß k ein importierender Knoten ist, zeigt man auf analoge Weise. Bei der Minimalwertbestimmung in Schritt 2 kann man sich also beschränken auf Knoten exportierender Regionen in E und importierender in I. Zu Schritt 2 (a) Es sei j ∈ E. Wir betrachten Nettoimporte im Sinne von Abschnitt 4.1.2. y_{j}(0) ist der Nettoimport von Region j, bevor der Export z_{ei} stattfindet. y_{j}(z_{ei}) ist der Nettoimport, nachdem z_{ei} erfolgt ist. Setze △y_{j} := y_{j}(z_{ei}) − y_{j}(0) ∀ j ∈ E ∪ I. Damit e exportieren und i importieren kann, ändern sich die Netzimporte in E bzw. in I in Höhe von z_{ei}:

− ∑_{j∈E} △y_{j} = z_{ei} = ∑_{k∈I} △y_{k}.  (6.3)

Definiere: {δ_{e} = π_{e}(z_{ei}) − π_{e}(0), \atop δ_{i} = π_{i}(z_{ei}) − π_{i}(0).}

π_{e}(0) ist Preis in r_{e} vor dem Export z_{ei}. π_{e}(z_{ei}) ist Preis in Region r_{e} nach dem Export z_{ei}. Da j ∈ E ⇒ z_{ej} > 0. ⟹^{(4.9)} π_{e}(0) − π_{j}(0) + t_{ej} = 0 und π_{e}(z_{ei}) − π_{j}(z_{ei}) + t_{ej} = 0 (Der Ast fällt nicht weg). ⇒ π_{e}(0) − π_{j}(0) = π_{e}(z_{ei}) − π_{j}(z_{ei}). Es ist also δ_{e} = π_{j}(z_{ei}) − π_{j}(0) ∀ j ∈ E. Analog zeigt man: δ_{i} = π_{k}(z_{ei}) − π_{k}(0) ∀ k ∈ I. Ist δ_{e} positiv, dann ist δ_{i} negativ. Mit (4.9) und (4.10) folgern wir δ_{e} = π_{j}(z_{ei}) − π_{j}(0) = −b_{j}(y_{j}(z_{ei}) − y_{j}(0)) = −b_{j}△y_{j} ∀ j ∈ E und δ_{i} = −b_{k}(y_{k}(z_{ei} − y_{k}(0)) = −b_{k}△y_{k} ∀ k ∈ I. Es folgt:

δ_{e} ∑_{j∈E} β_{j} = z_{ei} = −δ_{i} ∑_{k∈I} β_{k},  (6.4)

wobei β_{j} = \frac{1}{b_{j}}.

Zu Schritt 2 (b) π_{e}(z_{ei}) + t_{ei} = π_{i}(z_{ei}) gilt dann, wenn

δ_{i} − δ_{e} = π_{e}(0) + t_{ei} − π_{i}(0) =: \overline{t}_{ei}  (6.5)

gilt. Über Straight-Forward-Rechnungen und der Trivialität ∑_{j∈E} β_{j} + ∑_{k∈I} β_{k} = ∑_{k∈E∪I} β_{k} erhält man aus (6.4) und (6.5):

δ_{e} = −\overline{t}_{ei} \frac{∑_{k∈I} β_{k}}{∑_{k∈E∪I} β_{k}}  (6.6)

δ_{i} = \overline{t}_{ei} \frac{∑_{k∈E} β_{k}}{∑_{k∈E∪I} β_{k}} .  (6.7)

Und daraus dann:

z_{ei} = \overline{t}_{ei} ∑_{k∈E} β_{k} \frac{∑_{j∈I} β_{j}}{∑_{k∈E∪I} β_{k}} .  (6.8)

Zu Schritt 2 (c) Der mittels (6.8) bestimmte Wert von z_{ei} ist nicht richtig, wenn er größer ist als einer der beiden kritischen Werte von E bzw. I. Dieser kritische Wert gibt an, bei welcher Höhe von z_{ei} es einen Ast (j, k) in E bzw. I gibt mit z_{jk} = 0. Dazu betrachte einen Ast (j, k) ∈ I, der auf dem Pfad zwischen i und j die (e, i) entgegengesetzte Richtung hat ^{1}. T_{j} bezeichne den Teilbaum von I, der alle Pfade zusammenfaßt, die von Punkt j ausgehen, läßt man vorher (j, k) weg. Steigt z_{ei} an, dann reduziert sich z_{jk}. Es gilt die Erhaltung der Güterbewegung:

z_{jk} = − ∑_{n∈T_{j}} y_{n}.  (6.9)

Als Folgerung aus (4.9), (6.4) und (6.9) erhält man △z_{jk} =

= z_{jk}(z_{ei}) − z_{jk}(0) =^{(6.9)} − ∑_{n∈T_{j}}(y_{n}(z_{ei}) − y_{n}(0)) =^{(4.9)} − ∑_{n∈T_{j}} β_{n}(−π_{n}(z_{ei}) + π_{n}(0)) =

= ∑_{n∈T_{j}} β_{n}(π_{n}(z_{ei}) − π_{n}(0)) = δ_{i} ∑_{n∈T_{j}} β_{n} =^{(6.4)} −z_{ei} \frac{∑_{n∈T_{j}} β_{n}}{∑_{n∈I} β_{n}}.

1Ein Pfad zwischen 2 Punkten i und j ist eine Verbindung gleichgerichteter Wege mit Anfangspunkt i und Endpunkt j.

Demzufolge lautet der kritische Wert von z_{ei}, bei dem z_{jk} 0 wird:

z_{jk}^{krit} = z_{jk} \frac{∑_{n∈I} β_{n}}{∑_{n∈T_{j}} β_{n}} .  (6.10)

Analog zeigt man für einen (j, k) ∈ E, der (e, i) entgegengesetzt ist:

z_{jk}^{krit} = z_{jk} \frac{∑_{n∈E} β_{n}}{∑_{n∈T_{k}} β_{n}} .  (6.11)

Die Programmierung des Baum-Wald-Algorithmus verläuft also so, daß man die beiden kritischen Werte und z_{ei} nach (6.8) berechnet. In jeder neu gebildeten Handelsgruppe werden nur die Gleichgewichtspreise π_{j} auf den aktuellen Stand gebracht. Die Nettoimportmenge y_{j} läßt sich ggf. über (4.9) herleiten.

Schlußbemerkung

Die vorliegende Arbeit ist entstanden aus einem Vortrag über Lineare Komplementaritätsprobleme; dieser war Bestandteil eines Seminars über Partielle Differentialgleichungen.

Um LKP s geht es in der Theorie partieller Differentialgleichungen, wenn ein Anfangs-Randwertproblem mittels einer Variationsungleichung diskretisiert wird.

Da die Methode der finiten Elemente zur Lösung einer Differentialgleichung zunehmend an Bedeutung gewinnt, dürften auch iterative LKP -Verfahren in der Numerischen Mathematik immer wichtiger werden.

Geht man nach der Interpretation von [Neumann ], so läßt sich Netzwerkoptimierung der Spieltheorie zuordnen. Denn von mehreren Spielern ( = Regionen) versucht jeder, seinen Gewinn zu maximieren unter der Prämisse, daß dies jedem bekannt ist. Es kommt hinzu, daß jeder Spieler nicht nur Handelsware verkauft, um Gewinn zu erzielen, sondern u.U. auch Handelsware kauft, um Bedarf zu decken. Insofern ist das Netzwerkproblem eine Idee näher an der Realität als vergleichbare Modelle der Spieltheorie.

Es gibt noch eine ganze Reihe von Möglichkeiten, das Netzwerkproblem zu diskutieren und zu interpretieren. So könnte man beispielsweise weitere Restriktionen an die effektive Nachfrage x_{i}(π_{i}) stellen, um das Modell realitätsnäher zu machen. Möglich ist es auch, Ergebnisänderungen durch unterschiedliche Datenvorgaben näher zu untersuchen und anhand empirischer Werte zu analysieren.

Dies hätte aber den Rahmen dieser Arbeit gesprengt.

Anhang

Symbolverzeichnis

Analytische Bezeichnungen ∥ · ∥ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . euklidische Norm des IR^{n}, d.h. ∥x∥ = \sqrt{x^{t}x} ∀ x ∈ IR^{n}.

∥ · ∥_{∞} . . . . . . . . . . . . . . . . . . . . . . . . . . . . Maximumnorm des IR^{n}, d.h. ∥x∥_{∞} = max_{i∈{1,...,n}}{x_{i}} ∀ x ∈ IR^{n}.

| · | . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . L_{1}-Norm: |y| = ∑_{i=1}^{n} |y_{i}| ∀ y ∈ IR^{n}.

IR_{+}^{n} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x ∈ IR_{+}^{n} ⇔ x ∈ IR^{n} ∧ x ≥ 0.

[ · ]^{+} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Für x ∈ IR^{n} ist [ x ]_{i}^{+} = \{ \begin{matrix}x_{i} & , wenn x_{i} ≥ 0 \\ 0 & , sonst \end{matrix} , i = 1, ..., n.

∇f(z) . . . . . . . . . . . . . . . . . . . Gradient von f an der Stelle z, d.h. ∇f(z) = (\frac{∂}{∂z_{1}} f(z), ..., \frac{∂}{∂z_{n}} f(z))^{t}.

∇_{z}f(z, y) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gradient von f bzgl. z an der Stelle (z, y).

Untermatrizen und Untervektoren I, J ⊆ {1, ..., n} seien Indexmengen. IR^{n×m} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Menge der (n, m)-Matrizen mit reellen Einträgen. C^{n×m} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Menge der (n, m)-Matrizen mit komplexen Einträgen. |J| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Anzahl der Elemente der Menge J. J^{C} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Das Komplement von J bzgl. {1, ..., n}. I_{L} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . zur Definition siehe S.63. B_{JJ} ∈ IR^{|J|×|J|} . . . . . . . . . Untermatrix von B, die die mit J indizierten Zeilen und Spalten enthält. B_{IJ} ∈ IR^{|I|×|J|} . . . . . . Untermatrix von B, mit den mit I indizierten Zeilen/ J indizierten Spalten. B_{J·} ∈ IR^{|J|×n} . . . . . . . . . . . . . . . . . . . . . . . Untermatrix von B, die die mit J indizierten Zeilen enthält. B_{·J} ∈ IR^{n×|J|} . . . . . . . . . . . . . . . . . . . . . . Untermatrix von B, die die mit J indizierten Spalten enthält.

Bezeichnungen im Marktgleichgewichtsmodell π_{i} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Preis für ein Stück Handelsware in Region i. π_{i}^{∗} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gleichgewichtspreis in Region i. x_{i}(π_{i}) . . . . . . . . . . . . . . . . . . . . . . . . Anzahl der effektiv nachgefragten Stücke Handelsware in Region i. s_{i}(π_{i}) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Anzahl der angebotenen Stücke Handelware in Region i. d_{i}(π_{i}) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Anzahl der nachgefragten Stücke Handelsware in Region i. t_{ij} . . . . . . . . . . . Transportkosten, um ein Stück Handelware von Region i nach Region j zu bringen. z_{ij} . . . . . . . . . Anzahl der Stücke Handelware, die von Region i nach Region j transportiert werden. z_{ij}^{∗} . . . . . . . . . . . . . . Von Region i nach Region j transportierte Handelsware im Marktgleichgewicht. a_{i} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gleichgewichtspreise ohne Im/Exporte. b_{i} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Elastizität der Nachfrage in Region i.

Abkürzungen im fortlaufenden Text BS-Methode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Basis-Splitting-Methode. bzgl. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . bezüglich. bzw. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . beziehungsweise. d.h. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . daher. FA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Funktionalanalysis. ggf. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . gegebenenfalls. i.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . im folgenden. i.B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . im Besonderen. KKT-Bedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Karush-Kuhn-Tucker-Bedingungen. m.E. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . mit Einschränkung. N.B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Nebenbedingungen. Œ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ohne Einschränkung. u.U. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . unter Umständen. w.z.b.w. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . was zu beweisen war. z.B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . zum Beispiel. z.T. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . zum Teil.

Literaturquellen

[Alt ] Hans-Wilhelm Alt, Lineare Funktionalanalysis, 2.Auflage, Springer

[Cottle ] R.W. Cottle, J.-S. Pang, R.E. Stone, The linear complementarity problem, Academic press 1992

[Cryer ] C.W. Cryer, The solution of a quadratic programming problem using SOR, SIAM J. Contr. 9, pp.385-392

[Daniel ] J.W. Daniel, On pertubations in systems of linear inequalities, SIAM J.Num.An. 10 (1973), pp.299-307

[Everett ] H.Everett, Generalized Lagrangemultiplier method for solving problems of allocations of resources, OR 11 (1963), pp.399-417

[Farkas ] J.Farkas, Theorie der einfachen Ungleichungen, J. reine u. angew. Mathematik 124, S.1-24

[Fuhrmann ] W. Fuhrmann, Makroökonomik - Zur Theorie interdependenter Märkte, Oldenbourg 3.Auflage (1991)

[Frank ] M.Frank, P.Wolfe, An algorithm for quadratic programming, Nav.Res.Log.Quat. 3 (1956) pp.95-110

[Gale ] D. Gale, The theory of linear economic models , McGraw-Hill, 1960

[Glassey ] C.R. Glassey, A quadratic network optimization model for equilibrium single commodity trade flows, Math. Progr. 14 (1978), pp.98-117

[Hildreth ] C. Hildreth, A quadratic programming procedure, Nav. Res. Log. Quart. 4 (1957), pp.79-85

[Jones ] P.C. Jones, R.Saigal, M.Schneider, Computing nonlinear network equilibria, Math.Progr. 31 (1985), pp.57-66

[Keller ] H.B. Keller, On the solution of singular and semidefinite linear systems by iteration, SIAM J. Num. Analysis Ser.B 2, pp.281-291 [Künzi ] H.P. Künzi, W. Krelle, R. von Randow,

Nichtlineare Programmierung, 2.Auflage, Springer 1979

[Lent ] A.Lent, Y.Censor, Extensions of Hildreth‘s row-action method for quadratic programming, SIAM J. Contr. Opt. 18 pp. 444-453

[Luo ] Z.-Q. Luo, P. Tseng, On the convergence of a matrix splitting algorithm for the symmetric monotone linear complementarity problem, SIAM J. Contr. Opt. 29 pp.1037-1060

[Manga ] O.L. Mangasarian, Nonlinear programming, McGraw-Hill 1969

[Matth 1 ] L.Matthiesen, An algorithm based on a sequence of linear complementarity problems applied to a walrasian model: an example, Math. Progr. 37 (1987) pp. 1-18

[Matth 2 ] L.Matthiesen, Computation of economic equilibria by a sequence of linear complementarity problems, Math. Progr. St. 23 pp. 144-162

[Neumann ] K.Neumann, Operations Research Verfahren Band I, Carl Hanser 1975

[Pang ] J.S. Pang, More results on the convergence of iterative methods for the symmetric linear complementarity problems, J. Opt. Th. Appl. 49, pp.107-134

[Schwarz ] H.R. Schwarz, Numerische Mathematik, Teubner 2.Aufl. 1993

[Varga ] R.S. Varga, Matrix Iterative Analysis, Prentice-Hall 1962

Programmlistings

PROJAC.C

/*  Programm zum projezierten Jacobi-Verfahren
               2  2  1       -1
LKP(q,M), mit M= 2  3  1    q=  2
               1  1  1       -2                       */
#include<stdio.h>
#include<math.h>
#define eps 0.0001
main()
{  float z[3]={0.0, 0.0 , 0.0}, zneu[3], test;
int i, j=0;
do
{  zneu[0]=0.5-z[1]-0.5*z[2];
zneu[1]=-(2./3)-(2./3)*z[0]-(1./3)*z[2];
zneu[2]=2-z[0]-z[1];
test=0;
for (i=0;i<3;i++)
 { if (zneu[i]<0)
    { zneu[i]=0.;
      printf("\n  Der %d-te Vektor ist kleiner als Null",i+1);}}
for (i=0;i<3;i++)
 test=test+(z[i]-zneu[i])*(z[i]-zneu[i]);
test=sqrt(test);
for (i=0;i<3;i++)
 z[i]=zneu[i];
j++;
if (j>100)
 { printf("\n Es klappt wohl nicht so mit der Iteration!");
   break;}
}
while (test>eps);
printf(" Eine Loesung gefunden. Sie lautet:\n");

   for (i=0;i<3;i++)
    printf(" %f\n",z[i]);
   printf("\n Anzahl der Iterationsschritte: %d",j);
}

EHV1.C

/* Programmlisting EHV1.C. Verfahren zur Loesung einer Netzwerkstruktur mit dem
   erweiterten Hildreth-Verfahren.
   Das Problem hat die Dimension 3 und dient dazu, die durchgefuehrten
   Transformationen zu ueberpruefen.                                    */
#include<stdio.h>
#include<math.h>
#define eps 0.001
#define n 3
void ehvvek();
void ehvmat();
void preise();
double omega=1.1;       /*  konstante Folge von Relaxationsparametern   */
double a[n]={2., 9., 1.};                /*  Der vorgegebene Datensatz   */
double t[n*(n-1)]={.5 ,.5 ,1. ,1. ,1. ,1.};
double b[n]={.04, .16, .25};
double Q[n*(n-1)][n], c[n*(n-1)];
main()
{  double x[n], altz[n*(n-1)],
          z[n*(n-1)]={.0, .0, .0, .0, .0, .0},  /* Initialisierung des
                                                   Startvektors z(0)  */
          r, norm;
   int    count=0;
   register i, k;
/* Zunaechst bestimmt man aus dem Datensatz fuer a, b, t die beim erweiterten
   Hildreth-Verfahren verarbeitete Matrix Q und den Vektor c. Diese berechnen
   wir in den Funktionen ehvmat() und ehvvek() (siehe unten)              */

ehvmat(Q);
ehvvek(c);
/* Berechne den Startvektor x(0) aus x(0)=-Qtz(0)  */
for (i=0; i<n; i++)
 { norm=0.0;
   for (k=0; k<n*(n-1); k++)
    norm=norm+Q[k][i]*z[k];
   x[i]=-norm;
 }
/*   Hier beginnt der Hildreth-Algorithmus. Jeder Iterationsschritt besteht
 aus n(n-1) Einzelschritten. Abbruchbedingung ist, dass die Abweichung
 der Norm im x-Vektor kleiner als eps ist.                             */
do
{
for (i=0; i<n*(n-1); i++)
{ norm=0.0;
  r=0.0;
/*  Berechnung von r  */
for (k=0; k<n; k++)
 { r-=Q[i][k]*x[k];
   norm+=Q[i][k]*Q[i][k];}
r=(omega*(r+c[i]))/norm;
if (r>=z[i]) r=z[i];
/*   neuer Iterationsvektor x   */
for (k=0; k<n; k++) x[k]+=r*Q[i][k];
/*   neuer Iterationsvektor z   */
altz[i]=z[i];
z[i]-=r;
count++;
}
/*   Berechnung der Norm der Abweichung ---> Abbruchbedingung   */
norm=0.0;
for (k=0; k<n*(n-1); k++) norm+=(altz[k]-z[k])*(altz[k]-z[k]);
norm=sqrt(norm);

}
while ((norm>eps)&&(count<1000));
if (count>=1000) printf("\n   Das waren mehr als 1000 Iterationen !");
else
 {  printf("\n   Der Loesungsvektor x lautet:\n\n");
    for (i=0; i<n; i++) printf("   %10.4f\n",x[i]);
 }
printf("\n   Es wurden %d Iterationen durchgefuehrt.\n",count);
/*  Schliesslich: Berechnung und Ausgabe der Gleichgewichtspreise in
    der Funktion preise().                                           */
preise(x);
return;
}
/*  Die Funktion ehvmat() bestimmt die Matrix Q im EHV        */
void ehvmat(mat)
double mat[n*(n-1)][n];
{ int k, l, K, L;
for (k=0; k<n*(n-1); k++)
  { K=k/(n-1);
    L=k%(n-1);
    for (l=0; l<n; l++)
      { mat[k][l]=0.0;
        if (l==K)              mat[k][l]=-sqrt(b[K]);
        if ((l<K) && (L==l))   mat[k][l]= sqrt(b[L]);
        if ((l>K) && (L+1==l)) mat[k][l]= sqrt(b[L+1]);
      }
  }
return;
}
/*  Die Funktion beschvek() bestimmt den Vektor c im EHV      */
void ehvvek(vektor)
double vektor[n*(n-1)];
{  int index=0;
register k, l;
for (k=0; k<n; k++)

     { for (l=0; l<n-1; l++)
         { if (l<k) vektor[index]=a[k]-a[l]+t[index];
           else   vektor[index]=a[k]-a[l+1]+t[index];
           index++;
         }
      }
    return;
}
/*   Berechnung der Gleichgewichtspreise         */
void preise(vektor)
double vektor[n];
{  register k;
   printf("\n   Der %d-Vektor der Gleichgewichtspreise lautet:\n\n",n);
   for (k=0; k<n; k++)
     { vektor[k]=sqrt(b[k])*vektor[k]+a[k];
       printf("   %10.4f\n",vektor[k]);}
   return;
}

EHV2.C

Das Listing entspricht fast wortwörtlich dem Listing EHV1.C. Lediglich die Daten der Vektoren t, a und b werden als globale Größen in der Funktion dsatz() bestimmt.

...
main()
{  ...
/* Initialisierung des Startvektors z(0)    */
   for (i=0; i<n*(n-1); i++) z[i]=0.0;
   dsatz(a,b,t);
   ...
}
/*   Berechnung des vorgegebenen Datensatzes b, a, t       */
void dsatz(v,vek,vektor)
double vektor[n*(n-1)], v[n], vek[n];

{  register i, j;
   int index=0;
   double mod=0.0, stuetze;
   for (i=0; i<n; i++)
   /*   Zunaechst die Transportkosten   */
    {  for (j=0; j<n; j++)
        { if (i!=j)
           { stuetze=i-j;
             vektor[index]=stuetze/5;
             if (vektor[index]<0) vektor[index]*=-1;
             index++;
           }
        }
    /*    Jetzt die Preiselastizitaeten      */
       vek[i]=1./((i+1)*(i+1));
    /*   Jetzt die Preise ohne Im-/Exporte   */
      if (((i+1)%10)!=0)
       { v[i]=(mod+1)/(mod+2);
         mod++;}
     else
      { v[i]=2.;
        mod=0.;}
    }
   return;
}

PSOR1.C

#define n 3              /* Dimension des LKP festlegen  */
#define eps 0.001       /* Genauigkeit  */
#include <stdio.h>
#include <math.h>
void mataus();
void lkpmat();
void psor();
void beschmat();

void xpiblock();
void qk();
double t[n*(n-1)];
double omega=1.,     /* Setze den Relaxationsparameter */
norm;
double Q[n*n][n*n], lkpvek[n*n], itvek[n*n], diag[n];
int count=0;
main()
{  double   hilfsvek[n];
register lauf,i;
/* Initialisierung des Startvektors als 0     */
for (lauf=0;lauf<n*n;lauf++) itvek[lauf]=0.0;
/* Initialisierung des Transportvektors   */
for (lauf=0;lauf<n*(n-1);lauf++) t[lauf]=lauf+1;
/* xpiblock() und qk() verarbeiten die ersten n Iterationskomponenten */
for (lauf=0; lauf<n; lauf++) hilfsvek[lauf]=itvek[lauf];
do
{
/* Berechnung des Vektors q, der im LKP(q,M) verarbeitet wird. */
qk(hilfsvek);
xpiblock(hilfsvek);
lkpmat(Q,diag);
psor(itvek,Q,lkpvek);
/* Bestimme die Norm der Abweichung zwischen neuer und alter
   Iteration im Rahmenprogramm                               */
count++;
norm=.0;
for (i=0;i<n;i++)
 norm=norm+(hilfsvek[i]-itvek[i])*(hilfsvek[i]-itvek[i]);
norm=sqrt(norm);
 /*   Zuordnung des neuen Iterationsvektors   */

for (lauf=0; lauf<n; lauf++) hilfsvek[lauf]=itvek[lauf];
}
while ((norm>eps)&&(count<1000));
printf("\n   Der Loesungsvektor lautet:\n\n");
for(lauf=0; lauf<n*n; lauf++) printf(" %10.4f\n", itvek[lauf]);
printf("\n   Es wurden %d Iterationen durchgefuehrt!", count);
return;
}
/* Der PSOR-Algorithmus. Uebergeben wird der Startvektor z, die LKP-Matrix M
und der LKP-Vektor q; geliefert wird der Loesungsvektor zstern            */
void psor(z,M,q)
double z[n*n], q[n*n], M[n*n][n*n];
{ double y[n*n],
      altz[n*n],            /* Kontrollvektor zur Ueberpruefung der  */
                            /* Genauigkeit                           */
      sorvek, r;
register i,j,k;
do
{
/* Ein Iterationsschritt im PSOR-Algorithmus */
/* Er besteht aus n Einzelschritten          */
for (i=0;i<n*n;i++)
 {  r=.0;
    for (j=0;j<n*n;j++)
     r=r+M[i][j]*z[j];
    r=-q[i]-r;
    sorvek=z[i]+omega*(r/M[i][i]);
    altz[i]=z[i];
    if (sorvek<0)
     z[i]=0;
    else
     z[i]=sorvek;
    y[i]=-r+M[i][i]*(z[i]-altz[i]);
 }
/* Bestimme die Norm der Abweichung zwischen neuer und alter */
/* Iteration                                                 */

norm=.0;
for (i=0;i<n*n;i++)
 norm=norm+(altz[i]-z[i])*(altz[i]-z[i]);
norm=sqrt(norm);
/* Abbruchbedingung: Norm der Abweichung kleiner als eps */
}
while (norm>eps);
}
/*  Die Funktion beschmat() bestimmt die Beschraenkungsmatrix in Abhaengigkeit
 von n, die zur Berechnung der LKP-Matrix M benoetigt wird               */
void beschmat(mat)
double mat[n*(n-1)][n];
{ int k, l, K, L;
 for (k=0; k<n*(n-1); k++)
   { K=k/(n-1);
     L=k%(n-1);
     for (l=0; l<n; l++)
       { mat[k][l]=0.0;
         if (l==K)              mat[k][l]= 1;
         if ((l<K) && (L==l))   mat[k][l]= -1;
         if ((l>K) && (L+1==l)) mat[k][l]= -1;
       }
   }
return;
}
/*  In der Funktion lkpmat() wird die Matrix M, die im LKP(q,M) verarbeitet
 werden wird, bestimmt aus der Beschraenkungsmatrix A und dem n-Vektor
 der partiellen Ableitungen der Nachfragefunktionen                     */
void lkpmat(mat,vek)
double mat[n*n][n*n], vek[n];
{ int i, j;
double A[n*(n-1)][n];
beschmat(A);
for (i=0; i<n*n;i++)
{ for (j=0; j<n*n; j++)
   { if (i<n)
      { if (j<n)

            { if (i==j) mat[i][j]=vek[i];
              else mat[i][j]=0.0;
            }
           else
            mat[i][j]=-A[j-n][i];
         }
        else
         { if (j<n) mat[i][j]=A[i-n][j];
           else
            { if ((j-n)==(i-n)) mat[i][j]=9.0;
              else mat[i][j]=0.0;
            }
         }
      }
   }
 return;
}
/*   Die Funktion xpiblock() bestimmt die Diagonalwerte der linken oberen
     Dreiecksmatrix aus der aktuellen Iteration zur Berechnung der neuen
     Iteration.
     Die Nachfragefunktion lautet hier:
     x_i(pi_i)=1/(pi_i+1)-1                                               */
void xpiblock(vek)
double vek[n];
{ int i;
  for (i=0; i<n; i++)
   diag[i]=1./((vek[i]+1)*(vek[i]+1));
  return;
}
/*  In der Funktion qk() berechnen wir den aktuellen LKP-Vektor aus den
    vorgegebenen Werten                                                 */
void qk(vek)
double vek[n];
{ int i;
  /* Die ersten n Komponenten haengen ab von der Nachfragefunktion. */
  for(i=0; i<n; i++)
   lkpvek[i]=i+1-((2*vek[i]+1)/((vek[i]+1)*(vek[i]+1)));
  /* Die restlichen n(n-1) Komponenten sind der Vektor t.  */

  for(i=n;i<n*n;i++) lkpvek[i]=t[i-n];
  return;
}

PSOR2.C

Das Listing entspricht fast wortwörtlich dem Listing PSOR1.C. Da eine andere Nachfragefunktion benutzt wird, sind lediglich die Funktionen qk() und xpiblock() unterschiedlich.

/*   Die Funktion xpiblock() bestimmt die Diagonalwerte der linken oberen
     Dreiecksmatrix aus der aktuellen Iteration zur Berechnung der neuen
     Iteration.
     Die Nachfragefunktion lautet hier:
     x_i(pi_i)=exp(-pi_i)-1/(2i)                                           */
void xpiblock(vek)
double vek[n];
{ int i;
  for (i=0; i<n; i++)
   diag[i]=exp(-vek[i]);
  return;
}
/*  In der Funktion qk() berechnen wir den aktuellen LKP-Vektor aus den
    vorgegebenen Werten                                                 */
void qk(vek)
double vek[n];
{ int i;
  double itrans;
  /* Die ersten n Komponenten hngen ab von der Nachfragefunktion. */
  for(i=0; i<n; i++)
   { itrans=i+1;
     lkpvek[i]=-exp(-vek[i])*(1+vek[i])+(1./(2*itrans));}
  /* Die restlichen n(n-1) Komponenten sind der Vektor t.  */
  for(i=n;i<n*n;i++) lkpvek[i]=t[i-n];
  return;
}

Hiermit versichere ich, diese Diplomarbeit eigenständig und ohne unzulässige Hilfsmittel erstellt zu haben.

Münster, im Februar 1995