Benutzer:Stepri2005/Kurs:Optimierung/Der primale Affine-Scaling-Algorithmus

6.1 Zur Geschichte der linearen Optimierung

Der einzig relevante Algorithmus zur Lösung linearer Optimierungsprobleme war 40 Jahre lang der Mitte der 40er Jahre von Dantzig entwickelte und oben vorgestellte Simplexalgorithmus. Dieser Algorithmus macht sich zunutze, dass der Rand des Polyeders $Z_{P}$ endliche viele Ecken hat und dass ein lineares Optimierungsproblem der Form

${\begin{array}{lll}(P):&{\text{Minimiere}}&c^{T}x\\&{\text{u. d. N.}}&Ax=b,\\&&x\geq 0,\end{array}}$

sofern es überhaupt lösbar ist, auch eine solche Ecke als Lösung besitzt. Ausgehend von einer Startecke durchläuft der Simplexalgorithmus eine gewisse Anzahl von jeweils benachbarten Ecken so, dass der Zielfunktionswert von einer Ecke zur nächsten zumindest nicht größer wird und daher nach endlich vielen Schritten eine optimale Ecke erreicht wird.

Nun gibt es Polyeder, welche eine bezüglich $n$ exponentielle Anzahl von Ecken haben. Ein solches Polyeder ist der durch $2n$ Ungleichungen beschreibbare Einheitskubus

$\{x\in \mathbb {R} ^{n}{\big |}0\leq x_{i}\leq 1\quad (i=1,\ldots ,n)\},$

der offenbar $2^{n}$ Ecken besitzt, die nacheinander ohne Mehrfachkontakt durchlaufen werden können. Im Jahr 1972 zeigten Klee und Minty, dass der Simplexalgorithmus zur Lösung eines gewissen Problems, welches darin besteht, die Funktion $f(x):=-x_{n}$ über einem „etwas verbeulten“ Einheitskubus zu minimieren, bei Wahl einer bestimmten Startecke alle Ecken des zulässigen Gebietes dieses Problems abwandert und folglich $2^{n}-1$ Iterationen benötigt (siehe z. B. [BerTsi97]). Damit war bewiesen, dass der Simplexalgorithmus im ungünstigsten Fall eine in Bezug auf die Zahl der Variablen exponentielle Anzahl von Iterationen zur Lösung eines linearen Optimierungsproblems aufwenden muss.

Es stellte sich nun die Frage, ob ein derartiges „Worst-Case-Verhalten“ ein Charakteristikum des Simplexalgorithmus ist oder ob es eine Eigenschaft ist, die dem Problem der linearen Optimierung selbst innewohnt und damit für jeden anderen Algorithmus für lineare Optimierungsprobleme auch nachgewiesen werden kann. Die Beantwortung einer solchen Frage, d. h. die Erfassung des Aufwands, der benötigt wird, um ein Problem einer gegebenen Problemklasse mit vorgegebener Genauigkeit zu lösen, fällt in das Gebiet der in den 70er Jahren begründeten Komplexitätstheorie.

In der Komplexitätstheorie unterscheidet man die Worst Case Analysis, bei der es darum geht, den Aufwand zu bemessen, der zur Lösung des schwierigsten Problems vorgegebener Größe einer Problemklasse erforderlich ist, von der Average Case Analysis, die sich mit dem durchschnittlichen Aufwand zur Lösung eines Problems fester Größe bezüglich einer Problemklasse beschäftigt. Für die Praxis sind Ergebnisse einer Average Case Analysis sicher die aussagefähigeren. Sie sind aber auch sehr viel schwieriger zu gewinnen, weswegen meistens auch nur eine Worst Case Analysis betrieben wird.

Grob gesprochen sagt man, dass eine Problemklasse polynomiale Komplexität besitzt, wenn jedes Problem der Größe „$n$“ (z. B. der Variablenzahl $n$ bei einem Optimierungsproblem) durch einen Algorithmus in höchstens $p(n)$ „elementaren Rechenoperationen“ im Rahmen einer gewissen Genauigkeit gelöst werden kann, wobei $p$ ein Polynom ist. Verfügt insbesondere ein spezieller Algorithmus über diese Lösungseigenschaft, so sagt man, dass er eine polynomiale Laufzeit hat. Die polynomiale Komplexität einer Problemklasse folgt natürlich, wenn für sie ein Algorithmus mit polynomialer Laufzeit gefunden wurde.

Für die lineare Optimierung war durch das Beispiel von Klee und Minty nachgewiesen worden, dass der Simplexalgorithmus keine polynomiale, sondern nur eine exponentielle Laufzeit besitzt. (Ob dies für jede denkbare Pivotregel gilt, ist noch ungeklärt. Für die wichtigsten Pivotregeln ist dies aber gezeigt worden.) Die viele Jahre offene Frage im Rahmen der Komplexitätstheorie war also, ob das lineare Optimierungsproblem polynomiale oder exponentielle Komplexität hat oder, anders ausgedrückt, ob es Algorithmen zur Lösung linearer Optimierungsaufgaben mit polynomialer Laufzeit gibt. Dass die Beantwortung dieser Frage grundsätzlich von Interesse ist, zeigt das folgende Beispiel.

Beispiel 6.1

Angenommen, man hat die Wahl zwischen zwei Algorithmen, von denen der erste $10^{n}$ (exponentielle Laufzeit) und der zweite $n^{3}$ (polynomiale Laufzeit) Rechenoperationen benötigt, um ein Problem in $n$ Variablen einer bestimmten Problemklasse zu lösen. Hat man einen Computer zur Verfügung, der $10^{9}$ Rechenoperationen pro Sekunde durchführen kann, und ist man bereit, diesen $1\ 000$ Sekunden für die Lösung eines Problems arbeiten zu lassen, so errechnet man

$10^{n}=10^{9}\cdot 1\ 000\Rightarrow n=12,$

$n^{3}=10^{9}\cdot 1\ 000\Rightarrow n=10\ 000.$

Das bedeutet, dass der erste Algorithmus in dieser Zeit Probleme bis $n=12$ und der zweite Probleme bis $n=10\ 000$ lösen kann.

---

Im Jahr 1979 veröffentlichte Khachian einen neuen Algorithmus zur Lösung linearer Optimierungsprobleme, die sog. Ellipsoidmethode, welche eine polynomiale Laufzeit besitzt. Damit war geklärt, dass das Problem der linearen Optimierung polynomiale Komplexität hat. Leider stellte sich sehr schnell heraus, dass die auch im SPIEGEL im Dezember 1979 verbreitete Vermutung „Wahrscheinlich wird sich durch Khachian’s Erkenntnis die dagegen doch recht umständliche Simplexmethode weithin ersetzen lassen“ unzutreffend war und die Ellipsoidmethode dem Simplexalgorithmus in der Praxis fast immer deutlich unterlegen ist. Dies hängt unter Anderem damit zusammen, dass das Beispiel von Klee und Minty nicht das typische, sondern nur das in Extremfällen mögliche Verhalten des Simplexverfahrens aufzeigt.

Inzwischen vorliegende Untersuchungen über durchschnittliche Laufzeiten des Simplexalgorithmus für gewisse Problemklassen geben ein sehr viel günstigeres Bild dieses Verfahrens. Dieses Bild spiegelt die auf 50-jährigem Umgang mit dem Simplexalgorithmus beruhende Erfahrung wider, dass dieser im Schnitt nur $Km$ Iterationen zur Lösung des anfangs beschriebenen linearen Optimierungsproblems benötigt, wobei $K$ eine relativ kleine Konstante ist. (Es werden unterschiedliche Werte zwischen 1 und 6 angegeben. Dantzig vermutete 1979 in einer Arbeit über den Khachian-Algorithmus, dass der Simplexalgorithmus nach durchschnittlich höchstens $m\ln(n)$ Iterationen eine Lösung von $(P)$ findet.) Derartige Ergebnisse bzw. Beobachtungen zur durchschnittlichen Iterationszahl eines Verfahrens sind z. B. auch für die Entscheidung darüber wichtig, ob man bei großen Problemen das ursprüngliche primale Problem oder das dazu gehörige duale Problem lösen sollte.

Nachdem durch Khachian geklärt worden war, dass das Problem der linearen Optimierung polynomiale Komplexität besitzt, wuchs die Hoffnung, dass man eines Tages auch effizientere Algorithmen als den Simplexalgorithmus zu seiner Lösung finden würde. Aufgrund der Eigenschaften des Simplexalgorithmus war zu erwarten, dass ein solches Verfahren Iterierte erzeugen muss, welche im Inneren der zulässigen Menge liegen. Die wesentliche Frage, die es in diesem Zusammenhang zu lösen galt, war dann die nach einer geeigneten Richtung, die von einem gegebenen Punkt im Inneren des Gebietes zu einem inneren Punkt zeigt, der eine „wesentlich bessere“ Näherung in Bezug auf eine Lösung des Problems darstellt.

Im Jahr 1984 schlug Karmarkar auf einer Tagung ein neues Verfahren vor, welches eine polynomiale Laufzeit mit einem günstigeren Worst-Case-Verhalten als das von Khachian besitzt. Diesem Verfahren liegt die folgende Beobachtung zugrunde: Angenommen das zulässige Gebiet eines LOPs ist ein beschränktes Polyeder, d. h. ein Polytop und man befindet sich in einem Punkt nahe dem Zentrum dieses Polytops, dann ermöglicht die zulässige Richtung steilsten Abstiegs einen relativ großen Schritt in Richtung einer Lösung des Problems, während sie für nicht zentral liegende Punkte schnell in die Nähe des Randes des Polytops führen und daher ungünstig sein kann. Denn kann man dem Rand nicht mehr entkommen und folgt man diesem bzw. muss man diesem aufgrund der festgelegten Richtungswahl (im Inneren des Polytops) in Richtung einer Lösung folgen, so kann dies im ungünstigsten Fall zu einem ähnlichen Verhalten führen, wie es Klee und Minty für den Simplexalgorithmus gezeigt hatten.

Karmarkar hatte deshalb die Idee, in jedem Schritt eine projektive Transformation zu verwenden, die das Polytop so auf sich selbst abbildet, dass die aktuelle (transformierte) Iterierte zentral bezüglich des durch die Transformation erzeugten Polytops liegt. Er schlug dann weiter vor, für das mittels dieser Transformation erzeugte Problem bzw., da dieses nichtlinear ist, für ein dazu äquivalentes lineares Problem die zulässige Richtung steilsten Abstiegs zu bestimmen und damit für dieses Problem eine neue Näherung zu erzeugen, die wiederum im Inneren des Polytops liegt. Diese mit der inversen Transformation in den Ursprungsraum zurück transformierte Näherung dient dann als nächste Iterierte für das ursprüngliche Problem. Karmarkar konnte Konvergenz und polynomiale Laufzeit für sein Verfahren beweisen.

Im Unterschied zum Simplexalgorithmus arbeitet das Verfahren von Karmarkar also mit inneren Punkten des zulässigen Gebietes. Karmarkar’s Behauptung, dass sein Verfahren, welches ja eine polynomiale Laufzeit aufweist, überdies bei den größten von ihm gerechneten Problemen bis zu 50 mal schneller sei als der Simplexalgorithmus, löste große Erwartungen aus. Die Tatsache, dass er - vermutlich aus vertraglichen Gründen gegenüber seinem Arbeitgeber, den AT&T Bell Laboratories - nicht bereit war (s. Combinatorica 5, 1984), rechnerische Details anzugeben (ein guter Code zur Lösung linearer Optimierungsprobleme bringt heute sehr viel Geld ein), verursachte jedoch gleichzeitig große Verärgerung und Diskussionen bis hinein in die Tagespresse (wie z.B. Schrijver im CWI Newsletter 8, 1985, beschreibt). Mit der Zeit wurden aber die Details des Karmarkar-Verfahrens bekannt und stellte es sich heraus, dass es im Einzelfall zwar erheblich schneller als der Simplexalgorithmus sein, ein solches Verhalten aber nicht allgemein angenommen werden kann.

Karmarkar geht für seinen Algorithmus von einer speziellen Gestalt eines linearen Optimierungsproblems aus. Die Überführung eines gegebenen LOPs in diese Gestalt kann ein Problem erheblich vergrößern. Außerdem ist die von Karmarkar - im Rahmen der linearen Optimierung verwendete - projektive Transformation eine nichtlineare Abbildung, und hat diese Nachteile. Das Karmarkar-Verfahren selbst hat heute keine praktische Bedeutung mehr. Wesentlich war jedoch, wie Karmarkar die polynomiale Komplexität seines Verfahrens bewiesen hatte. Er hatte nämlich eine sog. Potentialfunktion verwendet und gezeigt, dass die polynomiale Komplexität seines Verfahrens folgt, wenn der Wert dieser Funktion in jeder Iteration um den gleichen Wert abnimmt (und damit gegen $-\infty$ strebt). Also lag es nahe, direkt oder indirekt eine solche Potentialfunktion zu verwenden und Verfahren zu konstruieren, deren Iterierte die entsprechende Funktion in jeder Iteration um denselben Wert reduzieren.

Die Auseinandersetzungen mit dem Karmarkar-Verfahren erwiesen sich für die Optimierung als äußerst fruchtbar. Sie lösten eine wahre Flut von Publikationen aus und führten schließlich zu einer großen Zahl neuer Verfahren zur Lösung linearer und konvexer Optimierungsprobleme, den sog. Innere-Punkte-Verfahren (interior-point methods). Dabei haben sich inzwischen einige unterschiedliche Verfahrensklassen herausgebildet, wobei sich jede dieser Klassen wiederum in primale, duale und primal-duale Verfahren untergliedert, in Bezug darauf, ob diese primal, dual oder primal und dual zulässige Iterierte erzeugen.

Zunächst einmal gibt es die Klasse der Affine-Scaling-Verfahren, die wie das Karmarkar-Verfahren auf einer Transformation beruhen, welche aber im Vergleich zu diesem linear ist. Überdies geht beispielsweise ein primaler Affine-Scaling-Algorithmus direkt von der Standard-Normalform eines LOPs aus. Ein erster solcher primaler Affine-Scaling-Algorithmus wurde 1986 unabhängig voneinander von Barnes und von Vanderbei, Meketon und Freedman entwickelt und stellt wohl die einfachste Innere-Punkte-Methode überhaupt dar. (Wie man später feststellte, war dieser Algorithmus bereits 1967 von Dikin vorgeschlagen worden). Leider zeigten aber Experimente (ein theoretischer Beweis steht noch aus), dass der Algorithmus, wenn er auf das Problem $(P)$ angewandt und in der Nähe einer Ecke von $Z_{P}$ gestartet wird, dem Rand von $Z_{P}$ im Inneren von $Z_{P}$ folgen und somit im Extremfall wie der Simplexalgorithmus eine exponentielle Anzahl von Iterationen zur Lösung von $(P)$ benötigen kann.

Die Potentialreduktionsverfahren (potential reduction methods) verwenden direkt eine Potentialfunktion für ein LOP und vermindern deren Wert in jedem Schritt wenigstens um eine feste Größe. Diese Verfahren gehen auf Ye, 1988, zurück. Für sie kann üblicherweise polynomiale Laufzeit nachgewiesen werden.

Es wurde erkannt, dass die Funktionswerte einer Potentialfunktion in jeder Iteration am ehesten dann um dieselbe Größe verringert werden können, wenn die Iterierten des zugehörigen Verfahrens in gewissem Sinne zentral bezüglich des Inneren des jeweiligen zulässigen Gebietes liegen. Diese Beobachtung hat zu einer weiteren, der wohl wichtigsten Klasse von Innere-Punkte-Verfahren geführt, den Pfadverfolgungsmethoden (path following methods). Solche Methoden (die erste geht auf Meggido, 1989, zurück) erzeugen Iterierte, welche näherungsweise dem sog. zentralen Pfad im Inneren des betreffenden zulässigen Gebietes folgen. Bei den wichtigen primal-dualen Pfadverfolgungsmethoden ist ein Punkt auf dem zugehörigen Pfad als eindeutige Lösung eines schwach nichtlinearen Gleichungssystems gegeben, so dass bei diesen Verfahren das Newton-Verfahren zur Lösung nichtlinearer Gleichungssysteme eine entscheidende Rolle spielt.

Pfadverfolgungsverfahren vereinigen ein ausgezeichnetes theoretisches und praktisches Verhalten und sind daher gegenwärtig wohl die Verfahren, welche man insbesondere zur Lösung sehr großer linearer Optimierungsprobleme verwenden sollte. Ihre Grundlagen gehen auf Arbeiten von Frisch, 1956, und Fiacco und McCormick ([FiMcCo68]) zurück.

Innere-Punkte-Verfahren können also dem Simplexalgorithmus bei großen linearen Problemen deutlich überlegen sein. Solche Verfahren wurden inzwischen auch für die Lösung konvexer Optimierungsprobleme entwickelt, und es gibt bereits zahlreiche Vorschläge, wie die ihnen zugrunde liegenden Ideen für die Lösung nichtlinearer Probleme genutzt werden können. Man kann jetzt schon sagen, dass diese Entwicklung auch das Gebiet der numerischen nichtlinearen Optimierung, welches schon weitgehend abgeschlossen zu sein schien, erheblich bereichert hat und möglicherweise auch noch weiter bereichern wird.

Das Karmarkar-Verfahren, ein primaler Affine-Scaling-Algorithmus und ein Potentialreduktionsverfahren werden ausführlicher z. B. in [Ree01] diskutiert. Im Folgenden wollen wir hier nur auf die wichtigen primal-dualen Pfadverfolgungsmethoden näher eingehen.

6.2 Definitionen

Wir gehen im Folgenden aus von dem linearen Optimierungsproblem in Normalform

${\begin{array}{lll}(P):&{\text{Minimiere}}&c^{T}x\\&{\text{u. d. N.}}&Ax=b,\\&&x\geq 0.\end{array}}$

Den zulässigen Bereich von Problem $(P)$ bezeichnen wir wieder mit

$Z_{P}:=\{x\in \mathbb {R} ^{n}{\big |}Ax=b,\quad x\geq 0\}$

und sein Inneres (vgl. Abschnitt 4.4) mit

$Z_{P}^{o}:=\{x\in \mathbb {R} ^{n}{\big |}Ax=b,\quad x>0\}.$

Das zu $(P)$ duale Problem lautet

${\begin{array}{lll}(D):&{\text{Maximiere}}&b^{T}y\\&{\text{u. d. N.}}&A^{T}y+s=c,\\&&s\geq 0.\end{array}}$

Es besitzt den zulässigen Bereich

$Z_{D}:=\left\{(y,s)\in \mathbb {R} ^{m+n}{\big |}A^{T}y+s=c,\quad s\geq 0\right\}$

mit Innerem

$Z_{D}^{o}:=\left\{(y,s)\in \mathbb {R} ^{m+n}{\big |}A^{T}y+s=c,\quad s>0\right\}.$

Für $x\in \mathbb {R} _{+}^{n}$ definieren wir durch

$X:=\operatorname {diag} (x_{1},...,x_{n})$

eine Diagonalmatrix, welche für $x>0$ positiv definit ist und die positiv definite Inverse

$X^{-1}=\operatorname {diag} \left({\frac {1}{x_{1}}},...,{\frac {1}{x_{n}}}\right)$

besitzt. Entsprechend sind $S,X_{k}$ und $S_{k}$ die durch $s,x^{k},s^{k}\geq 0$ erzeugten Diagonalmatrizen.

Schließlich machen wir grundsätzlich die für primal-duale Innere-Punkte-Verfahren typischen Annahmen

(A1) $\operatorname {Rang} (A)=m$,

(A2) $Z_{P}^{o}\neq \emptyset$,

(A3) $Z_{D}^{o}\neq \emptyset$.

Die Annahmen (A2) und (A3) implizieren nach dem starken Dualitätssatz insbesondere die Lösbarkeit von $(P)$ und $(D)$. In diesem Zusammenhang ist auch folgendes Resultat interessant (s. [Wri97], S. 26 ff.):

Satz 6.2

Es seien $Z_{P}\neq \emptyset$ und $Z_{D}\neq \emptyset .Danngilt::''(i)<math>Z_{D}^{o}\neq \emptyset \Leftrightarrow$ Die Lösungsmenge von $(P)$ ist nichtleer und beschränkt.

(ii) $Z_{P}^{o}\neq \emptyset \Leftrightarrow \{s\in \mathbb {R} ^{n}{\big |}(y,s)\in Z_{D}\ f{\ddot {u}}r\ ein\ y\in \mathbb {R} ^{m}\}$ ist nichtleer und beschränkt.

6.3 Existenz des zentralen Pfades

Da bekannt ist, wie lineare Gleichungssystem zu lösen sind, verursacht die Ungleichungsnebenbedingung $x\geq 0$ die Hauptschwierigkeit bei der Lösung des LOPs $(P)$. Eine Idee ist es daher zu versuchen, diese Bedingung als Nebenbedingung zu entfernen und sie in einer geeigneten Form wenigstens näherungsweise in der Zielfunktion mit zu berücksichtigen. Ähnlich kann man mit der Bedingung $s\geq 0$ in dem dualen Problem $(D)$ verfahren.

Man betrachtet dazu die von einem Parameter $\tau >0$ abhängenden Hilfsprobleme

${\begin{array}{lll}(P_{\tau }):&{\text{Minimiere}}&V_{\tau }(x):=c^{T}x-\tau \sum \limits _{i=1}^{n}\ln(x_{i})\\&{\text{u. d. N.}}&Ax=b,\\&&x>0\end{array}}$

und

${\begin{array}{lll}(D_{\tau }):&{\text{Maximiere}}&W_{\tau }(y,s):=b^{T}y+\tau \sum \limits _{i=1}^{n}\ln(s_{i})\\&{\text{u. d. N.}}&ATy+s=c,\\&&s>0.\end{array}}$

Die Bedingungen $x>0$ und $s>0$ müssen wir in die Problemformulierungen mit aufnehmen, da die jeweilige Zielfunktion für andere Werte nicht definiert ist. Sie sind aber im Rahmen von Innere-Punkte-Verfahren einfach zu behandeln. Die Probleme $(P_{\tau })$ und $(D_{\tau })$ besitzen also gerade $Z_{P}^{o}$ und $Z_{D}^{o}$ als zulässige Gebiete, welche nach unseren Annahmen (A2) und (A3) nichtleer sind.

Die Zielfunktion von $(P_{\tau })$ gewinnt man dabei aus der Zielfunktion $\sum _{i=1}^{n}c_{i}x_{i}$ von $(P)$, indem man zu letzterer die logarithmische Barrierefunktion

$-\sum _{i=1}^{n}\ln(x_{i})$

mit Barriereparameter $\tau$ hinzu addiert. Für $\tau >0$ erzeugt diese im Inneren des zulässigen Gebietes $Z_{P}$ von $(P)$ eine Barriere, welche verhindert, dass man bei der Lösung von $(P_{\tau })$ dem Rand von $\mathbb {R} _{+}^{n}$ zu nahe kommt. Denn für $x_{i}\to 0+$ hat man $c_{i}x_{i}-\tau \ln(x_{i})\to +\infty$, so dass eine Lösung $x_{\tau }^{*}$ von $(P_{\tau })$, sofern eine solche existiert, nicht „zu nahe“ am Rand von $\mathbb {R} _{+}^{n}$ liegen wird. Dabei ist zu vermuten, dass der Barriereterm bei der Lösung von $(P_{\tau })$ für $\tau \to 0+$ eine immer kleinere Rolle spielt und eine Folge $\{x_{\tau }^{*}\}$, sofern eine solche definiert ist, für $\tau \to 0+$ gegen eine Lösung von $(P)$ konvergiert. Analoges lässt sich für das Problem $(D_{\tau })$ sagen.

Es gilt nun

${\begin{matrix}\nabla V_{\tau }(x)=c-\tau X^{-1}e,&\nabla ^{2}V_{\tau }(x)=\tau X^{-2},\\\nabla W_{\tau }(y,s)={\begin{pmatrix}b\\\tau S^{-1}e\end{pmatrix}},&\nabla ^{2}W_{\tau }(y,s)={\begin{pmatrix}0&0\\0&-\tau S^{-2}\end{pmatrix}}.\end{matrix}}$

Für $x>0$ ist die Matrix $\nabla ^{2}V_{\tau }(x)$ eine positiv definite und für $s>0$ die Matrix $\nabla ^{2}W_{\tau }(y,s)$ eine negativ semidefinite Diagonalmatrix. Also ist $V_{\tau }$ strikt konvex auf $Z_{P}^{o}$ und $W_{\tau }$ konkav auf $Z_{D}^{o}$. Insbesondere hat daher Problem $(P_{\tau })$ nach Satz 2.37 höchstens eine Lösung $x_{\tau }^{*}$.

Die Probleme $(P)$ und $(D)$ besitzen beide das System $(S)$ aus Abschnitt 4.6.1 als notwendige und hinreichende Optimalitätsbedingungen (siehe auch das System $(S^{*})$ in Abschnitt 4.6.3). Wir wollen als nächstes zeigen, dass auch $(P_{\tau })$ und $(D_{\tau })$ durch ein gemeinsames System von Optimalitätsbedingungen charakterisierbar sind. Dazu überlege man sich zunächst:

Bemerkung 6.3

Man überlegt sich leicht: Enthält das Problem $({\mathcal {P}})$ aus Abschnitt 3.2 zusätzlich die Nebenbedingung $x\in U$, wobei $U\subseteq \mathbb {R} ^{n}$ eine offene Menge ist, so bleibt die Aussage des Korollars 3.14 gültig, wenn man diese Nebenbedingung mit in die KKT-Bedingungen aufnimmt.

Unter Verwendung von Bemerkung 6.3 schließen wir für $U:=\{x\in \mathbb {R} ^{n}{\big |}x>0\}$, dass $x$ genau dann Lösung von $(P_{\tau })$ ist, wenn $(x,\lambda )$ für ein $\lambda$ das folgende System löst:

(6.2) $c-\tau X^{-1}e+A^{T}\lambda =0,$

(6.3) $Ax=b,$

(6.4) $x>0.$

Wir setzen nun $y:=-\lambda$ und $s:=\tau X^{-1}e$, wobei wir für letztere Beziehung mehrere äquivalente Schreibweisen angeben können:

(6.5) $s=\tau X^{-1}e\Leftrightarrow Xs=\tau e\Leftrightarrow XSe=\tau e\Leftrightarrow x_{i}s_{i}=\tau ~(i=1,\ldots ,n)\Leftrightarrow Sx=\tau e\Leftrightarrow x=\tau S^{-1}e.$

Offenbar gilt somit $x>0$ genau dann, wenn $s>0$ ist, so dass wir die (redundante) Bedingung $s>0$ auch noch mit in das obige System aufnehmen können. Das System (6.2)–(6.4) besitzt also genau dann eine Lösung $(x,\lambda )$, wenn das System

${\begin{array}{lr}(S_{\tau }^{*}):&Ax=b,\\&A^{T}y+s=c,\\&XSe=\tau e,\\&x>0,s>0\end{array}}$

eine Lösung $(x,y,s)$ hat.

Problem $(D_{\tau })$ können wir in ein konvexes Minimierungsproblem mit konvexer Zielfunktion $-W_{\tau }(y,s)$ umschreiben. Gemäß Bemerkung 6.3 löst demnach $(y,s)$ das Problem $(D_{\tau })$ genau dann, wenn $(y,s,\rho )$ für ein $\rho \in \mathbb {R} ^{n}$ den folgenden Bedingungen genügt:

(6.6) $-b+A\rho =0,$

(6.7) $-\tau S^{-1}e+\rho =0,$

(6.8) $A^{T}y+s=c,$

(6.9) $s>0.$

Fassen wir (6.6) und (6.7) zu

(6.10) $\tau AS^{-1}e=b$

zusammen, setzen wir $x:=\tau S^{-1}e$, was mit $s>0$ auch $x>0$ impliziert und nutzen wir (6.5) aus, so geht das System (6.8)–(6.10) ebenfalls in das System $(S_{\tau }^{*})$ über.

Insbesondere existiert also zu einer Lösung $x_{\tau }^{*}$ von $(P_{\tau })$, sofern es ein solches $x_{\tau }^{*}$ gibt, eine Lösung $(y,s)$ von $(D_{\tau })$, so dass $(x_{\tau }^{*},y,s)$ eine Lösung von $(S_{\tau }^{*})$ ist. Wir werden mit Satz 6.5 zeigen, dass $(P_{\tau })$ und $(D_{\tau })$ eindeutige Lösungen haben und dass somit gilt:

(6.11) $x_{\tau }^{*}$ löst $(P_{\tau }),(y_{\tau }^{*},s_{\tau }^{*})$ löst $(D_{\tau })\Leftrightarrow (x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*})$ löst $(S_{\tau }^{*})$.

Satz 6.5 im voraus verwendend können wir also unsere Überlegungen zu folgendem Ergebnis zusammenfassen:

Satz 6.4

Sei $\tau >0$.

(i) Es gilt:

$(P_{\tau })$ ist lösbar $\Longleftrightarrow (S_{\tau }^{*})$ ist lösbar $\Longleftrightarrow (D_{\tau })$ ist lösbar.

(ii) $x_{\tau }^{*}$ und $(y_{\tau }^{*},s_{\tau }^{*})$ sind genau dann Lösungen von $(P_{\tau })$ und $(D_{\tau })$, wenn $(x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*})$ eine Lösung von $(S_{\tau }^{*})$ ist.

Das System $(S_{\tau }^{*})$ ähnelt offenbar dem System $(S^{*})$ aus Abschnitt 4.6.3. Während $(S^{*}),(P)$ und $(D)$ keine eindeutigen Lösungen haben müssen, können wir beweisen, dass $(S_{\tau }^{*}),(P_{\tau })$ und $(D_{\tau })$ unter den Voraussetzungen (A1)–(A3) eindeutig lösbar sind. Insbesondere hat man dann damit, dass $(P_{\tau }),(D_{\tau })$ und $(S_{\tau }^{*})$ Lösungen besitzen.

Satz 6.5

Die Probleme $(P_{\tau })$ und $(D_{\tau })$ und das System $(S_{\tau }^{*})$ besitzen für jedes $\tau >0$ eindeutige Lösungen.

Beweis.

Wir zeigen zunächst, dass $(P_{\tau })$ eine Lösung hat. Dazu sei $({\bar {y}},{\bar {s}})\in Z_{D}^{o}$ fest gewählt. Für $x\in Z_{P}^{o}$ gilt dann nach dem schwachen Dualitätssatz

$0\leq x^{T}{\bar {s}}=c^{T}x-b^{T}{\bar {y}}.$

Damit schließen wir

(6.12) $V_{\tau }(x)=c^{T}x-\tau \sum _{i=1}^{n}\ln(x_{i})=\sum _{i=1}^{n}(x_{i}{\bar {s}}_{i}-\tau \ln(x_{i}))+b^{T}{\bar {y}}$

$=\sum _{i=1}^{n}\left(x_{i}{\bar {s}}_{i}-\tau \ln(x_{i})-\tau +\tau \ln({\frac {\tau }{{\bar {s}}_{i}}})\right)+b^{T}{\bar {y}}+n\tau -\tau \sum _{i=1}^{n}\ln({\frac {\tau }{{\bar {s}}_{i}}}).$

Sei nun

$g_{i}(t):={\bar {s}}_{i}t-\tau \ln(t)-\tau +\tau \ln({\frac {\tau }{{\bar {s}}_{i}}}).$

Dann folgt

$g'_{i}(t)={\bar {s}}_{i}-{\frac {\tau }{t}},\quad g''_{i}(t)={\frac {\tau }{t^{2}}}.$

Also ist $g_{i}$ strikt konvex für $t>0$ und besitzt $g_{i}$ den eindeutigen Minimalpunkt ${\hat {t}}:=\tau /{\bar {s}}_{i}>0$ mit $g_{i}({\hat {t}})=0$, was $g_{i}(t)\geq 0$ für alle $t>0$ impliziert. Schließlich hat man

$g_{i}(t)\to \infty \Leftrightarrow t\to 0$ oder $t\to \infty$.

Für ein $x^{0}\in Z_{P}^{o}$ definieren wir jetzt

$\gamma _{0}:=\sum _{i=1}^{n}g_{i}(x_{i}^{0})\geq 0.$

Aus den Eigenschaften von $g_{i}$ folgt die Existenz von $\alpha _{i}>0,\beta _{i}\geq \alpha _{i}$, so dass gilt:

$\{t>0{\big |}g_{i}(t)\leq \gamma _{0}\}=[\alpha _{i},\beta _{i}].$

Für die Niveaumenge

$N_{\tau }(x^{0}):=\left\{x\in Z_{P}^{o}{\big |}V_{\tau }(x)\leq V_{\tau }(x^{0})\right\}=\left\{x\in Z_{P}^{o}{\big |}\sum _{i=1}^{n}g_{i}(x_{i})\leq \sum _{i=1}^{n}g_{i}(x_{i}^{0})\right\}$

können wir daher schließen:

$N_{\tau }(x^{0})\subseteq \{x\in Z_{P}^{o}{\big |}g_{i}(x_{i})\leq \gamma _{0}~(i=1,\ldots ,n)\}\subseteq [\alpha _{1},\beta _{1}]\times \ldots \times [\alpha _{n},\beta _{n}].$

Letzteres zeigt, dass $N_{\tau }(x^{0})$ beschränkt ist.

Für den Nachweis der Abgeschlossenheit von $N_{\tau }(x^{0})$ gehen wir von einer Folge $\{x^{k}\}$ mit $x^{k}\in N_{\tau }(x^{0})$ und $\lim _{k\to \infty }x^{k}={\hat {x}}$ aus. Offenbar ist dann ${\hat {x}}\in Z_{P}$. Es ist sogar ${\hat {x}}\in Z_{P}^{o}$. Denn wäre ${\hat {x}}_{\ell }=0$ für ein $\ell \in \{1,\ldots ,n\}$, so hätte man $\lim _{k\to \infty }x_{\ell }^{k}={\hat {x}}_{\ell }=0$ und damit den Widerspruch

$V_{\tau }(x^{0})\geq V_{\tau }(x^{k})=c^{T}x^{k}-\tau \sum _{i=1}^{n}\ln(x_{i}^{k})\to +\infty \quad (k\to \infty ).$

Stetigkeit von $V_{\tau }$ auf $Z_{P}^{o}$ impliziert schließlich ${\hat {x}}\in N_{\tau }(x^{0})$. Aus den Sätzen 2.37 und 2.40 kann man daher folgern, dass das Problem $(P_{\tau })$ genau eine Lösung besitzt.

Da $(P_{\tau })$ lösbar ist, folgt nun nach dem vor Satz 6.4 Gezeigten die Existenz einer Lösung $({\tilde {x}}_{\tau },{\tilde {y}}_{\tau },{\tilde {s}}_{\tau })$ von $(S_{\tau }^{*})$, wobei ${\tilde {x}}_{\tau }:=x_{\tau }^{*}$ die eindeutige Lösung von $(P_{\tau })$ ist. Die Bedingung $XSe=\tau e$ legt ${\tilde {s}}_{\tau }$ als ${\tilde {s}}_{\tau }:=\tau {\tilde {X}}_{\tau }^{-1}e$ mit

${\tilde {X}}_{\tau }^{-1}:=\operatorname {diag} \left({\frac {1}{{\tilde {x}}_{\tau ,1}}},\ldots ,{\frac {1}{{\tilde {x}}_{\tau ,n}}}\right)$

auf eindeutige Weise fest. Weiter ist ${\tilde {y}}_{\tau }$ Lösung der Gleichung $A^{T}y+{\tilde {s}}_{\tau }=c$ und damit der Gleichung $AA^{T}y=A(c-{\tilde {s}}_{\tau })$, welche wegen (A1) (vgl. Lemma 2.21) eindeutig lösbar ist mit

${\tilde {y}}_{\tau }:=(AA^{T})^{-1}A(c-{\tilde {s}}_{\tau }).$

Also hat $(S_{\tau }^{*})$ genau eine Lösung.

Schließlich gilt nach dem vor Satz 6.4 Gezeigten, weil $(P_{\tau })$ eine Lösung hat, dass dies auch $(D_{\tau })$ tut. Wie wir dort hergeleitet haben, gibt es insbesondere zu Lösungen $(y_{\tau }^{1},s_{\tau }^{1})$ und $(y_{\tau }^{2},s_{\tau }^{2})$ von $(D_{\tau })$ Punkte $x_{\tau }^{1},x_{\tau }^{2}\in \mathbb {R} ^{n}$, so dass $(x_{\tau }^{1},y_{\tau }^{1},s_{\tau }^{1})$ und $(x_{\tau }^{2},y_{\tau }^{2},s_{\tau }^{2})$ das System $(S_{\tau }^{*})$ lösen. Die eindeutige Lösbarkeit von $(S_{\tau }^{*})$ impliziert also die von $(D_{\tau })$.

q.e.d.

Die somit wohldefinierte Menge

$C:=\{(x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*}){\big |}(x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*}){\mbox{ löst }}(S_{\tau }^{*}){\mbox{ für ein }}\tau >0\}$

nennt man den primal-dualen zentralen Pfad (primal-dual central path). Entsprechend bezeichnet man die Menge

$C_{P}:=\{x_{\tau }^{*}{\big |}(x_{\tau }^{*}){\mbox{ löst }}(P_{\tau }){\mbox{ für ein }}\tau >0\}=\{x_{\tau }^{*}{\big |}(x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*}){\mbox{ löst }}(S_{\tau }^{*}){\text{ mit einem }}(y_{\tau }^{*},s_{\tau }^{*}){\mbox{ für ein }}\tau >0\}$

als den primalen zentralen Pfad (primal central path).

Wenn für den primalen oder primal-dualen zentralen Pfad Konvergenz

$\lim _{\tau \to 0+}(x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*})=(x^{*},y^{*},s^{*})$

gezeigt werden kann, so muss der Grenzwert $(x^{*},y^{*},s^{*})$, wie man sofort durch Grenzwertbildung in $(S_{\tau }^{*})$ erkennt, Lösung des Systems

${\begin{array}{lr}(S):&Ax=b,\\&A^{T}y+s=c,\\&XSe=0,\\&x\geq 0,s\geq 0\end{array}}$

sein und müssen damit $x^{*}$ und $(y^{*},s^{*})$ das Problem $(P)$ bzw. $(D)$ lösen. Wir können also $(S)$ als Grenzwertproblem zur Folge der Probleme $(S_{\tau }^{*})$ für $\tau \to 0+$ ansehen. Es deutet sich damit an, dass Lösungen von $(S_{\tau }^{*})$ für $\tau \to 0+$ zu Lösungen von $(P)$ und $(D)$ führen.

Wie man durch Betrachtung der jeweiligen Optimalitätsbedingungen erkennt, ist das Problem $(P_{\tau })$ äquivalent mit dem Problem, das man erhält, wenn man die Zielfunktion $V_{\tau }$ von $(P_{\tau })$ durch $\tau$ dividiert. Diesem äquivalenten Problem kann man für $\tau =\infty$ formal das Problem

(6.13) ${\begin{array}{lr}{\text{Minimiere}}&-\sum \limits _{i=1}^{n}\ln(x_{i})\\{\text{u. d. N.}}&Ax=b,\\&x>0\end{array}}$

zuordnen, das in der Theorie der Inneren-Punkte-Verfahren eine besondere Bedeutung hat. Im Fall, dass Problem (6.13) einen eindeutigen Minimalpunkt besitzt, bezeichnet man diesen als das analytische Zentrum (analytic center) von $Z_{P}$. Es gilt z. B. (man beachte, dass wir Satz 6.5 unter den Voraussetzungen (A1), (A2) und (A3) für ein gegebenes $c$ bewiesen hatten):

Satz 6.6

Es sei $Z_{P}$ beschränkt. Dann besitzt Problem (6.13) eine eindeutige Lösung.

Beweis.

Um Satz 6.5 für $c:=0$ und $\tau :=1$ anwenden zu können, müssen wir das Erfülltsein von (A3) für $c:=0$ zeigen. Wir bezeichnen die Probleme $(P)$ und $(D)$ für den Fall $c:=0$ mit $({\bar {P}})$ und $({\bar {D}})$ und ihre zulässigen Bereiche und deren Innere mit $Z_{\bar {P}}=Z_{P},Z_{\bar {P}}^{o}=Z_{P}^{o}$ sowie $Z_{\bar {D}},Z_{\bar {D}}^{o}$. Da die Lösungsmenge von $({\bar {P}})$ nichtleer, gleich $Z_{P}$ und beschränkt ist, hat man nach den Sätzen 4.17 und 6.2, dass auch $({\bar {D}})$ eine Lösung besitzt und auch $Z_{\bar {D}}^{o}\neq \emptyset$ gilt.

q.e.d.

Wir geben dazu ein Beispiel.

Beispiel 6.7

Wir betrachten das folgende Problem:

(6.14) ${\begin{array}{ll}{\text{Minimiere}}&x_{2}\\{\text{u. d. N.}}&x_{1}+x_{2}+x_{3}=1,\\&(x_{1},x_{2},x_{3})^{T}\geq 0\end{array}}$

und dazu das $(S_{\tau }^{*})$ entsprechende KKT-System

$x_{1}+x_{2}+x_{3}=1,$

${\begin{pmatrix}1\\1\\1\end{pmatrix}}y+{\begin{pmatrix}s_{1}\\s_{2}\\s_{3}\end{pmatrix}}={\begin{pmatrix}0\\1\\0\end{pmatrix}},$

$x_{i}s_{i}=\tau \quad (i=1,2,3),$

$(x_{1},x_{2},x_{3})^{T}>0,$

$(s_{1},s_{2},s_{3})^{T}>0.$

Auflösung der dritten Zeile nach den $s_{i}$ und Einsetzung in die zweite Zeile führt auf das System

(6.15) $x_{1}+x_{2}+x_{3}=1,$

(6.16) ${\begin{pmatrix}1\\1\\1\end{pmatrix}}y+\tau {\begin{pmatrix}1/x_{1}\\1/x_{2}\\1/x_{3}\end{pmatrix}}={\begin{pmatrix}0\\1\\0\end{pmatrix}},$

(6.17) $(x_{1},x_{2},x_{3})^{T}>0.$

Aus der ersten und dritten Zeile von (6.16) erhält man $y=-{\frac {\tau }{x_{1}}}=-{\frac {\tau }{x_{3}}}$, also wegen $\tau >0$ insbesondere $x_{1}=x_{3}$. Gleichung (6.15) nach $x_{2}$ aufgelöst und in die zweite Zeile von (6.16) eingesetzt ergibt dann $y=1-{\frac {\tau }{1-2x_{1}}}$. Es folgt

$1-{\frac {\tau }{1-2x_{1}}}=-{\frac {\tau }{x_{1}}}\Leftrightarrow x_{1}(1-2x_{1}-\tau )=-\tau (1-2x_{1})\Leftrightarrow 2x_{1}^{2}-(1-3\tau )x_{1}-\tau =0.$

Wegen $x>0$ gelangt man mit letzterer Gleichung zu:

$x_{\tau ,1}^{*}:={\frac {1-3\tau +{\sqrt {9\tau ^{2}+2\tau +1}}}{4}}=x_{\tau ,3}^{*},\quad x_{\tau ,2}^{*}:={\frac {1-3\tau +{\sqrt {9\tau ^{2}+2\tau +1}}}{2}}.$

Die Folge $\{x_{\tau }^{*}\}$ konvergiert für $\tau \to 0+$ gegen den Vektor $x^{*}:=\left({\frac {1}{2}},0,{\frac {1}{2}}\right)^{T}$, der nach den obigen Überlegungen Lösung von $(P)$ ist.

Ersetzt man den Vektor $c$ auf der rechten Seite in (6.16) durch den Nullvektor und setzt man $\tau =1$, so erhält man aus (6.15)–(6.17) im Hinblick auf das hier (6.13) entsprechende Problem das KKT-System

$x_{1}+x_{2}+x_{3}=1,$

${\begin{pmatrix}1\\1\\1\end{pmatrix}}y+{\begin{pmatrix}1/x_{1}\\1/x_{2}\\1/x_{3}\end{pmatrix}}={\begin{pmatrix}0\\0\\0\end{pmatrix}},$

$(x_{1},x_{2},x_{3})^{T}>0.$

Dieses hat die Lösung $x_{a}:=\left({\frac {1}{3}},{\frac {1}{3}},{\frac {1}{3}}\right)^{T}$ und $y_{a}:=-3$. Der Vektor $x_{a}$ ist also das analytische Zentrum des zulässigen Bereichs von Problem (6.14). Man kann zeigen, dass $\lim _{\tau \to +\infty }x_{\tau }^{*}=x_{a}$ gilt. (Letzteres ist allgemein richtig, wenn $Z_{P}$ beschränkt ist; vgl. [Sai97], S. 225.)

Die Lösungsmenge von Problem (6.14) ist offenbar gegeben durch

${\mathcal {L}}_{B}:=\left\{(x_{1},0,x_{3})^{T}\in \mathbb {R} ^{3}{\big |}x_{1}+x_{3}=1,\quad x_{1}\geq 0,\quad x_{3}\geq 0\right\}.$

Das analytische Zentrum von ${\mathcal {L}}_{B}$ errechnet sich aus dem System

$x_{1}+x_{3}=1,$

${\begin{pmatrix}1\\1\end{pmatrix}}z+{\begin{pmatrix}1/x_{1}\\1/x_{3}\end{pmatrix}}={\begin{pmatrix}0\\0\end{pmatrix}},$

$x_{1}>0,\quad x_{2}=0,\quad x_{3}>0,$

welches $x^{*}$ mit $z^{*}:=-2$ als einzige Lösung besitzt. Also ist die Lösung $x^{*}$ von Problem (6.14), die man erhält, wenn man dem primalen Pfad folgt, das analytische Zentrum der Lösungsmenge des Problems. (Letzteres Resultat lässt sich allgemein beweisen; siehe [Sai97], S. 225.) Ferner ist offenbar

$x^{*}=\left({\frac {1}{2}},0,{\frac {1}{2}}\right)^{T},\quad (y^{*},s^{*}):=(0,(0,1,0)^{T})$

ein strikt komplementäres Lösungspaar für das Problem (6.14) und das dazu duale Problem, wie man durch Einsetzen von $x^{*}$ und $(y^{*},s^{*})$ in das oben aufgeführte, dem System $(S_{\tau }^{*})$ entsprechende KKT-System mit $\tau =0$ erkennt.

---

Analog kann man auch für $Z_{D}$ ein analytisches Zentrum definieren.

6.4 Ein allgemeiner Rahmen für primal-duale Verfahren

Der an Beispiel 6.7 gezeigte, aber hier nicht allgemein bewiesene Umstand, dass $\{x_{\tau }^{*}\}$ für $\tau \to 0+$ gegen eine Lösung von $(P)$ konvergiert (und zwar gegen das analytische Zentrum der Lösungsmenge von $(P)$, welche nach Satz 6.2 beschränkt ist), liefert die Motivation für eine auf Gonzaga ([Gon89]) zurückgehende primale Pfadverfolgungsmethode. Die Idee dabei ist es, ausgehend von Punkten $x^{0}\in Z_{P}^{o}$ und $(y^{0},s^{0})\in Z_{D}^{o}$ und einem $\tau _{0}>0$ in der $k$-ten Iteration einen Newton-Schritt für das Problem $(P_{\tau _{k}})$ durchzuführen und anschließend $\tau _{k}$ zu verkleinern. Ziel ist es, auf diese Weise dem primalen zentralen Pfad wenigstens näherungsweise zu folgen (daher: Pfadverfolgungsmethode).

Offenbar macht es keinen Sinn, $\tau _{k}>0$ im $k$-ten Schritt bei einer solchen Vorgehensweise festzuhalten und das Newton-Verfahren fortzuführen, da man ja an einer genauen Lösung von $(P_{\tau _{k}})$ für ein $\tau _{k}>0$ gar nicht interessiert ist. Bei dem beschriebenen Vorgehen ignoriert man die Nebenbedingungen $x>0$ und $s>0$ und sichert man deren Erfülltsein für $x^{k+1}$ und $s^{k+1}$ dadurch, dass man sich von $x^{k}>0$ und $s^{k}>0$ nur einen hinreichend kleinen Schritt weg bewegt. Auf diese primale Pfadverfolgungsmethode wollen wir hier aber nicht näher eingehen, sondern wir wollen uns gleich den effizienteren primal-dualen Methoden zuwenden.

In Abschnitt 6.3 haben wir gesehen, dass das System

${\begin{array}{lr}(S_{\tau }^{*}):&Ax=b,\\&A^{T}y+s=c,\\&XSe=\tau e,\\&x>0,\quad s>0,\end{array}}$

für festes $\tau >0$ eine eindeutige Lösung $(x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*})$ besitzt und dass $x_{\tau }^{*}\in Z_{P}^{o}$ und $(y_{\tau }^{*},s_{\tau }^{*})\in Z_{D}^{o}$ die Probleme $(P_{\tau })$ und $(D_{\tau })$ lösen. Die Menge aller solcher Lösungen für $\tau >0$ definiert den (primal-dualen) zentralen Pfad. Ferner haben wir festgestellt, dass im Fall der Konvergenz

(6.18) $\lim _{\tau \to 0+}(x_{\tau }^{*},y_{\tau }^{*},s_{\tau }^{*})=(x^{*},y^{*},s^{*})$

der Grenzwert $(x^{*},y^{*},s^{*})$ Lösung von $(S)$ ist und damit $x^{*}$ und $(y^{*},s^{*})$ Lösungen von $(P)$ bzw. $(D)$ sind. Also liegt es nahe, Verfahren zu konstruieren, die für eine Nullfolge $\{\tau _{k}\}$ Näherungen $(x^{k},y^{k},s^{k})$ von $(x_{\tau _{k}}^{*},y_{\tau _{k}}^{*},s_{\tau _{k}}^{*})$ erzeugen, so dass

(6.19) $\left\|(x^{k},y^{k},s^{k})-(x_{\tau _{k}}^{*},y_{\tau _{k}}^{*},s_{\tau _{k}}^{*})\right\|\to 0\quad (k\to \infty )$

und damit $(x^{k},y^{k},s^{k})\to (x^{*},y^{*},s^{*})$ für $k\to \infty$ gilt.

Im Hinblick auf die Lösung des Systems $(S_{\tau }^{*})$ für gegebenes $\tau$ betrachten wir das Gleichungssystem

(6.20) $Ax=b,$

(6.21) $A^{T}y+s=c,$

(6.22) $XSe=\tau e.$

In diesem sind die Gleichungen (6.20) und (6.21) linear und ist die Gleichung (6.22) schwach nichtlinear. Addiert man die einzelnen Zeilen der Gleichung (6.22), so folgt daraus für $\tau$ die Beziehung

$\tau ={\frac {x^{T}s}{n}}.$

Das System (6.20)–(6.22) schreiben wir jetzt in der Form

(6.23) $F_{\tau }(x,y,s):={\begin{bmatrix}Ax-b\\A^{T}y+s-c\\XSe-\tau e\end{bmatrix}}=0,$

so dass $F_{\tau }(x,y,s)$ für $x\in Z_{P}^{o}$ und $(y,s)\in Z_{D}^{o}$ durch

$F_{\tau }(x,y,s)={\begin{pmatrix}0\\0\\XSe-\tau e\end{pmatrix}}$

gegeben ist. Die Jacobi-Matrix zu $F_{\tau }$ lautet

$J(x,y,s):={\begin{pmatrix}A&0&0\\0&A^{T}&I\\S&0&X\end{pmatrix}}$

und ist offenbar von $y$ und $\tau$ unabhängig.

Das System $(S_{\tau }^{*})$ lässt sich für gegebenes $\tau$ im allgemeinen nicht exakt lösen und eine exakte Lösung von $(S_{\tau }^{*})$ wird ja auch, wie oben für ein primales Verfahren erläutert wurde, gar nicht benötigt. Die Idee bei den primal-dualen Pfadverfolgungsmethoden ist es, für gegebenes $\tau _{k}>0$ und einer durch $x^{k}\in Z_{P}^{o}$ und $(y^{k},s^{k})\in Z_{D}^{o}$ gegebenen Näherung $(x^{k},y^{k},s^{k})$ mit einer geeigneten Schrittweite einen Newton-Schritt in Richtung des auf dem zentralen Pfad liegenden Punktes $(x_{\tau _{k}}^{*},y_{\tau _{k}}^{*},s_{\tau _{k}}^{*})$ durchzuführen. Dabei soll $\{\tau _{k}\}$ gegen Null konvergieren und insbesondere durch die Schrittweitenwahl garantiert sein, dass die Glieder der Folge $\{(x^{k},y^{k},s^{k})\}$ in einer geeignet gewählten Umgebung des zentralen Pfades bleiben und die Konvergenz (6.19) folgt.

Das Newton-Verfahren wird bekanntlich zur Bestimmung einer Nullstelle einer einmal stetig differenzierbaren Funktion $f:\mathbb {R} ^{n}\to \mathbb {R} ^{n}$, also eines $x\in \mathbb {R} ^{n}$ mit $f(x)=0$, verwendet. Die Iterationsvorschrift des Newton-Verfahrens lautet

(6.24) $x^{k+1}:=x^{k}-J_{f}^{-1}(x^{k})f(x^{k}),$

wobei

$J_{f}(x^{k}):=\left({\frac {\partial f_{i}}{\partial x_{j}}}(x^{k})\right)_{i,j=1,\ldots ,n}$

die als nichtsingulär angenommene Jacobi-Matrix von $f$ in $x^{k}$ sei. Die in (6.24) vorkommende sog. Newton-Richtung

$d^{k}:=-J_{f}^{-1}(x^{k})f(x^{k})$

gewinnt man sinnvollerweise als Lösung des linearen Gleichungssystems

(6.25) $J_{f}(x^{k})d=-f(x^{k}).$

Denn die Berechnung der Inversen einer $(n\times n)$-Matrix, wie hier die der Matrix $J_{f}(x^{k})$, erfordert bekanntlich mehr Rechenoperationen.

Für $\tau _{k}$ und $(x^{k},y^{k},s^{k})$ mit $x^{k}\in Z_{P}^{o},(y^{k},s^{k})\in Z_{D}^{o}$ ergibt sich die Newton-Richtung für das System (6.23) als Lösung des linearen Gleichungssystems (vgl. (6.25))

$J(x^{k},y^{k},s^{k})(\Delta x,\Delta y,\Delta s)=-F_{\tau _{k}}(x^{k},y^{k},s^{k}),$

welches ausgeschrieben lautet:

(6.26) ${\begin{pmatrix}A&0&0\\0&A^{T}&I\\S_{k}&0&X_{k}\end{pmatrix}}{\begin{pmatrix}\Delta x\\\Delta y\\\Delta s\end{pmatrix}}={\begin{pmatrix}0\\0\\-X_{k}S_{k}e+\tau _{k}e\end{pmatrix}}.$

Lemma 6.8

Die Matrix $J(x^{k},y^{k},s^{k})$ in (6.26) ist für $x^{k}>0$ und $s^{k}>0$ nichtsingulär.

Beweis.

Mit $(u,v,w)\in \mathbb {R} ^{n+m+n}$ lautet das System $J(x^{k},y^{k},s^{k})(u,v,w)=0$ ausgeschrieben

(6.27) $Au=0,$

(6.28) $A^{T}v+w=0,$

(6.29) $S_{k}u+X_{k}w=0.$

Gleichung (6.28) liefert mit (6.27)

$0=u^{T}A^{T}v+u^{T}w=u^{T}w$

und Gleichung (6.29) impliziert

(6.30) $w=-X_{k}^{-1}S_{k}u.$

Also hat man

$u^{T}X_{k}^{-1}S_{k}u=0,$

woraus $u=0$ wegen der positiven Definitheit von $X_{k}^{-1}S_{k}$ folgt. Nach (6.30) ist damit $w=0$. Schließlich folgt mit (6.28) unter Verwendung unserer Standardvoraussetzung $\operatorname {Rang} (A)=m$, dass $v=0$ ist.

q.e.d.

Das System (6.26) besitzt also eine eindeutige Lösung $(\Delta x^{k},\Delta y^{k},\Delta s^{k})$. Geeigneterweise wird $\tau _{k}$ als Produkt der Form $\tau _{k}:=\sigma _{k}\mu _{k}$ gewählt, wobei

$\mu _{k}:={\frac {(x^{k})^{T}s^{k}}{n}}$

ein Dualitätsmaß, nämlich der arithmetische Mittelwert der Produkte $x_{i}^{k}s_{i}^{k}$ und $\sigma _{k}\in [0,1]$ ein Zentrierungsparameter ist. Für $\sigma _{k}:=1$ macht man offenbar einen Newton-Schritt in Richtung der Lösung des Systems $F_{\mu _{k}}(x,y,s)=0$, also in Richtung des zentralen Pfades, auf dem alle Produkte $x_{i}s_{i}$ denselben Wert $\mu _{k}$ haben. Man nennt daher eine solche Richtung auch zentrierende Richtung (centering direction). Eine solche Richtung bringt typischerweise nur wenig oder keinen Fortschritt bezüglich der Verkleinerung von $\mu _{k}$, lässt aber erwarten, dass man in der nächsten Iteration einen größeren Schritt bezüglich einer wohldefinierten Umgebung des zentralen Pfades gehen kann. Im anderen Extremfall $\sigma _{k}:=0$ macht man einen Newton-Schritt in Richtung einer Lösung des Systems (6.20)–(6.22) für $\tau =0$, was zwar im allgemeinen eine Verkleinerung von $\mu _{k}$ zur Folge hat, aber auch bedeuten kann, dass in der nächsten Iteration nur ein kleiner Schritt innerhalb der gewählten Umgebung des zentralen Pfades möglich ist. Die meisten Algorithmen verwenden daher Werte $\sigma _{k}\in (0,1)$, um einerseits Nähe zum zentralen Pfad zu erzwingen und andererseits $\mu _{k}$ zu verkleinern.

Diese Überlegungen führen nun zu folgendem Modellalgorithmus, durch den sich eine Reihe von primal-dualen Pfadverfolgungsmethoden mittels spezieller Parameterwahl beschreiben lassen.

Modellalgorithmus 6.9

(0) Wähle $x^{0}\in Z_{P}^{o},(y^{0},s^{0})\in Z_{D}^{o}$ und $\sigma _{0}\in [0,1]$. Setze $k:=0$.

(1) Berechne $\mu _{k}:=\left((x^{k})^{T}s^{k}\right)/n$ und bestimme die Lösung $(\Delta x^{k},\Delta y^{k},\Delta s^{k})$ des linearen Gleichungssystems

(6.31) ${\begin{pmatrix}A&0&0\\0&A^{T}&I\\S_{k}&0&X_{k}\end{pmatrix}}{\begin{pmatrix}\Delta x\\\Delta y\\\Delta s\end{pmatrix}}={\begin{pmatrix}0\\0\\-X_{k}S_{k}e+\sigma _{k}\mu _{k}e\end{pmatrix}}.$

(2) Bestimme eine geeignete Schrittweite $\alpha _{k}>0$ und ein $\sigma _{k+1}\in [0,1]$.

(3) Setze

$(x^{k+1},y^{k+1},s^{k+1}):=(x^{k},y^{k},s^{k})+\alpha _{k}(\Delta x^{k},\Delta y^{k},\Delta s^{k}).$

(4) Setze $k:=k+1$ und gehe nach (1).

Bemerkung 6.10

In der angegebenen Weise erzeugt der Modellalgorithmus 6.9 unendliche Folgen $\{x^{k}\}$ und $\{(y^{k},s^{k})\}$, was im Hinblick auf die folgenden Untersuchungen angenehmer ist. Für die Implementierung eines speziellen Verfahrens hat man zwischen Schritt (0) und (1) folgende Abbruchbedingung einzuführen:

Falls $0\leq \mu _{k}\leq \varepsilon$ ist, stop!

Im Schritt (0) ist dann außerdem noch

$\mu _{0}:={\frac {(x^{0})^{T}s^{0}}{n}}$

zu berechnen und $\varepsilon >0$ vorzugeben. Falls $\lim _{k\to \infty }\mu _{k}=0$ bzw. $\mu _{k}\leq \varepsilon$ für alle $k\geq K$ mit einem $K\in \mathbb {N} _{0}$ gezeigt werden kann, bricht dann das Verfahren nach endlich vielen Iterationen mit einem $\varepsilon$-optimalen Lösungspaar $x^{k}$ und $(y^{k},s^{k})$ von $(P)$ und $(D)$ ab (vgl. Korollar 4.19).

Man beachte, dass die Konvergenz von $\{\mu _{k}\}$ gegen 0 nicht die Konvergenz der Folgen $\{x^{k}\}$ und $\{(y^{k},s^{k})\}$ impliziert. Da man $x^{k}\in Z_{P}^{o}$ und $(y^{k},s^{k})\in Z_{D}^{o}$ für alle $k$ hat, können wir jedoch ähnlich wie (6.18) schließen: Gilt $\lim _{k\to \infty }\mu _{k}=0$, so ist jeder Häufungspunkt $(x^{*},y^{*},s^{*})$ der Folge $\{(x^{k},y^{k},s^{k})\}$ eine Lösung des Systems $(S)$ und sind damit $x^{*}$ und $(y^{*},s^{*})$ Lösungen von $(P)$ und $(D)$. Die Existenz von Häufungspunkten ist im Einzelfall zu beweisen.

Wir zeigen nun als erstes, dass die gewichtete Dualitätslücke $\mu _{k}$ bei geeigneter Wahl der Schrittweite und der Konstanten $\sigma _{k}$ in jedem Schritt verkleinert werden kann. Dazu machen wir von folgenden Notationen Gebrauch, wobei $\alpha >0$ eine beliebige Schrittweite ist:

(6.32) $\left(x^{k}(\alpha ),y^{k}(\alpha ),s^{k}(\alpha )\right):=(x^{k},y^{k},s^{k})+\alpha (\Delta x^{k},\Delta y^{k},\Delta s^{k}),$

(6.33) $\mu _{k}(\alpha ):={\frac {\left(x^{k}(\alpha )\right)^{T}s^{k}(\alpha )}{n}}.$

Lemma 6.11

Die Lösung $(\Delta x^{k},\Delta y^{k},\Delta s^{k})$ des linearen Gleichungssystems (6.31) besitzt die folgenden beiden Eigenschaften:

(i) $(\Delta x^{k})^{T}\Delta s^{k}=0$,

(ii) $\mu _{k}(\alpha )=[1-\alpha (1-\sigma _{k})]\mu _{k}$.

Beweis.

(i) Die eindeutige Lösung $(\Delta x^{k},\Delta y^{k},\Delta s^{k})$ des Gleichungssystems (6.31) genügt offenbar den Gleichungen

$A\Delta x^{k}=0,\quad A^{T}\Delta y^{k}+\Delta s^{k}=0.$

Multiplikation der zweiten Gleichung von links mit $(\Delta x^{k})^{T}$ und Anwendung der ersten Gleichung impliziert

$0=(\Delta x^{k})^{T}A^{T}\Delta y^{k}+(\Delta x^{k})^{T}\Delta s^{k}=(\Delta x^{k})^{T}\Delta s^{k}.$

(ii) Die dritte Blockzeile des Gleichungssystems (6.31) führt zu

$S_{k}\Delta x^{k}+X_{k}\Delta s^{k}=-X_{k}S_{k}e+\sigma _{k}\mu _{k}e.$

Summation der $n$ Komponenten dieser Gleichung liefert unter Verwendung der Definition von $\mu _{k}$

$(s^{k})^{T}\Delta x^{k}+(x^{k})^{T}\Delta s^{k}=-(1-\sigma _{k})(x^{k})^{T}s^{k}.$

Zusammen mit Aussage (i) ergibt sich daraus

$\left(x^{k}(\alpha )\right)^{T}s^{k}(\alpha )=\left(x^{k}+\alpha \Delta x^{k}\right)^{T}\left(s^{k}+\alpha \Delta s^{k}\right)$

$=(x^{k})^{T}s^{k}+\alpha \left((s^{k})^{T}\Delta x^{k}+(x^{k})^{T}\Delta s^{k}\right)+\alpha ^{2}(\Delta x^{k})^{T}\Delta s^{k}=(x^{k})^{T}s^{k}[1-\alpha (1-\sigma _{k})].$

q.e.d.

Der folgende Satz besagt nun, dass der Modellalgorithmus 6.9 eine polynomiale Laufzeit hat, wenn die Abnahme des Dualitätsmaßes in jedem Schritt in gewissem Sinne gleichmäßig ist und dieses Maß im Startpunkt nicht zu groß ist. Dabei bedeutet

$\lceil a\rceil :$ kleinste ganze Zahl größer oder gleich $a$.

Satz 6.12

Es sei $\varepsilon \in (0,1)$ beliebig gegeben und die Startvektoren $x^{0}\in Z_{P}^{o}$ und $(y^{0},s^{0})\in Z_{D}^{o}$ mögen mit einem $\kappa >0$ die Bedingung

(6.34) $\mu _{0}={\frac {(x^{0})^{T}s^{0}}{n}}\leq {\frac {1}{\varepsilon ^{\kappa }}}$

erfüllen. Genügen die durch die Iterierten $(x^{k},y^{k},s^{k}),k\in \mathbb {N} _{0}$ des Modellalgorithmus 6.9 erzeugten $\mu _{k}$ für gewisse Konstanten $\delta >0$ und $\omega >0$ mit $\delta /n^{\omega }<1$ der Ungleichung

(6.35) $0\leq \mu _{k+1}\leq \left(1-{\frac {\delta }{n^{\omega }}}\right)\mu _{k},$

so gilt

$0\leq \mu _{k}\leq \varepsilon ,\quad k\geq K$

für

$K:=\left\lceil (1+\kappa ){\frac {n^{\omega }}{\delta }}|\ln(\varepsilon )|\right\rceil .$

Beweis.

Aus (6.35) folgt für $k\in \mathbb {N} _{0}$

$\ln(\mu _{k+1})\leq \ln \left(1-{\frac {\delta }{n^{\omega }}}\right)+\ln(\mu _{k}).$

Mehrfache Anwendung dieser Formel liefert zusammen mit (6.34)

$\ln(\mu _{k})\leq k\ln \left(1-{\frac {\delta }{n^{\omega }}}\right)+\ln(\mu _{0})\leq k\ln \left(1-{\frac {\delta }{n^{\omega }}}\right)+\kappa \ln \left({\frac {1}{\varepsilon }}\right).$

Verwenden wir weiter die Abschätzung $1+a\leq e^{a}$ bzw. $\ln(1+a)\leq a$ für $a>-1$, so bekommen wir

$\ln(\mu _{k})\leq k\left(-{\frac {\delta }{n^{\omega }}}\right)+\kappa \ln \left({\frac {1}{\varepsilon }}\right).$

Das Konvergenzkriterium $\mu _{k}\leq \varepsilon$ ist also dann erfüllt, wenn

$k\left(-{\frac {\delta }{n^{\omega }}}\right)+\kappa \ln \left({\frac {1}{\varepsilon }}\right)\leq \ln(\varepsilon )$

gilt bzw., nach $k$ aufgelöst,

$k\geq (1+\kappa ){\frac {n^{\omega }}{\delta }}\ln \left({\frac {1}{\varepsilon }}\right)=(1+\kappa ){\frac {n^{\omega }}{\delta }}|\ln(\varepsilon )|.$

Daraus folgt die Behauptung.

q.e.d.

Die Hauptmühe im Zusammenhang mit Innere-Punkte-Verfahren macht der Nachweis des Erfülltseins der Bedingung (6.35). Wir werden im nächsten Abschnitt ein Verfahren vorstellen, für das man diese Bedingung verifizieren kann.

6.5 Ein zulässiges Verfahren

Wir beschreiben nun hintereinander zwei konkrete primal-duale Innere-Punkte-Verfahren zur Lösung linearer Optimierungsaufgaben. Die Idee dieser Verfahren besteht darin, dem primal-dualen zentralen Pfad in Richtung einer Lösung des Problems zu folgen. Bei den zulässigen Verfahren (feasible methods) müssen dafür der primale und duale Startvektor zulässig sein und sind damit alle Iterierten primal und dual zulässig. Bei den nichtzulässigen Verfahren (infeasible methods) kann man dagegen mit primal und dual nichtzulässigen Punkten starten, was diese besonders attraktiv macht. Im Folgenden werden wir zunächst ein zulässiges Verfahren untersuchen, welches sich dem Modellalgorithmus 6.9 unterordnet. Anschließend werden wir dieses zu einem nichtzulässigen Verfahren verallgemeinern.

Da jeder Punkt des zentralen Pfades Lösung eines nichtlinearen Gleichungssystems ist, kann man im allgemeinen dem zentralen Pfad nur approximativ folgen. Man definiert sich daher eine geeignete (und bisweilen recht großzügig angelegte) Umgebung des zentralen Pfades, aus der Iterierte noch akzeptiert werden. Insbesondere wird so erzwungen, dass die Iterierten dem Rand des nichtnegativen Orthanten nicht „zu nahe“ zu kommen.

Eine in diesem Zusammenhang verwendete Umgebung des zentralen Pfades ist die 2-Norm-Umgebung

${\mathcal {U}}_{2}(\vartheta ):=\{(x,y,s){\big |}x\in Z_{P}^{o},\quad (y,s)\in Z_{D}^{o},\quad \|XSe-\mu e\|\leq \vartheta \mu \}$

mit einem $\vartheta \in (0,1)$ und

$\mu :={\frac {x^{T}s}{n}}.$

Ein typischer Wert für $\vartheta$ in der Praxis ist $\vartheta :=0.5$.

Im Fall $\vartheta =0$ stimmt die Umgebung ${\mathcal {U}}_{2}(\vartheta )$ offenbar mit dem zentralen Pfad überein. Aber auch für $\vartheta >0$ erlaubt die Wahl von ${\mathcal {U}}_{2}(\vartheta )$ keine großen Abweichungen vom zentralen Pfad bzw. nur relativ kleine Schrittweiten $\alpha _{k}>0$. Aus praktischer Sicht sind daher Verfahren, welche eine 2-Norm-Umgebung verwenden, weniger empfehlenswert. Sie besitzen aber häufig schöne theoretische Konvergenzeigenschaften (s. [Wri97]).

Eine gebräuchliche und hier verwendete Umgebung des zentralen Pfades ist die einseitige $\infty$-Norm-Umgebung

(6.36) ${\mathcal {U}}_{-\infty }(\gamma ):=\{(x,y,s){\big |}x\in Z_{P}^{o},\quad (y,s)\in Z_{D}^{o},\quad x_{i}s_{i}\geq \gamma \mu \quad (i=1,\ldots ,n)\}$

mit einem $\gamma \in (0,1)$. Ihr Name resultiert aus den Implikationen

$\|XSe-\mu e\|_{\infty }\leq (1-\gamma )\mu \Leftrightarrow \pm (x_{i}s_{i}-\mu )\leq (1-\gamma )\mu \quad (i=1,\ldots ,n)\Rightarrow x_{i}s_{i}\geq \gamma \mu \quad (i=1,\ldots ,n).$

Mit diesen Implikationen kann man auch die Beziehung ${\mathcal {U}}_{-\infty }(\gamma )\supseteq {\mathcal {U}}_{2}(\vartheta )$ für $\gamma \leq 1-\vartheta$ schließen. Insbesondere stimmt ${\mathcal {U}}_{-\infty }(\gamma )$ für $\gamma =0$ mit dem Inneren $Z_{P}^{o}\times Z_{D}^{o}$ des primal-dual zulässigen Bereichs überein. Ein typischer, in der Praxis verwendeter Wert für $\gamma$ ist $\gamma :=10^{-3}$. Allgemein kann man sagen, dass Verfahren, die eine $\infty$-Norm-Umgebung des zentralen Pfades verwenden, meist größere Schrittweiten ermöglichen und daher in der Praxis schneller konvergieren als Verfahren mit einer 2-Norm-Umgebung. Jedoch sind ihre theoretisch nachgewiesenen Konvergenzeigenschaften häufig etwas schlechter als die letzterer.

Das folgende zulässige Verfahren erzeugt Iterierte in ${\mathcal {U}}_{-\infty }(\gamma )$ für ein vorgegebenes $\gamma \in (0,1)$ und ist aufgrund der Schrittweitenwahl ein sog. Long-Step-Verfahren. Wir verweisen in diesem Zusammenhang auf die in (6.32) und (6.33) eingeführten Notationen und auf Bemerkung 6.10.

Algorithmus 6.13

(0) Wähle $\gamma \in (0,1),\sigma _{\min }$ und $\sigma _{\max }$ mit $0<\sigma _{\min }<\sigma _{\max }<1$ und $(x^{0},y^{0},s^{0})\in {\mathcal {U}}_{-\infty }(\gamma )$. Setze $k:=0$.

(1) Wähle ein $\sigma _{k}\in [\sigma _{\min },\sigma _{\max }]$ und berechne $\mu _{k}:=\left((x^{k})^{T}s^{k}\right)/n$ sowie die Lösung $(\Delta x^{k},\Delta y^{k},\Delta s^{k})$ des linearen Gleichungssystems

(6.37) ${\begin{pmatrix}A&0&0\\0&A^{T}&I\\S_{k}&0&X_{k}\end{pmatrix}}{\begin{pmatrix}\Delta x\\\Delta y\\\Delta s\end{pmatrix}}={\begin{pmatrix}0\\0\\-X_{k}S_{k}e+\sigma _{k}\mu _{k}e\end{pmatrix}}.$

(2) Berechne

$\alpha _{k}:=\max \left\{\alpha \in [0,1]{\big |}\left(x^{k}(\alpha ),y^{k}(\alpha ),s^{k}(\alpha )\right)\in {\mathcal {U}}_{-\infty }(\gamma )\right\}.$

(3) Setze

$(x^{k+1},y^{k+1},s^{k+1}):=(x^{k},y^{k},s^{k})+\alpha _{k}(\Delta x^{k},\Delta y^{k},\Delta s^{k}).$

(4) Setze $k:=k+1$ und gehe nach (1).

Algorithmus 6.13 ist offenbar vom Typ des Modellalgorithmus 6.9, wobei aber durch die spezielle Wahl des Zentrierungsparameters $\sigma _{k}$ insbesondere die beiden Extremfälle $\sigma _{k}=0$ und $\sigma _{k}=1$ ausgeschlossen sind. Die Berechnung der Schrittweite $\alpha _{k}$ im Schritt (2) ist zwar prinzipiell möglich (der Leser möge dies verifizieren), sie erfordert jedoch einen gewissen numerischen Aufwand, so dass man meist eine kleinere, aber einfach zu bestimmende Schrittweite verwendet, für welche sich die Konvergenzeigenschaften des Verfahrens nicht ändern (s. hierzu Bemerkung 6.16).

Die Konvergenz von Algorithmus 6.13 wollen wir hier nicht im Detail beweisen (wir verweisen dafür auf [Wri97] und [Ree01]). Man kann insbesondere zeigen, dass die Schrittweiten $\alpha _{k}$ nach unten von Null weg beschränkt sind.

Lemma 6.14

Sei $(x^{k},y^{k},s^{k})\in {\mathcal {U}}_{-\infty }(\gamma )$. Dann hat man

$(x^{k}(\alpha ),y^{k}(\alpha ),s^{k}(\alpha ))\in {\mathcal {U}}_{-\infty }(\gamma )$

für alle $\alpha \in [0,{\bar {\alpha }}_{k}]$ mit

(6.38) ${\bar {\alpha }}_{k}:=2^{3/2}\gamma {\frac {\sigma _{k}}{n}}{\frac {1-\gamma }{1+\gamma }},$

wobei ${\bar {\alpha }}_{k}\leq 0.5$ für jedes $n\in \mathbb {N}$ ist.

Mit letzterem Ergebnis können wir nun die gewünschte Ungleichung vom Typ (6.35) beweisen.

Satz 6.15

Seien $(x^{k},y^{k},s^{k}),k\in \mathbb {N} _{0}$ die durch Algorithmus 6.13 erzeugten Iterierten. Dann gilt für eine von $k$ unabhängige Konstante $\delta >0$

$0\leq \mu _{k+1}\leq \left(1-{\frac {\delta }{n}}\right)\mu _{k}.$

Beweis.

Für $k\in \mathbb {N} _{0}$ schließt man aus Lemma 6.14 $\alpha _{k}\geq {\bar {\alpha }}_{k}$, so dass mit Lemma 6.11 (ii) folgt:

(6.39) $\mu _{k+1}=\mu _{k}(\alpha _{k})=(1-\alpha _{k}(1-\sigma _{k}))\mu _{k}\leq \left(1-{\frac {2^{3/2}}{n}}\gamma {\frac {1-\gamma }{1+\gamma }}\sigma _{k}(1-\sigma _{k})\right)\mu _{k}.$

Da die quadratische Funktion $\sigma \to \sigma (1-\sigma )$ konkav ist, nimmt sie ihr Minimum auf dem Intervall $[\sigma _{\min },\sigma _{\max }]$ in einem der Randpunkte an. Daher gilt für alle $\sigma _{k}\in [\sigma _{\min },\sigma _{\max }]$

$\sigma _{k}(1-\sigma _{k})\geq \min\{\sigma _{\min(}1-\sigma _{\min }),\sigma _{\max(}1-\sigma _{\max })\}.$

Mit

$\delta :=2^{3/2}\gamma {\frac {1-\gamma }{1+\gamma }}\min\{\sigma _{\min(}1-\sigma _{\min }),\sigma _{\max(}1-\sigma _{\max })\}$

folgt demnach die Behauptung.

q.e.d.

Bemerkung 6.16

Wählt man ${\bar {\alpha }}_{k}$ aus (6.38) anstelle von $\alpha _{k}$ als Schrittweite in Algorithmus 6.13, so gilt die Aussage von Satz 6.15 für das selbe $\delta$, wie man leicht aus dem Beweis des Satzes ersieht.

Verbindung der Sätze 6.15 und 6.12 liefert das entscheidende Konvergenzresultat für Algorithmus 6.13:

Satz 6.17

Seien $(x^{k},y^{k},s^{k}),k\in \mathbb {N} _{0}$ die durch Algorithmus 6.13 erzeugten Iterierten, wobei der Startvektor $(x^{0},y^{0},s^{0})\in {\mathcal {U}}_{-\infty }(\gamma )$ der Bedingung

(6.40) $0\leq \mu _{0}={\frac {(x^{0})^{T}s^{0}}{n}}\leq {\frac {1}{\varepsilon ^{\kappa }}}$

für ein $\kappa >0$ genüge. Dann hat man

$0\leq \mu _{k}\leq \varepsilon ,\quad k\geq K$

für

$K:=\left\lceil (1+\kappa ){\frac {n}{\delta }}|\ln(\varepsilon )|\right\rceil .$

Aufgrund von Bemerkung 6.16 erhält man das Konvergenzresultat von Satz 6.17 auch dann, wenn man in Algorithmus 2 in jeder Iteration anstelle von $\alpha _{k}$ die leicht berechenbare Konstante ${\bar {\alpha }}_{k}$ als Schrittweite verwendet. Weiter kann man zeigen (vgl. Bemerkung 6.10), dass die Folge $\{(x^{k},y^{k},s^{k})\}$ einen Häufungspunkt $(x^{*},y^{*},s^{*})$ besitzt und dass für jeden solchen Häufungspunkt die Vektoren $x^{*}$ und $(y^{*},s^{*})$ ein strikt komplementäres Lösungspaar von $(P)$ und $(D)$ bilden (s. [Wri97], S. 100 ff., und beachte, dass man $(y^{k},s^{k})\in Z_{D}^{o}$ hat und daher mit $\{s^{k}\}$ auch $\{y^{k}\}$ einen Häufungspunkt besitzt).

Für weitere Verfahren vom Typ des Modellalgorithmus 6.9 verweisen wir auf [Wri97].

6.6 Ein nichtzulässiges Verfahren

Wir beschreiben nun abschließend eine nichtzulässige Pfadverfolgungsmethode (infeasible path following method), welche man als eine Modifikation des Algorithmus 6.13 auffassen kann. Dieses Verfahren lässt also primale und duale Startpunkte zu, die nicht im Inneren der jeweiligen zulässigen Gebieten liegen müssen. Solche Methoden gelten heute als die effizientesten und besten Innere-Punkte-Verfahren.

Wir definieren dazu die Residuen

$r_{b}:=Ax-b,\quad r_{c}:=A^{T}y+s-c,$

wobei wir aus praktischen Gründen deren Abhängigkeit von $x$ bzw. $(y,s)$ nicht kennzeichnen, und wir bezeichnen weiter mit $r_{b}^{k}$ und $r_{c}^{k}$ die Residuen in Iterierten $x^{k}$ und $(y^{k},s^{k})$. Für Algorithmus 6.13 galt $r_{b}^{k}=0$ und $r_{c}^{k}=0$ für alle $k$, weshalb wir auch von einem zulässigen Verfahren sprachen. Jedoch ist für jenes Verfahren die Bestimmung eines Startvektors $(x^{0},y^{0},s^{0})$, der in ${\mathcal {U}}_{-\infty }(\gamma )$ liegt und zusätzlich der Bedingung (6.40) genügt, nicht unproblematisch.

Bei dem in diesem Abschnitt vorgestellten Verfahren müssen der Startpunkt und die Iterierten $x^{k}$ und $(y^{k},s^{k})$ nicht mehr innere Punkte des primalen bzw. dualen zulässigen Bereichs sein. Folglich werden hier die Voraussetzungen (A2) und (A3) aus Abschnitt 6.2 nicht benötigt. Statt dessen setzen wir neben der Rangbedingung (A1) in diesem Abschnitt nur voraus:

(A4) Das Problem $(P)$ besitzt eine Lösung.

Nach Satz 4.17 impliziert (A4), dass auch das Problem $(D)$ und das System $(S)$ lösbar sind.

Um für nichtzulässige Punkte, also Punkte mit $r_{b}\neq 0$ und $r_{c}\neq 0$, eine Umgebung des zentralen Pfades zu definieren, müssen wir die im vorigen Unterabschnitt eingeführte Menge ${\mathcal {U}}_{-\infty }(\gamma )$ geeignet erweitern. Wir definieren daher

(6.41) ${\mathcal {U}}_{-\infty }(\gamma ,\beta ):=\left\{(x,y,s){\big |}\|(r_{b},r_{c})\|\leq {\frac {\|(r_{b}^{0},r_{c}^{0})\|}{\mu _{0}}}\beta \mu ,\quad x>0,\quad s>0,\quad x_{i}s_{i}\geq \gamma \mu \quad (i=1,\ldots ,n)\right\},$

wobei $\gamma \in (0,1)$ und $\beta \geq 1$ gegebene Konstanten sind und sich $\mu _{0}$ und $(r_{b}^{0},r_{c}^{0})$ aus dem gewählten Startpunkt $(x^{0},y^{0},s^{0})$ ergeben. Die Forderung $\beta \geq 1$ ist hierbei nötig, damit der Startvektor $(x^{0},y^{0},s^{0})$ die in ${\mathcal {U}}_{-\infty }(\gamma ,\beta )$ gegenüber ${\mathcal {U}}_{-\infty }(\gamma )$ zusätzlich auftretende Bedingung

$\|(r_{b},r_{c})\|\leq {\frac {\|(r_{b}^{0},r_{c}^{0})\|}{\mu _{0}}}\beta \mu$

erfüllt und demnach in ${\mathcal {U}}_{-\infty }(\gamma ,\beta )$ liegt. Mit dieser Bedingung wird gemessen, inwieweit die linearen Gleichungen $Ax=b$ und $A^{T}y+s=c$ in $x$ bzw. $(y,s)$ verletzt sind.

Die Definition (6.41) einer Umgebung des zentralen Pfades garantiert für jede Folge $\{(x^{k},y^{k},s^{k})\}$ von Iterierten aus ${\mathcal {U}}_{-\infty }(\gamma ,\beta )$, dass mit $\mu _{k}\to 0$ auch $r_{b}^{k}\to 0$ und $r_{c}^{k}\to 0$ für $k\to \infty$ folgt und damit, ähnlich wie für zulässige Verfahren (vgl. Bemerkung 6.10), jeder Häufungspunkt der Folge $\left\{(x^{k},y^{k},s^{k})\right\}$ das System $(S)$ erfüllt bzw. Lösungen der Probleme $(P)$ und $(D)$ liefert. Also ist für ein Verfahren, welches die Umgebung ${\mathcal {U}}_{-\infty }(\gamma ,\beta )$ verwendet, das Hauptziel, die Konvergenz $\mu _{k}\to 0$ für $k\to \infty$ zu zeigen. Als Abbruchbedingung für ein solches Verfahren kann man dann wieder die Abfrage $0\leq \mu _{k}\leq \varepsilon$ für ein vorgegebenes $\varepsilon >0$ verwenden. Das zugehörige Iteriertenpaar $x^{k}$ und $(y^{k},s^{k})$ wird jedoch im allgemeinen nicht $\varepsilon$-optimal für $(P)$ und $(D)$ im Sinne von Korollar 4.19 sein, da $x^{k}$ und $(y^{k},s^{k})$ die Gleichungssysteme $Ax=b$ und $A^{T}y+s=c$ nicht exakt erfüllen müssen.

Der Algorithmus, den wir nun betrachten wollen, lautet unter Verwendung der Definitionen (6.32) und (6.33) folgendermaßen:

Algorithmus 6.18

(0) Wähle $\beta \geq 1,\gamma \in (0,1),\sigma _{\min }$ und $\sigma _{\max }$ mit $0<\sigma _{\min }<\sigma _{\max }<0.5$ sowie $x^{0}>0$ und $(y^{0},s^{0})$ mit $s^{0}>0$, so dass $(r_{b}^{0},r_{c}^{0})\neq (0,0)$ und $x_{i}^{0}s_{i}^{0}\geq \gamma \mu _{0}$ $(i=1,\ldots ,n)$ gelten. Setze $k:=0$.

(1) Für ein $\sigma _{k}\in [\sigma _{\min },\sigma _{\max }]$ und $\mu _{k}:=\left((x^{k})^{T}s^{k}\right)/n$ bestimme die Lösung $(\Delta x^{k},\Delta y^{k},\Delta s^{k})$ des linearen Gleichungssystems

(6.42) ${\begin{pmatrix}A&0&0\\0&A^{T}&I\\S_{k}&0&X_{k}\end{pmatrix}}{\begin{pmatrix}\Delta x\\\Delta y\\\Delta s\end{pmatrix}}={\begin{pmatrix}-r_{b}^{k}\\-r_{c}^{k}\\-X_{k}S_{k}e+\sigma _{k}\mu _{k}e\end{pmatrix}}.$

(2) Berechne

(6.43) $\alpha _{k}:=\max \left\{\alpha \in [0,1]{\big |}\left(x^{k}(\alpha ),y^{k}(\alpha ),s^{k}(\alpha )\right)\in {\mathcal {U}}_{-\infty }(\gamma ,\beta ),{\frac {(x^{k}(\alpha ))^{T}s^{k}(\alpha )}{n}}\leq (1-0.01\alpha )\mu _{k}\right\}.$

(3) Setze

$(x^{k+1},y^{k+1},s^{k+1}):=(x^{k},y^{k},s^{k})+\alpha _{k}(\Delta x^{k},\Delta y^{k},\Delta s^{k}).$

(4) Setze $k:=k+1$ und gehe nach (1).

Algorithmus 6.18 ähnelt offenbar dem Algorithmus 6.13 weitgehend. Die Voraussetzung $(r_{b}^{0},r_{c}^{0})\neq (0,0)$ ist im Folgenden formal notwendig. Wäre sie nicht erfüllt, d. h. wäre $(r_{b}^{0},r_{c}^{0})=(0,0)$, dann wären die Umgebungen ${\mathcal {U}}_{-\infty }(\gamma )$ und ${\mathcal {U}}_{-\infty }(\gamma ,\beta )$ in (6.36) bzw. (6.41) identisch und wäre Algorithmus 6.18 aufgrund der Schrittweitenwahl nur eine Variante von Algorithmus 6.13.

Die Bedingung $x_{i}^{0}s_{i}^{0}\geq \gamma \mu _{0}$ $(i=1,\ldots ,n)$ an den Startpunkt ist für genügend kleines $\gamma \in (0,1)$ immer erfüllbar und stellt sicher, dass $(x^{0},y^{0},s^{0})\in {\mathcal {U}}_{-\infty }(\gamma ,\beta )$ gilt. Da die Residuen $r_{c}^{k}$ und $r_{b}^{k}$, anders als bei dem Modellalgorithmus 6.9 und Algorithmus 6.13, nicht mehr notwendig gleich Null sind, ergibt sich hier das lineare Gleichungssystem (6.42) aus der Anwendung des Newton-Verfahrens auf das ursprüngliche System (6.23).

Aufgrund der Wahl des Startpunktes $(x^{0},y^{0},s^{0})$ (es ist $\beta \geq 1$) und der Wahl der Schrittweite in Schritt (2) des Verfahrens liegen alle durch Algorithmus 6.18 erzeugten Iterierten $(x^{k},y^{k},s^{k})$ in der Umgebung ${\mathcal {U}}_{-\infty }(\gamma ,\beta )$ des zentralen Pfades. Weiter impliziert die zweite Bedingung in (6.43) an die Wahl der Schrittweite $\alpha _{k}$ eine Abnahme des Dualitätsmaßes $\mu _{k}$.

Für Algorithmus 6.18 hat man nun das folgende Konvergenzresultat (s. [Wri97], [Ree01]).

Satz 6.19

Sei $\{(x^{k},y^{k},s^{k})\}$ eine durch Algorithmus 6.18 erzeugte Folge. Dann gilt für eine von $k$ unabhängige Konstante ${\bar {\alpha }}\in (0,1)$:

(i) $0\leq \mu _{k+1}\leq (1-0.01{\bar {\alpha }})\mu _{k},\quad k\in \mathbb {N} _{0},$

(ii) $\left\|(r_{b}^{k},r_{c}^{k})\right\|\leq {\frac {\|(r_{b}^{0},r_{c}^{0})\|}{\mu _{0}}}\beta \mu _{k},\quad k\in \mathbb {N} _{0}.$

Bemerkung 6.20

Bestimmt man die Schrittweite $\alpha _{k}$ von Algorithmus 6.18 nicht als den maximalen Wert aus der Menge $[0,1]$, sondern für ein vorgegebenes $\rho \in (0,1)$ als den maximalen Wert aus der Menge $\{1,\rho ,\rho ^{2},\ldots \}$, so dass die beiden Bedingungen in (6.43) erfüllt sind (letztere Zahl lässt sich leicht bestimmen), so gilt Satz 6.19 analog mit dem Wert $\rho {\bar {\alpha }}$ anstelle von ${\bar {\alpha }}$ (vgl. [Ree01]).

Offenbar impliziert Aussage (i) von Satz 6.19 die Konvergenz $\lim _{k\to \infty }\mu _{k}=0$. Aussage (ii) garantiert damit die Konvergenzen $\lim _{k\to \infty }r_{b}^{k}=\lim _{k\to \infty }r_{c}^{k}=0$ für die Residuen.

Darüber hinaus zeigt Satz 6.19, dass die Folge $\{\mu _{k}\}$ $q$-linear und die Folgen $\{r_{b}^{k}\}$ und $\{r_{c}^{k}\}$ $r$-linear gegen Null bzw. den Nullvektor konvergieren. Dabei heißt eine Nullfolge positiver reeller Zahlen $\{\alpha _{k}\}$ $q$-linear konvergent, wenn eine Konstante $c\in (0,1)$ existiert, so dass $\alpha _{k+1}\leq c\alpha _{k}$ für alle $k$ gilt und $r$-linear konvergent, wenn sie durch eine $q$-linear konvergente Folge $\{\beta _{k}\}$ majorisiert wird, d. h. wenn man $\alpha _{k}\leq \beta _{k}$ für alle $k$ hat.

Für Algorithmus 6.18 kann man ohne großen Mehraufwand sogar polynomiale Komplexität beweisen, wobei man allerdings die Wahl des Startvektors $(x^{0},y^{0},s^{0})$ einzuschränken hat (vgl. Satz 6.2 bei [Wri97]). Schließlich kann man unter den zusätzlichen Voraussetzungen $Z_{P}^{o}\neq \emptyset$ und $Z_{D}^{o}\neq \emptyset$ auch zeigen, dass die Folge $\{(x^{k},y^{k},s^{k})\}$ einen Häufungspunkt besitzt und dass jeder solche Häufungspunkt ein strikt komplementäres Lösungspaar von $(P)$ und $(D)$ liefert (vgl. [Wri97], S. 121 ff.).