%% LyX 1.3 created this file. For more info, see http://www.lyx.org/. %% Do not edit unless you really know what you are doing. \documentclass[oneside,german]{book} \usepackage[T1]{fontenc} \usepackage[latin1]{inputenc} \usepackage{a4wide} \usepackage{fancyhdr} \pagestyle{fancy} \setcounter{secnumdepth}{3} \setcounter{tocdepth}{3} \setlength\parskip{\medskipamount} \setlength\parindent{0pt} \usepackage{makeidx} \makeindex \usepackage{amssymb} \makeatletter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands. \usepackage{amssymb} \usepackage{amsmath} \usepackage{amsthm} \theoremstyle{plain} \newtheorem{thm}{Satz}[chapter] \numberwithin{equation}{chapter} \numberwithin{figure}{chapter} \theoremstyle{definition} \newtheorem{remark}[thm]{Bemerkung} \theoremstyle{definition} \newtheorem*{example*}{Beispiel} \theoremstyle{definition} \newtheorem{case}{Fall} %%Numbering of Cases not keyed to sections \theoremstyle{plain} \newtheorem{lem}[thm]{Lemma} %%Delete [thm] to re-start numbering \usepackage{algorithmic} \renewcommand{\algorithmicrequire}{\textbf{Input:}} \renewcommand{\algorithmicensure}{\textbf{Output:}} \newcommand{\COMPUTE}{\item[\textbf{Compute:}]} \theoremstyle{definition} \newtheorem{algorithm}[thm]{Algorithmus} \theoremstyle{definition} \newtheorem*{notation*}{Notation} \theoremstyle{definition} \newtheorem{conclusion}[thm]{Folgerung} \theoremstyle{plain} \newtheorem{cor}[thm]{Korollar} %%Delete [thm] to re-start numbering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands. \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\rank}{rank} \fancyhead[LO]{\slshape \leftmark} \fancyhead[RE]{\slshape \rightmark} \fancyhead[RO,LE]{\thepage} \fancyfoot[C]{} \usepackage{babel} \makeatother \begin{document} \title{Moderne Methoden der nichtlinearen endlichdimensionalen Optimierung} \author{Mitschrift% \footnote{Gefundene Fehler bitte an vigerske@mathematik.hu-berlin.de.% } zur Vorlesung von Prof. Griewank} \date{Sommersemester 2004} \maketitle \tableofcontents{} \chapter*{Allgemeine Voraussetzung\markboth{}{0. ALLGEMEINE VORAUSSETZUNG}} \begin{eqnarray*} & \min & f\left(x\right)\\ & \textrm{s.d.} & h\left(x\right)=0\\ & & c\left(x\right)\leq0\end{eqnarray*} wobei $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, $h:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, $c:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}$ und $f,g,h\in C^{2}\left(V\right)$ oder $C^{1}\left(V\right)$ für $V\subseteq\mathbb{R}^{n}$. \chapter{Freie Optimierungsprobleme } $m=p=0$ \section{\label{sub:AbstiegsminimierungStrahlsuche}Abstiegsminimierung mit Strahlsuche} \begin{eqnarray*} x_{k+1} & = & x_{k}+\alpha_{k}s_{k}\\ \textrm{mit $s_{k}$} & = & -B_{k}^{-1}g_{k}\end{eqnarray*} wobei $g_{k}=\nabla f\left(x_{k}\right)$, $B_{k}^{T}=B_{k}\approx\nabla^{2}f\left(x_{k}\right)$.\index{Abstiegsminimierung} \noindent Abstiegswinkel $\varphi_{k}:=\angle\left(s_{k},-g_{k}\right)$ für $B_{k}\succ0$, d.h. pos. definit\[ \cos\left(s_{k},-g_{k}\right)=\frac{-g_{k}^{T}s_{k}}{\left\Vert g_{k}\right\Vert \left\Vert s_{k}\right\Vert }\geq\frac{1}{\sqrt{\kappa\left(B_{k}\right)}}\qquad\kappa\left(B_{k}\right)=\left\Vert B_{k}\right\Vert \left\Vert B_{k}^{-1}\right\Vert \] \noindent Globale Konvergenz im Sinne von\[ \limsup_{k\rightarrow\infty}\left\Vert \nabla f\left(x_{k}\right)\right\Vert =0\qquad\textrm{oder}\qquad\lim_{k\rightarrow\infty}f\left(x_{k}\right)=-\infty\] folgt aus Kombination von Zoutendijkbedingung\index{Zoutendijkbedingung}\[ D:=\sum_{k=0}^{\infty}\cos^{2}\varphi_{k}\geq\sum_{k=0}^{\infty}\frac{1}{\kappa\left(B_{k}\right)}\] (d.h. $\sup\left\{ \kappa\left(B_{k}\right)\right\} <\infty$ reicht) mit effektiver Strahlsuche\index{Strahlsuche}. \noindent Falls zudem \[ \sup\left\Vert B_{k}\right\Vert <\infty,\qquad\sup\left\Vert B_{k}^{-1}\right\Vert <\infty\] folgt sogar\[ \lim_{k\rightarrow\infty}x_{k}=x_{*},\] d.h. Konvergenz zu stationärem Punkt, wo Optimalitätsbedingung 1. Ordnung erfüllt ist. \noindent Falls zudem $\nabla^{2}f\left(x_{*}\right)\succ0$ muß $x_{*}$ lokaler Minimalpunkt sein und $B_{k}=\nabla^{2}f\left(x_{k}\right)$ führt zu Q-quadratischer Konvergenz, d.h. \[ \left\Vert x_{k+1}-x_{*}\right\Vert \leq c\left\Vert x_{k}-x_{*}\right\Vert ^{2}\qquad\textrm{für Konstante }c.\] Falls nur $B_{k}\approx\nabla^{2}f\left(x_{k}\right)$ mittels Sekantenaufdatierung ergibt sich immer noch Q-superlineare Konvergenz, d.h.\[ \lim_{k\rightarrow\infty}\frac{\left\Vert x_{k+1}-x_{*}\right\Vert }{\left\Vert x_{k}-x_{*}\right\Vert }=0\] \begin{remark} \noindent Unter geeigneten Voraussetzungen lassen sich quadratische Konvergenz für Newton und superlineare Konvergenz für quasi-Newton auch bei NLP-Lösern, z.B. SQP und IP (Interior Point), beweisen. \end{remark} \section{Abstiegsminimierung mit Trust Region-Verfahren} \textbf{Motivation:} Problem bei Strahlsuche: Was tun, wenn $\nabla^{2}f\left(x\right)$ negative Eigenwerte hat? \begin{enumerate} \item nur für Iterierte $x=x_{n}$ (auf Weg) \item auch am stationären Punkt $x_{*}=\lim x_{n}$, d.h. notwendige Optimalitätsbedingung 2. Ordnung ist verletzt. \end{enumerate} Antwort 1: Korrektur zu 1.: Ersetze $\nabla^{2}f\left(x_{n}\right)$ durch positiv definites $B_{k}$, z.B. $B_{k}=\nabla^{2}f\left(x_{n}\right)+\lambda_{k}I$ mit $\lambda_{k}>0$. Diese Modifikation reicht nicht aus, um 2. zu beseitigen. \begin{example*} ~\begin{eqnarray*} f\left(x,y\right) & = & \frac{1}{2}x^{2}-\frac{1}{2}y^{2}\\ g\left(x,y\right)=\nabla f\left(x,y\right) & = & \left(x,-y\right)\\ \nabla^{2}f\left(x,y\right) & = & \left(\begin{array}{cc} 1 & 0\\ 0 & -1\end{array}\right)\quad\textrm{ist überall indefinit}\end{eqnarray*} Vom Ausganspunkt $\left(x_{0},y_{0}\right)$ mit $x_{0}\neq0$, $y_{0}=0$ ergibt die quasi-Newton Iteration mit \[ B_{k}=\left(\begin{array}{cc} 1+\lambda_{k} & 0\\ 0 & -1+\lambda_{k}\end{array}\right)\] für $\lambda_{k}\in\left[2,100\right]$ für alle $k\geq0$ \begin{eqnarray*} y_{k} & = & 0\\ x_{k+1} & = & x_{k}-x_{k}\frac{1}{1+\lambda_{k}}=\frac{\lambda_{k}}{1+\lambda_{k}}x_{k}\end{eqnarray*} $\leadsto$ $\frac{x_{k+1}}{x_{k}}\in\left[\frac{2}{3},\frac{100}{101}\right]$ $\leadsto$ $\left(x_{*},y_{*}\right)=\lim_{k\rightarrow\infty}\left(x_{k},y_{k}\right)=\left(0,0\right)$. Also konvergiert die Folge $\left(x_{n},y_{n}\right)$ gegen einen stationären Punkt, der nicht minimal ist. \end{example*} Antwort zu (1) und (2) ist das \emph{Trust Region Modell}\index{Trust Region!Modell}. \textbf{Idee:} Wähle \emph{Trust Region Radius}\index{Trust Region!Radius} $\Delta_{k}>0$, so daß\[ f\left(x_{k}+s\right)\approx\Phi_{k}\left(s\right):=f\left(x\right)+g\left(x\right)^{T}s+\frac{1}{2}s^{T}B_{k}s\qquad\textrm{für }\left\Vert s\right\Vert _{2}\leq\Delta_{k}\] \begin{enumerate} \item \emph{Bestimme} $s_{k}$ als Lösung von\[ \min\Phi_{k}\left(s\right)\quad\textrm{s.d. }\left\Vert s\right\Vert \leq\Delta_{k}\] \item \emph{Vergleiche} tatsächliche Funktionsdifferenz (actual reduction)\index{reduction!actual}\[ \textrm{ared}_{k}=-\left[f\left(x_{k}+s_{k}\right)-f\left(x_{k}\right)\right]=f\left(x_{k}\right)-f\left(x_{k}+s_{k}\right)\] mit vorhergesagter Reduktion (predicted reduction)\index{reduction!predicted}\[ \textrm{pred}_{k}=-\left[\Phi_{k}\left(s_{k}\right)-f\left(x_{k}\right)\right]=f\left(x_{k}\right)-\Phi_{k}\left(s_{k}\right)\geq0\] durch Berechnung des Quotienten\[ r_{k}=\frac{\textrm{ared}_{k}}{\textrm{pred}_{k}}.\] \item \emph{Akzeptiere} $x_{k+1}=x_{k}+s_{k}$, falls $r_{k}>\eta$, sonst setze $x_{k+1}=x_{k}$ und $\Delta_{k+1}\leq c_{4}\Delta_{k}$ mit $c_{4}\in\left(0,1\right)$. \end{enumerate} \begin{remark} Im Unterschied zur Strahlsuche wird bei Rückweisung von $x_{k}+s_{k}$ nicht nur die Schrittlänge, sondern auch die Schrittrichtung neu berechnet. Unter der schwachen Voraussetzung $\sup_{k}\left\Vert B_{k}\right\Vert <\infty$, d.h. $B_{k}$'s dürfen indefinit und speziell singulär sein, folgt unter gleichen Voraussetzungen wie in \ref{sub:AbstiegsminimierungStrahlsuche} für Strahlsuche {}``globale Konvergenz'' im Sinne von $\liminf\left\Vert g_{k}\right\Vert =0$. Unter zusätzlichen Voraussetzungen gilt sogar $x_{*}=\lim_{k\rightarrow\infty}x_{k}$, $g\left(x_{*}\right)=f\left(x_{*}\right)=0$, $\nabla^{2}f\left(x_{*}\right)\succeq0$. \end{remark} \subsection{Lösung des Trust Region Modells (exakt)} Nichtkonstruktive Existenzaussage: Die quadratische Funktion\[ \Phi\left(s\right)\equiv f+g^{T}s+\frac{1}{2}s^{T}Bs\] hat auf der kompakten Kugel \[ \left\{ s\in\mathbb{R}^{n}\left|c\left(s\right)\leq0\right.\right\} \qquad\textrm{mit }c\left(s\right)=\frac{1}{2}\left\Vert s\right\Vert ^{2}-\frac{1}{2}\Delta^{2}\] mindestens einen globalen Minimalpunkt $s_{*}\in\argmin\left\{ \Phi\left(s\right)\left|c\left(s\right)\leq0\right.\right\} $. \textbf{Fragen:} \begin{itemize} \item Gibt es weitere lokale Minima? \item Wie kann man globales Minimum eindeutig berechnen? \end{itemize} Zur Beantwortung betrachte (theoretische) Variablentransformation\[ s=Qz\qquad g=Qh\] wobei $B=Q\Lambda Q^{T}$, $QQ^{T}=I$, $\Lambda$ diagonal, so daß\begin{eqnarray*} \varphi\left(z\right):=\Phi\left(Qz\right) & = & f+g^{T}Qz+\frac{1}{2}z^{T}Q^{T}Q\Lambda Q^{T}Qz\\ & = & f+h^{T}z+\frac{1}{2}z^{T}\Lambda z\\ & = & f+\sum_{i=1}^{n}\left(h_{i}z_{i}+\frac{1}{2}\lambda_{i}z_{i}^{2}\right)\\ \textrm{und }c\left(Qz\right) & = & \frac{1}{2}\left\Vert Qz\right\Vert ^{2}-\frac{1}{2}\Delta^{2}=c\left(z\right)\end{eqnarray*} Mit anderen Worten: O.B.d.A. kann $B$ als diagonal angenommen werden. Untersuchung des transformierten Problemes\[ \min\varphi\left(z\right)\quad\textrm{s.d. }c\left(z\right)\leq0\] KKT-Bedingungen für stationären Punkt $z\neq0$ ($\leadsto$ $\nabla c\left(z\right)=z\neq0$, d.h. LICQ erfüllt):\begin{eqnarray*} 0 & = & \nabla\varphi\left(z\right)+\lambda\nabla c\left(z\right)=h+\Lambda z+\lambda z=h+\left(\Lambda+\lambda I\right)z\\ 0 & \leq & \lambda\\ c\left(z\right) & \leq & 0\\ 0 & = & \lambda c\left(z\right)=\lambda\left(\frac{1}{2}\left\Vert z\right\Vert ^{2}-\frac{1}{2}\Delta^{2}\right)\end{eqnarray*} \begin{case} Es existiert ein lokaler Minimalpunkt $z$ im Inneren, d.h. $c\left(z\right)<0$. \noindent $\leadsto$ $\lambda=0$ (wegen Komplementarität)\begin{eqnarray*} \nabla\varphi\left(z\right)=h+\Lambda z & = & 0\qquad\left(\textrm{notwendige Optimalität 1. Ordnung}\right)\\ \nabla^{2}\varphi\left(z\right)=\Lambda & \succeq & 0\qquad\left(\textrm{notwendige Optimalität 2. Ordnung}\right)\end{eqnarray*} Da aus der positiven Semidefinitheit von $\Lambda$ die Konvexität von $\varphi\left(z\right)$ folgt, sind diese notwendigen Bedingungen sogar hinreichend für globale Optimalität von $z$ im $\mathbb{R}^{n}$. \noindent Also gibt es in der Trust Region keine nichtglobalen lokalen Minima. \noindent Letzteres gilt sowohl bezüglich $z$ wie auch für die ursprünglichen Variablen $s$. \end{case} \begin{case} Es existiert ein lokales Minima $z$ am Rand, d.h. $c\left(z\right)=0$. \noindent Notwendige Bedingungen 1. Ordnung:\[ h+\left(\Lambda+\lambda I\right)z=0\qquad\frac{1}{2}\left\Vert z\right\Vert ^{2}-\frac{1}{2}\Delta^{2}=0\] \noindent Notwendige Bedingung 2. Ordnung:\[ \nabla c\left(z\right)^{T}\zeta=z^{T}\zeta=0\quad\Rightarrow\quad\zeta^{T}\left(\nabla_{z}^{2}L\left(z,\lambda\right)\right)\zeta=\zeta^{T}\left(\Lambda+\lambda I\right)\zeta\geq0\] wobei $L\left(z,\lambda\right)=f+h^{T}z+\frac{1}{2}z^{T}\Lambda z+\lambda\frac{1}{2}\left(\left\Vert z\right\Vert ^{2}-\Delta^{2}\right)$ \noindent Hinreichend für lokale Optimalität von $z$ ist die strikte Ungleichung \[ \zeta^{T}\left(\Lambda+\lambda I\right)\zeta=\sum_{i=1}^{n}\left(\lambda_{i}+\lambda\right)\zeta_{i}^{2}>0\] für alle $\zeta\neq0$ mit $z^{T}\zeta=0$. \noindent Selbst das reicht nicht für globale Optimalität aus, wie an folgendem Beispiel ersichtlich. \begin{example*} \noindent $n=1$, $\varphi\left(z\right)=-\frac{1}{2}\left(z-\frac{1}{2}\right)^{2}$, $\Delta=1$ \noindent Hinreichende Optimalitätsbedingung ist sowohl für $z=-1$ wie auch $z=1$ erfüllt, da $z^{T}\zeta=0$ nur für $\zeta=0$ gilt. \end{example*} \noindent Allgemein $\left(n\geq1\right)$ existieren maximal zwei lokal minimale Werte. Das globale Minimum erfüllt sogar $\Lambda+\lambda I\succeq0$ und läßt sich wie folgt charakterisieren. \noindent \textbf{Zusatzannahme}: $h_{i}\neq0$ für $i=1,\ldots,n$ \noindent $\leadsto$ $z=-\left(\Lambda+\lambda I\right)^{-1}h$ \noindent $\leadsto$ $\lambda$ ist Nullstelle der Funktion \begin{eqnarray*} \psi\left(\lambda\right)=c\left(z\left(\lambda\right)\right) & = & \frac{1}{2}\left\Vert \left(\Lambda+\lambda I\right)^{-1}h\right\Vert ^{2}-\frac{1}{2}\Delta^{2}\\ & = & \frac{1}{2}\sum_{i=1}^{n}\frac{h_{i}^{2}}{\left(\lambda_{i}+\lambda\right)^{2}}-\frac{1}{2}\Delta^{2}\end{eqnarray*} \noindent Sei $\lambda_{1}<\lambda_{2}<\ldots<\lambda_{n}$. $\psi\left(\lambda\right)$ hat bis zu $2n$ Nullstellen von denen genau eine $\lambda=\lambda_{*}$ größer als $-\lambda_{1}$ sein muß, da $h_{1}\neq0$ vorausgesetzt ist. \noindent Für dieses $\lambda_{*}$ ist $\left(\Lambda+\lambda_{*}I\right)=\textrm{diag}\left(\lambda_{i}+\lambda\right)_{i=1,\ldots,n}$ sogar positiv definit, die hinreichenden Optimalitätsbedingungen 2. Ordnung sind also an $z_{*}=-\left(\Lambda+\lambda_{*}I\right)^{-1}h$ erfüllt. $z_{*}$ ist sogar globales Minimum, da für beliebiges $\zeta$ mit $c\left(z_{*}+\zeta\right)\leq0$ gilt \begin{eqnarray*} 0\geq c\left(z_{*}+\zeta\right) & = & \frac{1}{2}\left\Vert z_{*}+\zeta\right\Vert ^{2}-\Delta^{2}\\ & = & \frac{1}{2}\left\Vert z_{*}\right\Vert ^{2}+z_{*}^{T}\zeta+\frac{1}{2}\left\Vert \zeta\right\Vert ^{2}-\frac{1}{2}\Delta^{2}\\ & = & z_{*}^{T}\zeta+\frac{1}{2}\left\Vert \zeta\right\Vert ^{2}\\ \leadsto\,-z_{*}^{T}\zeta & \geq & \frac{1}{2}\left\Vert \zeta\right\Vert ^{2}\\ \leadsto\,\varphi\left(z_{*}+\zeta\right)-\varphi\left(z_{*}\right) & = & h^{T}\zeta+\frac{1}{2}\left(z_{*}+\zeta\right)^{T}\Lambda\left(z_{*}+\zeta\right)-\frac{1}{2}z_{*}^{T}\Lambda z_{*}\\ & = & h^{T}\zeta+z_{*}^{T}\Lambda\zeta+\frac{1}{2}\zeta^{T}\Lambda\zeta\\ & = & \zeta^{T}\left(h+\Lambda z_{*}\right)+\frac{1}{2}\zeta^{T}\Lambda\zeta\\ & = & \zeta^{T}\left(-\lambda I\right)z_{*}+\frac{1}{2}\zeta^{T}\Lambda\zeta\\ & \geq & \frac{1}{2}\lambda\left\Vert \zeta\right\Vert ^{2}+\frac{1}{2}\zeta^{T}\Lambda\zeta\\ & = & \frac{1}{2}\zeta^{T}\left(\Lambda+\lambda I\right)\zeta\\ & > & 0\end{eqnarray*} \end{case} \begin{thm} $s_{*}$ ist globales Minimum des Problems\[ \min\Phi\left(s\right)=f+g^{T}s+\frac{1}{2}s^{T}Bs\qquad\textrm{s.d. }\left\Vert s\right\Vert \leq\Delta\] gdw. für ein $\lambda_{*}\geq0$\begin{eqnarray*} Bs_{*}+g+\lambda_{*}s_{*} & = & 0,\\ \left(\Delta-\left\Vert s_{*}\right\Vert \right)\geq0 & = & \lambda_{*}\left(\Delta-\left\Vert s_{*}\right\Vert \right)\\ B+\lambda_{*}I & \succeq & 0\end{eqnarray*} \end{thm} \begin{proof} Falls $h=Q^{T}g$ keine verschwindenden Elemente hat, ergibt sich die Aussage aus vorgehenden Betrachtungen für den Fall 1 und 2. Anderfalls betrachte das gestörte Problem mit\[ h_{\varepsilon}=Q^{T}g+\varepsilon\mathbf{1}=Q^{T}\left(g+\varepsilon Q\mathbf{1}\right)\quad\textrm{für }\varepsilon=0,\mathbf{1}=\left(1,\ldots,1\right)^{T}\] Dann benutze Stetigkeitsargument, d.h. für alle hinreichend kleinen positiven $\varepsilon>0$ hat $h_{\varepsilon}$ keine verschwindenden Elemente, so daß es $s_{\varepsilon}$ und $\lambda_{\varepsilon}$ gibt mit\[ \left\Vert s_{\varepsilon}\right\Vert \leq\Delta\quad\textrm{und}\quad Bs_{*}+g+\lambda_{\varepsilon}s_{\varepsilon}=0.\] Da $\left\{ s_{\varepsilon}\right\} $ beschränkt ist, gibt es mindestens einen Häufungspunkt $s_{*}$. Falls es keinen nichtverschwindenden Häufungspunkt gibt, gilt $0=\lim_{\varepsilon\rightarrow0}s_{\varepsilon}$, so daß für hinreichend kleine $\varepsilon$ alle Minimalpunkte $s_{\varepsilon}$ unbeschränkt sind und deshalb $0=Bs_{\varepsilon}+g_{\varepsilon}$. Dann folgt $\lim_{\varepsilon\rightarrow0}g_{\varepsilon}=\lim_{\varepsilon\rightarrow0}g+\varepsilon Q\mathbf{1}=0$ $\Rightarrow$ $g=0$, $B\succeq0$. Also ist $s_{*}=0$ globales Minimum des ungestörten Problemes und der Satz gilt z.B. für $\lambda_{*}=0$. Falls ein nichtverschwindender Häufungspunkt $s_{*}=\lim_{j\rightarrow\infty}s_{\varepsilon_{j}}$ mit $\lim_{j\rightarrow\infty}\varepsilon_{j}=0$ existiert, so folgt für die entsprechenden $\lambda_{\varepsilon_{j}}$, daß\[ \limsup_{j\rightarrow\infty}\left|\lambda_{\varepsilon_{j}}\right|\leq\limsup_{j\rightarrow\infty}\frac{\left\Vert Bs_{\varepsilon_{j}}+g+\varepsilon_{j}Q^{T}\mathbf{1}\right\Vert }{\left\Vert s_{\varepsilon_{j}}\right\Vert }\leq\left\Vert B\right\Vert +\frac{\left\Vert g\right\Vert }{\left\Vert s_{*}\right\Vert }.\] Also besitzen auch die $\lambda_{\varepsilon_{j}}$ einen Häufungspunkt $\lambda_{*}$ der mit $s_{*}$ die obigen Bedingungen erfüllt. \end{proof} Falls $B\succ0$ und deshalb $s^{q}=-B^{-1}g$ wohl definiert ist, ergibt sich ein glatter Pfad\begin{eqnarray*} s\left(\lambda\right) & = & -\left(B+\lambda I\right)^{-1}g\\ \textrm{mit }s\left(0\right) & = & s^{q}=-B^{-1}g\\ \textrm{und }s\left(\lambda\right) & \approx & -\frac{1}{\lambda}g\qquad\textrm{für }\lambda\gg0\end{eqnarray*} Im normalen, nichtdegenerierten Fall löst $\lambda$ die Skalargleichung\begin{eqnarray*} 0=\varphi\left(\lambda\right) & = & \left\Vert \left(B+\lambda I\right)^{-1}g\right\Vert -\Delta\\ & = & \left[g^{T}\left(B-\lambda I\right)^{-2}g\right]^{\frac{1}{2}}-\Delta\\ & \approx & \frac{\left|g^{T}Qe_{1}\right|}{\lambda+\lambda_{1}}-\Delta\qquad\textrm{falls }\lambda+\lambda_{1}\gtrsim0,\textrm{ d.h. klein positiv}\end{eqnarray*} Für die Nullstellensuche ist die transformierte (annähernd lineare) Funktion\[ \psi\left(\lambda\right)=\frac{1}{\varphi\left(\lambda\right)+\Delta}-\frac{1}{\Delta}\approx\frac{\left(\lambda+\lambda_{1}\right)}{\left|g^{T}Qe_{1}\right|}-\frac{1}{\Delta}\] besser geeignet. Iterative Lösung mit Newton verlangt Ableitungen\[ \varphi'\left(\lambda\right)=\frac{d}{d\lambda}\varphi\left(\lambda\right)=\frac{\frac{1}{2}\left(-2\right)}{\varphi\left(\lambda\right)+\Delta}g^{T}\left(B+\lambda I\right)^{-3}g=-\frac{g^{T}\left(B+\lambda I\right)^{-3}g}{\varphi\left(\lambda\right)+\Delta}<0,\] denn wegen\[ 0=\frac{d}{d\lambda}\left[\left(B+\lambda I\right)^{-2}\left(B+\lambda I\right)\left(B+\lambda I\right)\right]=\frac{d}{d\lambda}\left[\left(B+\lambda I\right)^{-2}\right]\cdot\left(B+\lambda I\right)^{2}+2\left(B+\lambda I\right)^{-2}I\left(B+\lambda I\right)\] ist \[ \frac{d}{d\lambda}\left[\left(B+\lambda I\right)^{-2}\right]=-2\left(B+\lambda I\right)^{-3}.\] Entsprechend gilt\[ \psi'\left(\lambda\right)=-\frac{\varphi'\left(\lambda\right)}{\left[\varphi\left(\lambda\right)+\Delta\right]^{2}}=\frac{g^{T}\left(B+\lambda I\right)^{-3}g}{\left(\varphi\left(\lambda\right)+\Delta\right)^{3}}>0,\] falls $B+\lambda I\succ0$. Newton-Iteration (mit geeigneter Stabilisierung)\[ \lambda\leftarrow\lambda-\frac{\psi\left(\lambda\right)}{\psi'\left(\lambda\right)}=\lambda-\frac{\varphi\left(\lambda\right)\left(\varphi\left(\lambda\right)+\Delta\right)}{\Delta\cdot\varphi'\left(\lambda\right)}\] verlangt an Hauptarbeit die Lösung der linearen Systeme\[ \left(B+\lambda I\right)s=g\quad\textrm{und}\quad\left(B+\lambda I\right)t=s\] für eine Folge von Multiplikatoren $\lambda$. Für eine direkte Lösung bestehen folgende Möglichkeiten: \begin{enumerate} \item Die einmalige Eigenwertzerlegung $B=Q\Lambda Q^{T}$ zum Preis von $\approx4n^{3}$ Operationen ermöglicht die sukzessive Lösung des transformierten Problemes mit Aufwand $O\left(n\right)$ pro Iteration. \item Eine wiederholte Cholesky-Faktorisierung $B+\lambda I=LDL^{T}$ mit $D\succ0$ verlangt jedesmal einen Aufwand von $\approx\frac{1}{6}n^{3}$ pro Iteration. (Mor\'e Sorenson) \item \noindent Eine einmalige Reduktion von $B$ auf symmetrische Hessenberg-tridiagonale Form $B=QTQ^{T}$ mit tridiagonaler Matrix $T$ und $Q^{T}Q=I$ kostet $\frac{2}{3}n^{3}$ Operationen. (Ähnlichkeitstransformation auf Hessenbergform wird allgemein als erster Schritt bei der Berechnung von Eigenwerten eingesetzt ($QR$-Algorithmus). Die Berechnung von $Q$ erfolgt durch Householder Matrizen (elementare Reflektoren).) \noindent In transformierten Systemen gilt dann\begin{eqnarray*} \varphi\left(\lambda\right) & = & \left\Vert \tilde{s}\right\Vert -\Delta\quad\textrm{mit }\left(T+\lambda I\right)\tilde{s}=h=Q^{T}g\\ \varphi'\left(\lambda\right) & = & \frac{\tilde{s}^{T}\tilde{t}}{\varphi\left(\lambda\right)+\Delta}\quad\textrm{mit }\left(T+\lambda I\right)\tilde{t}=\tilde{s}\end{eqnarray*} $\tilde{s}$ und $\tilde{t}$ lassen sich mittels Cholesky-Verfahren mit Aufwand $O\left(n\right)$ berechnen. \end{enumerate} \begin{remark} Das Trust-Region-Problem wird heute nur bei Problemen mit mittlerer Größe, d.h. $\approx300$ Variablen, genau gelöst. Bei wirklich großen Problemen kommen folgende Vergröberungen zum Einsatz. \end{remark} \subsection{Vereinfachungen einschließlich Steihaug's CG-Variante} Zur Vermeidung jeglicher Gleichungslösung betrachte den \emph{Cauchy-Punkt}\index{Cauchy-Punkt}\[ s^{c}=-g\cdot\argmin_{0\leq\alpha\leq\Delta/\left\Vert g\right\Vert }\left(\Phi\left(-\alpha g\right)\right)=-g\cdot\argmin_{0\leq\alpha\leq\Delta/\left\Vert g\right\Vert }\left(-\alpha\left\Vert g\right\Vert ^{2}+\frac{1}{2}\alpha^{2}g^{T}Bg\right)\] Freier Minimalpunkt an Stelle $\alpha=\frac{\left\Vert g\right\Vert ^{2}}{g^{T}Bg}$. \begin{lem} Der Cauchy-Punkt $s=s^{c}$ verspricht in jedem Fall die Reduktion \[ \mathrm{pred}=f-\Phi\left(s\right)\geq\frac{1}{2}\left\Vert g\right\Vert \min\left(\Delta,\frac{\left\Vert g\right\Vert }{\left\Vert B\right\Vert }\right).\] \end{lem} Falls $g^{T}Bg\leq0$ folgt unmittelbar $s=-g\frac{\Delta}{\left\Vert g\right\Vert }$ und\[ f-\Phi\left(s\right)=-\left[\underbrace{-\left\Vert g\right\Vert ^{2}\frac{\Delta}{\left\Vert g\right\Vert }}_{g^{T}s}+\frac{1}{2}\underbrace{\frac{g^{T}Bg\Delta^{2}}{\left\Vert g\right\Vert ^{2}}}_{s^{T}Bs}\right]\geq\left\Vert g\right\Vert \Delta.\] Falls $g^{T}Bg>0$ und $\alpha=\frac{\left\Vert g\right\Vert ^{2}}{g^{T}Bg}\leq\Delta$ gilt, so ist\[ f-\Phi\left(s\right)=-\left[-\frac{\left\Vert g\right\Vert ^{4}}{g^{T}Bg}+\frac{1}{2}g^{T}Bg\frac{\left\Vert g\right\Vert ^{4}}{\left(g^{T}Bg\right)^{2}}\right]=\frac{1}{2}\frac{\left\Vert g\right\Vert ^{4}}{g^{T}Bg}\geq\frac{1}{2}\frac{\left\Vert g\right\Vert ^{2}}{\left\Vert B\right\Vert }=\frac{1}{2}\left\Vert g\right\Vert \left(\frac{\left\Vert g\right\Vert }{\left\Vert B\right\Vert }\right).\] Falls $g^{T}Bg>0$ und $\frac{\left\Vert g\right\Vert ^{2}}{g^{T}Bg}\geq\frac{\Delta}{\left\Vert g\right\Vert }$, d.h. $\frac{g^{T}Bg}{\left\Vert g\right\Vert ^{2}}\leq\frac{\left\Vert g\right\Vert }{\Delta}$, gilt, so ist $\alpha=\frac{\Delta}{\left\Vert g\right\Vert }$, $s=-\frac{g}{\left\Vert g\right\Vert }\Delta$ und\[ f-\Phi\left(s\right)=-\left[-\Delta\left\Vert g\right\Vert +\frac{1}{2}\frac{\Delta^{2}}{\left\Vert g\right\Vert ^{2}}g^{T}Bg\right]\geq\Delta\left\Vert g\right\Vert -\frac{1}{2}\frac{\Delta^{2}\left\Vert g\right\Vert }{\Delta}=\frac{1}{2}\Delta\left\Vert g\right\Vert .\] \textbf{Idee}:\index{Steihaug's CG-Variante} Falls der Cauchy $s^{c}$ noch im Inneren der Trust-Region liegt, entspricht er genau der ersten Iterierten $s_{1}$ des Conjugate-Gradient(CG)-Verfahrens, angewandt auf $Bs=-g$, gestartet mit $s_{0}=0$. Die Iteration kann fortgesetzt werden, bis zum ersten mal \begin{enumerate} \item eine Richtung $\zeta_{k}$ mit $\zeta_{k}B\zeta_{k}\leq0$ auftritt, oder \item ein normaler CG-Schritt $s_{k}-\frac{\nabla\Phi\left(s_{k}\right)\zeta_{k}}{\zeta_{k}B\zeta_{k}}$ die Trust-Region verlassen würde. \end{enumerate} In beiden Fällen wird $\alpha_{k}$ so gewählt, daß $\left\Vert s_{k}+\alpha_{k}\zeta_{k}\right\Vert =\Delta$ gilt. Falls weder 1. noch 2. auftreten, wird nach spätestens $n$ Schritten der Minimalpunkt $s_{*}=-B^{-1}g$ von $\Phi\left(s\right)$ erreicht. Letzeres folgt aus folgendem Lemma: \begin{lem} Das CG-Verfahren mit $s_{0}=0$ verläßt die Trust-Region höchstens genau einmal, da die Norm $\left\Vert s_{k}\right\Vert =\left\Vert s_{k}\right\Vert _{2}$ monoton wächst, solange $\zeta_{k}^{T}B\zeta_{k}>0$ ist. \end{lem} \begin{proof} Die Aussage folgt offenbar aus $s_{k}^{T}\zeta_{k}\geq0$ für alle $k=0,1,\ldots$. Wir zeigen dies durch Induktion über $k$: \[ s_{k+1}=s_{k}+\alpha_{k}\zeta_{k}\] Induktionsanfang $k=0$: $\zeta_{0}=-g$, $s_{0}=0$ $\Rightarrow$ $s_{0}^{T}\zeta_{0}=0\geq0$ Induktionsschluß: \[ s_{k+1}^{T}\zeta_{k+1}=s_{k+1}^{T}\left(-\nabla\Phi\left(s_{k+1}\right)+\underbrace{\beta_{k}}_{\geq0}\zeta_{k}\right),\] wobei $s_{k+1}$ eine Linearkombination der Schritt $\zeta_{j}$ für $j\leq k$ ist. Wegen der Expanding-Subspace-Eigenschaft ist nach dem CG-Theorem der aktuelle Gradient $\nabla\Phi\left(s_{k+1}\right)$ zu allen Suchrichtungen $\zeta_{j}$ mit $j\leq k$ orthogonal. Also gilt\begin{align*} s_{k+1}^{T}\zeta_{k+1} & =s_{k+1}^{T}\beta_{k}\zeta_{k}\\ & =\left(s_{k}+\alpha_{k}\zeta_{k}\right)^{T}\beta_{k}\zeta_{k}\\ & =\beta_{k}\underbrace{s_{k}^{T}\zeta_{k}}_{\geq0\textrm{ nach IV}}+\underbrace{\alpha_{k}\beta_{k}\left\Vert \zeta_{k}\right\Vert ^{2}}_{\geq0}\\ & \geq0.\qedhere\end{align*} \end{proof} \subsubsection*{Letzte Vereinfachung: Dog-Leg} Falls $B\succ0$, so betrachte die stückweise lineare Kurve, die den Cauchy-Punkt $s^{c}$ und das unbeschränkte Minimum $s=-B^{-1}g$ verbindet. Der Dog-Leg-Schritt\index{Dog-Leg} $s^{d}$ ist dann der Punkt auf dieser Kurve mit $\Phi\left(s^{d}\right)\leq\Phi\left(s^{c}\right)$. \subsection{Beweis der globalen Konvergenz} \begin{algorithm} Betrachte den Algorithmus (nach \cite{NoWr}): \begin{algorithmic}\STATE Initialisiere $\Delta_{0}>0$ und $\eta\in\left[0,1\right)$ sowie $x_{0}\in\mathbb{R}^{n}$. \FOR{$k=0,1,\ldots$} \STATE Wähle $s_{k}\in\mathbb{R}^{n}$, so daß\[ \textrm{pred}_{k}=f\left(x_{k}\right)-\Phi_{k}\left(s_{k}\right)\geq\frac{1}{2}\left\Vert g_{k}\right\Vert \cdot\min\left(\Delta_{k},\frac{\left\Vert g_{k}\right\Vert }{\left\Vert B_{k}\right\Vert }\right)\] \STATE Berechne \[ r_{k}=\frac{f\left(x_{k}\right)-f\left(x_{k}+s_{k}\right)}{\textrm{pred}_{k}}=\frac{\textrm{ared}_{k}}{\textrm{pred}_{k}}\] \IF{$r_{k}\geq\eta$}\STATE setze $x_{k+1}=x_{k}+s_{k}$ \ELSE \STATE $x_{k+1}=x_{k}$ \ENDIF \IF{ $r_{k}<\frac{1}{4}$}\STATE setze $\Delta_{k+1}=\frac{1}{4}\left\Vert s_{k}\right\Vert \leq\frac{1}{4}\Delta_{k}$ \ELSIF{$r_{k}>\frac{3}{4}$ und $\left\Vert s_{k}\right\Vert =\Delta_{k}$}\STATE $\Delta_{k+1}=2\Delta_{k}$ \ELSE \STATE $\Delta_{k+1}=\Delta_{k}$ \ENDIF \ENDFOR\end{algorithmic} \end{algorithm} \begin{thm} \label{thm:KonvergenzTrustRegion}Angenommen, die Menge $\left\{ x\in\mathbb{R}^{n}\left|f\left(x\right)\leq f\left(x_{0}\right)\right.\right\} $ ist für gegebenes $x_{0}\in\mathbb{R}^{n}$ beschränkt und $f$ ist auf dieser Niveaumenge Lipschitz-stetig differenzierbar. Falls zudem $\beta=\sup_{k>0}\left\Vert B_{k}\right\Vert <\infty$, so hat die erzeugte Iterationsfolge $x_{k}$ die Eigenschaft\[ \liminf_{k\rightarrow\infty}\left\Vert \nabla f\left(x_{k}\right)\right\Vert =0.\] Insbesondere konvergiert eine Teilfolge zu einem stationären Punkt $x_{*}\in\mathbb{R}^{n}$. \end{thm} \begin{proof} Entweder haben wir $r_{k}\geq\frac{1}{4}$ unendlich oft, oder $\Delta_{k}$ wird nur endlich oft nicht um einen Faktor $\frac{1}{4}$ reduziert, so daß $\lim_{k\rightarrow\infty}\Delta_{k}=0$ gilt. Im ersten Fall gilt unendlich oft \[ \textrm{ared}_{k}\equiv f\left(x_{k}\right)-f\left(x_{k+1}\right)\geq\frac{1}{4}\textrm{pred}_{k}=\frac{1}{8}\left\Vert g_{k}\right\Vert \cdot\min\left(\Delta_{k},\frac{\left\Vert g_{k}\right\Vert }{\beta}\right).\] Da die Summe aller akzeptierten $\textrm{ared}_{k}\geq0$ durch $f\left(x_{0}\right)-\min f\left(x\right)$ beschränkt ist, muß gelten\[ \liminf_{k\rightarrow\infty}\textrm{ared}_{k}=0.\] Daraus folgt\[ \liminf_{k\rightarrow\infty}\left\Vert g_{k}\right\Vert \cdot\min\left(\Delta_{k},\frac{\left\Vert g_{k}\right\Vert }{\beta}\right)=0\] und somit entweder\[ \liminf_{k\rightarrow\infty}\left\Vert g_{k}\right\Vert =0\] oder\[ \liminf_{k\rightarrow\infty}\Delta_{k}=0.\] Auszuschließen bleibt nur noch der zweite Fall, d.h. wir nehmen o.B.d.A. durch Weglassen von Zwischeniterationen an, daß\[ \lim_{k\rightarrow\infty}\Delta_{k}=0\quad\textrm{und}\quad r_{k}\leq\frac{1}{4}\textrm{ für alle }k.\] Für diese (Teil)folge gilt also\[ -g_{k}^{T}s_{k}-\frac{1}{2}L\left\Vert s_{k}\right\Vert ^{2}\leq\textrm{ared}_{k}\leq\frac{1}{4}\textrm{pred}_{k}\leq\frac{1}{4}\left(-g_{k}^{T}s_{k}+\underbrace{\frac{1}{2}\beta\left\Vert s_{k}\right\Vert ^{2}}_{\geq-\frac{1}{2}s_{k}^{T}Bs_{k}}\right)\] wobei $L$ die Lipschitzkonstante von $\nabla f\left(x\right)$sei. Damit folgt\begin{align*} \limsup_{k\rightarrow\infty}\frac{3}{4}\frac{\left(-g_{k}^{T}s_{k}\right)}{\Delta_{k}} & \leq\limsup_{k\rightarrow\infty}\frac{\left\Vert s_{k}\right\Vert ^{2}}{\Delta_{k}}\left(\frac{1}{8}\beta+\frac{1}{2}L\right)=0\\ \Rightarrow\,\limsup_{k\rightarrow\infty}\frac{\textrm{pred}_{k}}{\Delta_{k}} & \leq\limsup_{k\rightarrow\infty}\frac{-g_{k}^{T}s_{k}+\frac{1}{2}\beta\left\Vert s_{k}\right\Vert ^{2}}{\Delta_{k}}=0\\ \Rightarrow\,\limsup_{k\rightarrow\infty}\left\Vert g_{k}\right\Vert \min\left(1,\frac{\left\Vert g_{k}\right\Vert }{\Delta_{k}\beta}\right) & =0\\ \Rightarrow\,\limsup_{k\rightarrow\infty}\frac{\left\Vert g_{k}\right\Vert }{\Delta_{k}} & =0\\ \Rightarrow\,\liminf_{k\rightarrow\infty}\left\Vert g_{k}\right\Vert & =0.\qedhere\end{align*} \end{proof} \begin{thm} \label{thm:KonvergenzTrustRegion2}Unter den Voraussetzungen von Satz \ref{thm:KonvergenzTrustRegion} erzwingt $0<\eta<\frac{1}{4}$ sogar \[ \lim_{k\rightarrow\infty}g_{k}=0.\] Also ist jeder Häufungspunkt von $\left\{ x_{k}\right\} $ stationär. \end{thm} \begin{proof} Siehe \cite{NoWr}. \end{proof} Um auch die notwendigen Optimalitätsbedingungen 2. Ordnung an allen Häufungspunkten zu erzwingen, setze $B_{k}=\nabla^{2}f\left(x_{k}\right)$. Mit anderen Worten, wir betrachten die Trust-Region-Variante des Newton-Verfahrens. Dann folgt für die exakte Trust-Region-Lösung = Minimalpunkt $s_{k}$\begin{eqnarray*} -f\left(x_{k}\right)+\Phi_{k}\left(s_{k}\right) & = & \min_{s:\left\Vert s\right\Vert \leq\Delta}g_{k}^{T}s+\frac{1}{2}s^{T}\nabla^{2}f\left(x_{k}\right)s\\ & \leq & \min_{s:\left\Vert s\right\Vert \leq\Delta}\frac{1}{2}s^{T}\nabla^{2}f\left(x_{k}\right)s\\ & \leq & \min_{s:\left\Vert s\right\Vert \leq\Delta}\frac{1}{2}\left\Vert s\right\Vert ^{2}\lambda_{\min}\\ & \leq & \frac{1}{2}\Delta^{2}\lambda_{\min},\end{eqnarray*} wobei $\lambda_{\min}=\lambda_{\min}\left(\nabla^{2}f\left(x_{k}\right)\right)$ den kleinsten Eigenwert von $\nabla^{2}f\left(x_{k}\right)$ bezeichne. Also gilt \[ \textrm{pred}_{k}\geq-\frac{1}{2}\Delta^{2}\lambda_{\min}>0\quad\Leftrightarrow\quad\lambda_{\min}<0.\] \begin{thm} Unter den Voraussetzungen von Satz \ref{thm:KonvergenzTrustRegion2} impliziert \[ B_{k}=\nabla^{2}f\left(x_{k}\right)\quad\textrm{und}\quad s_{k}=\argmin_{s:\left\Vert s\right\Vert \leq\Delta}\Phi_{k}\left(s\right),\] daß \[ \nabla^{2}f\left(x_{*}\right)\succeq0\] für alle Häufungspunkte $x_{*}$ der Folge $\left\{ x_{k}\right\} $. \end{thm} \begin{proof} Siehe \cite{JaSt}. \end{proof} \begin{remark} Da $\left\Vert s\right\Vert =\left\Vert s\right\Vert _{2}$ sehr stark von der Skalierung der Variablen $e_{i}^{T}x$, d.h. der Einheiten in denen diese gemessen werden, abhängt, benutzt man oft Trust Regions der Form\[ \left\Vert s\right\Vert =\left\Vert Ds\right\Vert _{2},\] wobei $D=\textrm{diag}\left(d_{j}\right)_{j=1}^{n}$ mit $d_{j}>0$ oder\begin{eqnarray*} \left\Vert s\right\Vert & = & \left\Vert Ds\right\Vert _{1}=\sum_{i=1}^{n}\left|d_{i}s_{i}\right|\\ \textrm{oder }\left\Vert s\right\Vert & = & \left\Vert Ds\right\Vert _{\infty}=\max_{1\leq i\leq n}\left|d_{i}s_{i}\right|.\end{eqnarray*} Diese ergeben ein QP der Form $\min\Phi_{k}\left(s\right)$ s.d. $\left\Vert s\right\Vert \leq\Delta$, siehe Abschnitt \ref{sub:TrustRegionQP}. Allgemein ist es schwierig, eine geeignete Matrix $D$ zu finden. \end{remark} \chapter{Quadratische Programmierung} \section{Definition, Motivation, Beispiele} \[ \left(\textrm{QP}\right)\qquad\begin{array}{lllll} \min & f\left(x\right) & \equiv g^{T}x+\frac{1}{2}x^{T}Bx & \quad & \textrm{(quadr. Zielfunktion)}\smallskip\\ \textrm{s.d.} & a_{i}^{T}x=b_{i} & \textrm{für }i\in\mathcal{G} & & \textrm{(lineare Gleichheitsrestriktionen)}\smallskip\\ & a_{i}^{T}x\leq b_{i} & \textrm{für }i\in\mathcal{U} & & \textrm{(lineare Ungleichungen)}\end{array}\] Ein (QP)\index{QP}\index{Quadratische Programmierung} heißt konvex, falls $B\succeq0$ und streng konvex, falls sogar $B\succ0$. \textbf{Frage:} Warum nicht $a_{i}^{T}x\leq b_{i}$ und $-a_{i}^{T}x\leq-b_{i}$ statt $a_{i}^{T}x\leq b_{i}$? \textbf{Antwort:} Da beide Ungleichungen kolineare Gradienten $a_{i}$ und $-a_{i}$ haben, geht auf jeden Fall LICQ (Linear Independent Constraint Qualification) verloren. $\Rightarrow$ Unbestimmte Multiplikatoren \textbf{Frage:} Warum nicht $A_{\mathcal{G}}x=b_{\mathcal{G}}$ mit $A_{\mathcal{G}}=\left(a_{i}^{T}\right)_{i\in\mathcal{G}}$, $b_{\mathcal{G}}=\left(b_{i}\right)_{i\in\mathcal{G}}$ benutzen, um mit dem impliziten Funktionentheorem einen Teil der Variablen $x_{i}$ zu eliminieren? \textbf{Antwort:} Zerstört typischerweise Dünnbesetztheit, falls vorhanden! \subsection{Orthogonale Projektion} $B=I$, $g=-x_{*}$. Dann ist \[ f\left(x\right)=-x_{*}^{T}x+\frac{1}{2}\left\Vert x\right\Vert ^{2}=\frac{1}{2}\left\Vert x-x_{*}\right\Vert ^{2}-\underbrace{\frac{1}{2}\left\Vert x_{*}\right\Vert ^{2}}_{\textrm{konstant}}\] \[ M\equiv\left\{ x\in\mathbb{R}^{n}\left|a_{i}^{T}x=b_{i}\textrm{ für }i\in\mathcal{G},\, a_{i}^{T}x\leq b_{i}\textrm{ für }i\in\mathcal{U}\right.\right\} \] Die Lösung ist nur unwesentlich einfacher als allgemeines konvexes QP. Anwendung von projezierten Gradienten auf konvexe QP macht also wenig Sinn. \subsection{\label{sub:TrustRegionQP}Trust Regions bezüglich $l_{\infty}$ oder $l_{1}$ Norm.} $g$ und $B$ beliebig und für $i=1,\ldots,n$ $l_{\infty}$-Norm ($2n$ Ungleichungen):\[ \begin{array}{rcl} e_{i}^{T}x & \leq & b_{i}+\Delta\\ -e_{i}^{T}x & \leq & -b_{i}+\Delta\end{array}\quad\Leftrightarrow\quad\left|x_{i}-b_{i}\right|\leq\Delta\] $l_{1}$-Norm ($2^{n}$-Ungleichungen):\[ \sum_{i=1}^{n}d_{i}\left(e_{i}^{T}x-b_{i}\right)\leq\Delta\] für beliebiges $d=\left(d_{1},\ldots,d_{n}\right)\in\left\{ -1,1\right\} ^{n}$ Warnung: Sphärisches Trust Region Problem bzgl. $\left\Vert x\right\Vert =\left\Vert x\right\Vert _{2}$\[ \min f\left(x\right)\quad\textrm{s.d. }\left\Vert x\right\Vert _{2}^{2}\equiv\sum_{i=1}^{n}x_{i}^{2}\leq\Delta\] ist kein QP, da dieses keine quadratischen Nebenbedingungen zulässt. \subsection{QP als lokales Modell für $\min f\left(x\right)$, s.d. $c\left(x\right)=0$} In Newtonähnlicher Form bestimme Schritt $s$ am Punkt $x$ als Lösung des QP\[ \min\underbrace{\nabla f\left(x\right)^{T}s+\frac{1}{2}s^{T}Bs}_{\textrm{quadratisches Zielfunktional}}\quad\textrm{s.d. }\underbrace{c\left(x\right)+\nabla c\left(x\right)s=0}_{\textrm{linearisierte Nebenbedingungen}}.\] \textbf{Frage:} Können wir $B=\nabla^{2}f\left(x\right)$ setzen, d.h. das Zielfunktional durch Taylorreihe zweiter Ordnung von $f\left(x\right)$ selbst ersetzen? \textbf{Antwort 0:} Falls $m=\dim\left(c\right)=n$, d.h. wir haben $m=n$ aktive Beschränkungen, spielt $B$ keine Rolle, da lediglich Zulässigkeit durch $c\left(x\right)=0$ herzustellen ist. (Ja!) Linearisierte Nebenbedingungen ergeben quadratische Konvergenz unter den üblichen Voraussetzungen, d.h. $\nabla c\left(x\right)=c'\left(x\right)$ ist Lipschitz und regulär. \textbf{Antwort 1:} Wenn $m=0$, das Problem also unbeschränkt, ergibt $B\equiv\nabla^{2}f\left(x\right)$ Newton's Methode zur Lösung von $\nabla f\left(x\right)=0$ und damit wiederum quadratische Konvergenz. (Ja!) \textbf{Antwort 2:} Nein! Aus folgendem Beispiel ergibt sich, daß $B=\nabla^{2}f\left(x\right)$ im allgemeinen Fall nur lineare Konvergenz erzielen kann (Mit anderen Worten: Die 2. Ableitungen (Krümmungen) der Restriktionen können nicht vernachlässigt werden, sondern müssen in $B$ eingehen.). \begin{example*} ~\begin{eqnarray*} & \min & \frac{1}{2}\left(x_{1}-1\right)^{2}+\frac{1}{2}\left(x_{2}-1\right)^{2}\\ & \textrm{s.d.} & \sqrt{x_{1}^{2}+x_{2}^{2}}-1=0\end{eqnarray*} Lokales QP mit $B=\nabla^{2}f\left(x\right)$:\begin{eqnarray*} & \min & \frac{1}{2}\left(x_{1}+s_{1}-1\right)^{2}+\frac{1}{2}\left(x_{2}+s_{2}-1\right)^{2}\\ & \textrm{s.d.} & \sqrt{x_{1}^{2}+x_{2}^{2}}-1+\frac{x_{1}s_{1}+x_{2}s_{2}}{\sqrt{x_{1}^{2}+x_{2}^{2}}}=0\end{eqnarray*} Linearisierte Restriktionen entsprechen einer Tangente am Kreis. Leicht (elementar, aber etwas mühsam) zu zeigen: Die resulierenden Iterierten $x^{\left(k+1\right)}=x^{\left(k\right)}+s^{\left(k\right)}$ konvergieren nur linearen gegen $\left(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\right)$. \end{example*} \begin{remark} Auch lineare Konvergenz ist keinesfalls sicher, die linearisierten Nebenbedingungen können sogar inkonsistent sein, so daß das lokale QP eine leere zulässige Menge hat. Es muß daher \emph{relaxiert} werden. Diese Beobachtung betrifft nur die Restriktionen und gilt deshalb für beliebiges $B$. \end{remark} \subsubsection*{{}``Herleitung'' des {}``richtigen'' $B$} Wir wollen schnelle, möglichst quadratische, Konvergenz zur Lösung der KKT Gleichung\[ \nabla f\left(x\right)+\sum_{i=1}^{m}\lambda_{i}\nabla c_{i}\left(x\right)=0,\qquad c_{i}\left(x\right)=0\textrm{ für }i=1,\ldots,m,\] immer unter der Voraussetzung von Gleichheitsnebenbedingungen. Nach einem Satz von Dennis \& Mor\'e konvergiert eine Folge $\left(x_{k},\lambda_{k}\right)$ genau dann superlinear gegen eine Lösung $\left(x_{*},\lambda_{*}\right)$, wenn die Korrekturen $\left(\Delta x_{k},\Delta\lambda_{k}\right)=\left(x_{k}-x_{k-1},\lambda_{k}-\lambda_{k-1}\right)$ asymptotisch gegen die Newtonschritte $\left(s_{k},\sigma_{k}\right)$ konvergieren, d.h. \[ \lim_{k\rightarrow\infty}\frac{\left\Vert \Delta x_{k}-s_{k}\right\Vert +\left\Vert \Delta\lambda_{k}-\sigma_{k}\right\Vert }{\left\Vert \Delta x_{k}\right\Vert +\left\Vert \Delta\lambda_{k}\right\Vert }=0.\] Die $\left(s_{k},\sigma_{k}\right)$ sind Lösung von\[ \left[\begin{array}{cc} \nabla_{x}^{2}L\left(x_{k},\lambda_{k}\right) & \nabla c\left(x_{k}\right)^{T}\\ \nabla c\left(x_{k}\right) & 0\end{array}\right]\left[\begin{array}{c} s_{k}\\ \sigma_{k}\end{array}\right]=\left[\begin{array}{c} -\nabla f\left(x_{k}\right)-\sum_{i=1}^{m}\lambda_{i}\nabla c_{i}\left(x_{k}\right)\\ -c\left(x_{k}\right)\end{array}\right]=\left[\begin{array}{c} -\nabla_{x}L\left(x_{k},\lambda_{k}\right)\\ -c\left(x_{k}\right)\end{array}\right]\] ($L\left(x,\lambda\right)=f\left(x\right)+\sum_{i=1}^{m}\lambda_{i}c_{i}\left(x\right)$). Für $x=x_{k}$, $s=s_{k}$, $g=\nabla f\left(x_{k}\right)$ und $\hat{\lambda}=\lambda_{k}+\sigma_{k}$ ergibt sich äquivalenterweise mit $B=\nabla_{x}^{2}L\left(x,\lambda\right)$, $A=\nabla c\left(x\right)$\[ \left[\begin{array}{cc} B & A^{T}\\ A & 0\end{array}\right]\left[\begin{array}{c} s\\ \hat{\lambda}\end{array}\right]=\left[\begin{array}{c} -g\\ -c\left(x\right)\end{array}\right],\] wegen $A^{T}\sigma+A^{T}\lambda=A^{T}\hat{\lambda}$, d.h. $\sum_{i=1}^{n}\lambda_{i}\nabla c_{i}$ kann von rechts nach links gebracht werden und $\hat{\lambda}$ ersetzt $\sigma$. Falls $B\succ0$, so ist dies äquivalent zu der Aussage, daß $s$ Minimalpunkt des folgenden QP ist\begin{eqnarray*} & \min & g^{T}s+\frac{1}{2}s^{T}Bs\\ & \textrm{s.d.} & c+As=0\end{eqnarray*} D.h. $B=\nabla_{x}^{2}L\left(x,\lambda\right)$ für {}``geeignetes'' $\lambda\in\mathbb{R}^{m}$ ist richtiges $B$. D.h. die Krümmungen $\nabla^{2}c$ der Nebenbedingungen werden gewichtet in den quadratischen Term des lokalen QP Modelles eingebracht. \begin{remark} Bei SLP (Sukzessives LP) wird das Zielfunktional entweder überhaupt nicht verändert, oder linearisiert, Konvergenz ist deshalb fast immer nur linear. \end{remark} \begin{remark} Auch bei NLPs mit Ungleichungen schien lokale Modellierung mit linearisierten Restriktionen und $B=\nabla_{x}^{2}L$ bislang optimal. (Jetzt wird vermutet, daß die volle Lösung des lokalen QP mit Ungleichungen häufig zu aufwendig ist.) \end{remark} \section{Bezeichnungen und Schrittberechnung} \begin{notation*} \emph{Aktive Menge}\index{aktive Menge}\[ A\left(x\right)=\mathcal{G}\cup\left\{ i\in\mathcal{U}\left|a_{i}^{T}x=b_{i}\right.\right\} \] \end{notation*} \begin{notation*} Die \emph{Arbeitsmenge\index{Arbeitsmenge}} \[ W\subseteq A\left(x\right)\quad\textrm{mit }\mathcal{G}\subseteq W\] ergibt sich nicht eindeutig aus $x$, sondern dem Iterationsverlauf. $m:=\left|W\right|$. \end{notation*} \begin{notation*} $A_{W}\in\mathbb{R}^{m\times n}$, $b_{W}\in\mathbb{R}^{m}$,\begin{eqnarray*} A_{W} & := & \left\{ a_{i}\right\} _{i\in W}\\ b_{W} & := & \left(b_{i}\right)_{i\in W}\end{eqnarray*} Wegen $W\subseteq A\left(x\right)$ gilt $A_{W}x-b_{W}=0$. \end{notation*} \begin{notation*} $Z_{W}\in\mathbb{R}^{\left(n-m\right)\times n}$, so daß \[ \det\left(\left[\begin{array}{c} A_{W}\\ Z_{W}\end{array}\right]\right)\neq0,\quad A_{W}Z_{W}^{T}=0,\, Z_{W}A_{W}^{T}=0.\] $Z_{W}$ existiert, vorausgesetzt alle auftretenden $A_{W}$ haben vollen Rang $m$. \end{notation*} \begin{notation*} $Y_{W}\in\mathbb{R}^{m\times n}$, so daß \[ \det\left(\left[\begin{array}{c} Y_{W}\\ Z_{W}\end{array}\right]\right)\neq0\] $\Rightarrow$ $\det\left(A_{W}Y_{W}^{T}\right)\neq0$. \end{notation*} \subsection{Nullraumverfahren zur Berechnung von $s_{W}$} Der Vektor\[ s_{W}=\argmin\left\{ \left(g+Bx\right)s+\frac{1}{2}s^{T}Bs\left|A_{W}s=\hat{b}\right.\right\} \] löst nach KKT die Gleichung\begin{equation} \left[\begin{array}{cc} B & A_{W}^{T}\\ A_{W} & 0\end{array}\right]\left[\begin{array}{c} s\\ \hat{\lambda}\end{array}\right]=\left[\begin{array}{c} -\left(g+Bx\right)\\ \hat{b}\end{array}\right]\label{eq:NullraumverfahrenKKT}\end{equation} Sei $\hat{g}=g+Bx$. Die Zerlegung\[ s_{W}=Y_{W}^{T}s_{Y}+Z_{W}^{T}s_{Z}\] ergibt durch Einsetzen in die zweite Gleichung von (\ref{eq:NullraumverfahrenKKT}) folgendes Gleichungssystem zur Bestimmung von $s_{Y}$:\begin{equation} A_{W}Y_{W}^{T}s_{Y}+\underbrace{A_{W}Z_{W}^{T}}_{=0}s_{Z}=\hat{b}.\label{eq:Nullraum_sY}\end{equation} Einsetzen in die erste Gleichung von (\ref{eq:NullraumverfahrenKKT}) ergibt\[ BY_{W}^{T}s_{Y}+BZ_{W}^{T}s_{Z}+A_{W}^{T}\hat{\lambda}_{W}=-\hat{g}.\] Multiplikation von links mit $Z_{W}$ ergibt ein Gleichungssystem zur Bestimmung von $s_{Z}$, \begin{equation} Z_{W}BZ_{W}^{T}s_{Z}=-Z_{W}\hat{g}-Z_{W}BY_{W}^{T}s_{Y}-\underbrace{Z_{W}A_{W}^{T}}_{=0}\hat{\lambda}_{W},\label{eq:Nullraum_sZ}\end{equation} während man durch Multiplikation von links mit $Y_{W}$ ein Gleichungssystem zur Bestimmung von $\hat{\lambda}_{W}$ erhält:\begin{equation} Y_{W}A_{W}^{T}\hat{\lambda}_{W}=-Y_{W}\hat{g}-Y_{W}BY_{W}^{T}s_{Y}-Y_{W}BZ_{W}^{T}s_{Z}.\label{eq:Nullraum_lW}\end{equation} Gelöst werden müssen also die linearen Systeme (\ref{eq:Nullraum_sY}), (\ref{eq:Nullraum_sZ}) und (\ref{eq:Nullraum_lW}) mit den Matrizen $A_{W}Y_{W}^{T}\in\mathbb{R}^{m\times m}$, $Z_{W}BZ_{W}^{T}\in\mathbb{R}^{\left(n-m\right)\times\left(n-m\right)}$ und $Y_{W}A_{W}^{T}\in\mathbb{R}^{m\times m}$. Dazu bieten sich QR und eine LV basierte Faktorisierungvariante an. \subsubsection*{QR} \begin{eqnarray*} A_{W} & = & \underbrace{\left[\begin{array}{cc} * & 0\end{array}\right]}_{U_{W}}\left[\begin{array}{c} Y_{W}\\ Z_{W}\end{array}\right]\qquad U_{W}\textrm{ eine untere Dreiecksmatrix}\\ \Leftrightarrow\, A_{W}^{T} & = & \underbrace{\left[\begin{array}{cc} Y_{W}^{T} & Z_{W}^{T}\end{array}\right]}_{Q}\underbrace{\left[\begin{array}{c} *\\ 0\end{array}\right]}_{R}\qquad R\textrm{ eine obere Dreiecksmatrix}\end{eqnarray*} $Q^{T}Q=QQ^{T}=I$ $\Rightarrow$ $Y_{W}Z_{W}^{T}=0$, $Y_{W}Y_{W}^{T}=I_{m}$, $ZZ^{T}=I_{n-m}$ Für die Matrizen in den linearen Systemen (\ref{eq:Nullraum_sY}), (\ref{eq:Nullraum_sZ}) und (\ref{eq:Nullraum_lW}) ist dann \begin{itemize} \item $A_{W}Y_{W}^{T}=U_{W}$ unterhalbs dreiecksformig und $\det\left(A_{W}\right)=\det\left(U_{W}\right)$, \item $Y_{W}A_{W}^{T}=U_{W}^{T}$ oberhalb dreiecksformig, und \item $Z_{W}BZ_{W}^{T}$ {}``normalerweise'' positiv definit und hat deshalb eine Cholesky-Faktorisierung $Z_{W}^{T}BZ_{W}=LDL^{T}$ mit $D\succ0$ diagonal. \end{itemize} \textbf{Vorteil}: Numerische Stabilität ohne jegliche Pivotierung. \textbf{Nachteil}: Der Speicherbedarf für $U_{W}$, $Y_{W}$, $Z_{W}$ und $B$ ist $\frac{1}{2}m^{2}+n^{2}+\frac{1}{2}n^{2}=\frac{1}{2}m^{2}+\frac{3}{2}n^{2}$, während der Speichbedarf von $A$ und $B$ höchstens $mn+\frac{1}{2}n^{2}$ beträgt, wobei sich die Dünnbesetztheit in $A_{W}$ nicht auf die Faktoren $U_{W}$, $Y_{W}$, $Z_{W}$ und schon garnicht auf die Projektionen $Z_{W}BZ_{W}^{T}$ usw. überträgt. \subsubsection*{LV} \[ P_{W}AP_{W}'=L_{W}\left[\begin{array}{cc} U_{W} & N_{W}\end{array}\right]\] mit unterer Dreiecksmatrix $L_{W}$, oberer Dreiecksmatrix $U_{W}$, Zeilenpermutationen $P_{W}$ und Spaltenpermutationen $P_{W}'$.\[ \Rightarrow\,\left[\begin{array}{c} Y_{W}\\ Z_{W}\end{array}\right]=\left[\begin{array}{cc} I_{m} & 0\\ -N_{W}^{T}U_{W}^{-T} & I_{n-m}\end{array}\right].\] Probe:\begin{eqnarray*} A_{W}Z_{W}^{T} & = & L_{W}\left[\begin{array}{cc} U_{W} & N_{W}\end{array}\right]\left[\begin{array}{c} -U_{W}^{-1}N_{W}\\ I\end{array}\right]=L_{W}\cdot0=0\\ A_{W}Y_{W}^{T} & = & L_{W}U_{W}\textrm{ ist regulär}\end{eqnarray*} Ausdruck für $Z_{W}BZ_{W}^{T}$ ist\[ \left[\begin{array}{cc} -N_{W}^{T}U_{W}^{-T} & I\end{array}\right]\left[\begin{array}{cc} B_{11} & B_{12}\\ B_{12}^{T} & B_{22}\end{array}\right]\left[\begin{array}{c} -U_{W}^{-1}N_{W}\\ I\end{array}\right]=B_{22}-N_{W}^{T}U_{W}^{-T}B_{12}-B_{12}^{T}U_{W}^{-1}N_{W}+N_{W}^{T}U_{W}^{-T}B_{11}U_{W}^{-1}N_{W}\] \textbf{Vorteil}: $Y_{W}=I$ ist trivial und $Z_{W}$ braucht nicht explizit berechnet zu werden. Der Speicherbedarf von $U_{W}$, $L_{W}$, $N_{W}$ und $B$ beträgt $4m+\frac{1}{2}n^{2}$ und ist im vollbesetzen Fall minimal. Im dünnbesetzten Fall läßt sich durch geeignete Pivotisierung der Speicherplatzbedarf eventuell drastisch reduzieren. Das gilt dann auch für die Anzahl der arithmetischen Operationen. \textbf{Nachteil}: Wenn wegen Änderungen in $W$ Zeilen von $A_{W}$ herausfallen oder hinzukommen, muß schon aus Gründen der Stabilität eventuell aufwendig nachpivotisiert werden (siehe \cite{KiSc}). \begin{remark} Vorteil beider Optionen (QR und LV) ist, daß Hinzunahme oder Entfernung eines Elementes in $W$ mit Aufwand der Ordnung $n^{2}$ auf Ebene der Faktorisierungen nachvollzogen werden können. Der Gesamtaufwand für die QP-Lösung ist deshalb $O\left(n^{2}\cdot\textrm{Schrittzahl}\right)$. Die Schrittzahl kann theoretisch exponentiell wachsen (z.B. wie $2^{n}$, siehe Klee-Winty). Folgendes Verfahren ist eine Verallgemeinerung des Simplex-Verfahrens. \end{remark} \section{Aktive Mengen Methode für konvexes QP} \begin{algorithm} [Skizze]Ausgehend von $x_{0}\in M$ mit $W_{0}\subset A\left(x_{0}\right)$ erzeuge eine Folge \[ x_{k+1}=x_{k}+\alpha_{k}s_{k}\] wobei \[ s_{k}=s_{W_{k}}\left(x_{k}\right)\] und jeweils gilt $\alpha_{k}<1$, falls der volle Schritt aus $M$ führt. Dann wird $m_{k+1}=\left|W_{k+1}\right|>\left|W_{k}\right|=m_{k}$, da mindestens eine Restriktion aktiv wird. $\alpha_{k}=1$, falls der Schritt in $M$ bleibt. Dann werden die Multiplikatoren $\hat{\lambda}_{k}$ geprüft. Gilt $e_{i}^{T}\hat{\lambda}_{k}\geq0$ für alle $i\in W_{k}\cap\mathcal{U}$, so ist ein Optimum erreicht. Andernfalls wird ein solches $i$ aus $W$ entfernt. \end{algorithm} Zu zeigen: \begin{enumerate} \item Die Berechnung von $\alpha_{k}$ (einfach). \item $s_{k}$ ist Abstiegsrichtung, falls $s_{k}^{T}A_{W_{k}}=b_{k}$ $\Leftrightarrow$ $\hat{b}_{k}=0$. \item $W_{k+1}=W_{k}\setminus\left\{ i\right\} $ ergibt eine Abstiegsrichtung, wenn $e_{i}^{T}\hat{\lambda}_{k}<0$. \end{enumerate} \begin{enumerate} \item ~\[ \alpha_{k}=\min_{\substack{i\in U\setminus W_{k}\\ a_{i}^{T}s_{k}>0} }\left(1,\frac{b_{i}-a_{i}^{T}x}{a_{i}^{T}s_{k}}\right)\] \item ~\[ \left[\begin{array}{cc} B & A_{W_{k}}^{T}\\ A_{W_{k}} & 0\end{array}\right]\left[\begin{array}{c} s_{k}\\ \hat{\lambda}_{k}\end{array}\right]=-\left[\begin{array}{c} \hat{g}_{k}\\ 0\end{array}\right],\] $\hat{g}_{k}=g+Bx_{k}$ der aktuelle Gradient. Multiplikation von links mit $\left[\begin{array}{cc} s_{k}^{T} & 0\end{array}\right]\in\mathbb{R}^{n+n}$ ergibt\[ s_{k}^{T}Bs_{k}+\underbrace{s_{k}^{T}A_{W_{k}}}_{=0}\hat{\lambda}_{k}=-\hat{g}_{k}^{T}s_{k},\] d.h. $\hat{g}_{k}^{T}s_{k}=-s_{k}^{T}Bs_{k}$. Unter der Voraussetzung $B\succ0$ oder zumindest $z_{k}Bz_{k}^{T}\succ0$ für alle vorkommenden $W_{k}$ ist $\hat{g}_{k}^{T}s_{k}<0$. \item $x=x_{k+1}$ ist Lösung des Problems\[ \min\left\{ g^{T}x+\frac{1}{2}x^{T}Bx\left|A_{W}x=b_{W}\right.\right\} \] und erfüllt deshalb die KKT-Bedingung \[ g+Bx+A_{W}^{T}\hat{\lambda}_{W}=0\] mit Rang$\left(A_{W}\right)=m=\left|W_{k}\right|$. Falls ein $\hat{\lambda}_{i}=e_{i}^{T}\hat{\lambda}_{W}<0$ ist, setze $\hat{W}=W\setminus\left\{ i\right\} $ und finde eine Lösung $s$ mit $A_{\hat{W}}s=0$, aber $a_{i}^{T}s=-1$. Dann gilt \[ \left(g+Bx\right)^{T}s=-\hat{\lambda}_{W}^{T}A_{W}s=-\left(-\hat{\lambda}_{i}\right)=\hat{\lambda}_{i}<0,\] d.h. eine Lösung des durch $\hat{W}$ definierten Problemes führt zu einer weiteren Reduktion. \end{enumerate} Es gibt zwei Arten von Schritten, je nachdem ob eine bindende (blockierende) Restriktion erreicht wird oder nicht:\[ x_{k}+s_{k}\in\argmin\left\{ f\left(x\right)\left|a_{i}^{T}x=b_{i},i\in W_{k}\right.\right\} .\] Falls $\left\{ j\right\} $ bindend ist, so wird $W_{k+1}=W_{k}\cup\left\{ j\right\} $ und ein neuer Schritt berechnet. Diese Erweiterung von $W_{k}$ kann höchstens einmal in Folge passieren, da die Gradienten $\left\{ a_{i}\right\} _{i\in W_{k}}\subset\left\{ a_{i}\right\} _{i\in A\left(x_{k}\right)}$ nach Voraussetzung linear unabhängig sind. Wird ein $x_{k+1}=x_{k}+s_{k}$ mit $\alpha_{k}=1$ erreicht, so ist $x_{k+1}$ entweder \underbar{die} Lösung des konvexen QP (da $\det B\neq0$) oder es gilt:\begin{eqnarray*} \nabla f\left(x\right)+\lambda^{T}\nabla c\left(x\right) & = & 0\\ g\left(x\right)+B^{T}x+\lambda^{T}A_{W}^{T} & = & 0\qquad\left(\textrm{mit }\lambda_{i}\geq0\textrm{ falls }g_{i}\left(x\right)\leq0\right),\end{eqnarray*} \emph{aber} mindestens ein Lagrange-Multiplikator $\lambda_{i}$ mit $i\in\mathcal{U}$ ist negativ. Dann wird $W_{k+1}=W_{k}\setminus\left\{ i\right\} $ gesetzt und der nächste Schritt berechnet. \textbf{Problem}: Zumindest theoretisch ist es möglich, daß wiederholt $\alpha_{k}=0$ ist und somit keine Reduktion in $f$ erzielt wird. Dabei können Restriktionen wiederholt in die aktive Menge $W_{k}$ eingeführt und wieder rausgenommen werden. Dieses {}``cycling'' Phänomen kann schon in LP auftreten. Es wird meist von Algorithmen ignoriert und durch Rundungsfehler unterbrochen. \section{QP Algorithmus nach W\&N} \begin{algorithm} \label{alg:QPnachWN}~ \begin{algorithmic} \STATE Bestimme zulässigen Anfangspunkt $x_{0}$, setze $W_{0}\subseteq A\left(x_{0}\right)$ mit $\mathcal{G}\subseteq W_{0}$. \FOR{$k=0,1,\ldots$} \STATE Berechne $s_{k}=s_{W_{k}}$ und $\hat{\lambda}_{k}=\hat{\lambda}_{W_{k}}$\[ \left[\begin{array}{cc} B & A_{W_{k}}^{T}\\ A_{W_{k}} & 0\end{array}\right]\left[\begin{array}{c} s_{W_{k}}\\ \hat{\lambda}_{W_{k}}\end{array}\right]=-\left[\begin{array}{c} g+Bx_{k}\\ \left(b-Ax_{k}\right)_{W_{k}}\end{array}\right].\] \IF{$s_{k}=0$} \STATE Betrachte $\hat{\lambda}_{k}$. \IF{$e_{i}^{T}\hat{\lambda}_{k}\geq0$ für alle $i\in W_{k}\cap\mathcal{U}$} \STATE Stop mit $x_{*}=x_{k}$ \ELSE \STATE $i_{*}\equiv\argmin\left\{ e_{i}^{T}\hat{\lambda}_{k}\left|i\in W_{k}\cap\mathcal{U}\right.\right\} $, $W_{k+1}=W_{k}\setminus\left\{ i_{*}\right\} $, $x_{k+1}=x_{k}$\ENDIF \ELSE \STATE Berechne $\alpha_{k}$ wie oben spezifiziert \STATE Setze $x_{k+1}=x_{k}+\alpha_{k}s_{k}$ \IF{$\alpha_{k}<1$} \STATE Setze $W_{k+1}=W_{k}\cup\left\{ i\right\} $ für bindende Restriktion $i\in\mathcal{U}\setminus W_{k}$ mit $a_{i}^{T}x_{k+1}=b_{i}$ \ENDIF \ENDIF \ENDFOR \end{algorithmic} \end{algorithm} \begin{thm} Vorausgesetzt, alle $A_{k}\equiv A_{W_{k}}$ haben vollen Rank und alle $\alpha_{k}>0$ für $k=1,2,\ldots$. Dann konvergiert Algorithmus \ref{alg:QPnachWN} in endlich vielen Schritten zum Minimalpunkt des streng konvexen QP\begin{eqnarray*} & \min & f\left(x\right)\equiv g^{T}x+\frac{1}{2}x^{T}Bx\qquad\left(B\succ0\right)\\ & \textrm{s.d.} & a_{i}^{T}x=b_{i}\qquad\textrm{für }i\in\mathcal{G}\\ & & a_{i}^{T}x\leq b_{i}\qquad\textrm{für }i\in\mathcal{U}\end{eqnarray*} \end{thm} \begin{proof} Wegen vorausgesetzten $\alpha_{k}>0$ bilden die $f\left(x_{k}\right)$ für $k=1,\ldots,n$ eine streng monoton fallende Folge. Mindestens jede $n$-te Iterierte ist eindeutiger Minimalpunkt eines Problemes $\min f\left(x\right)$ s.d. $a_{i}^{T}x=b_{i}$ für $i\in W$. Es gibt aber maximal $2^{\left|\mathcal{U}\right|}$ unterschiedliche Mengen $W$. Also ist die Gesamtzahl der Schritte durch $n\cdot2^{\left|\mathcal{U}\right|}$ beschränkt und damit endlich. \end{proof} \begin{remark} Wegen des Klee\&Minty Beispiels kann die Schrittzahl des Simplex-Algorithmus und damit auch des Algorithmus \ref{alg:QPnachWN} tatsächlich exponentiell mit der Zahl der Variablen bzw. Beschränkungen wachen. In der Praxis erwartet man eher $O\left(n\left|\mathcal{U}\right|\right)$ Schritte. \end{remark} \textbf{Frage}: Wie läßt sich ein zulässiger Anfangspunkt $x_{0}\in\mathbb{R}$ bestimmen? \textbf{Antwort}: Phase I Prozedur, d.h. Lösung eines Hilfs-LP oder -QP\begin{eqnarray*} & \min & f\left(z\right)=\sum_{i=1}^{m}z_{i}^{q}\qquad m=\left|\mathcal{G}\right|+\left|\mathcal{U}\right|,\, q\in\left\{ 1,2\right\} \\ & \textrm{s.d.} & z_{i}\geq0\qquad i=1,\ldots,m\\ & & a_{i}^{T}x_{i}+z_{i}\textrm{sign}\left(b_{i}\right)=b_{i}\qquad i\in\mathcal{G}\\ & & a_{i}^{T}x_{i}-z_{i}\leq b_{i}\qquad i\in\mathcal{U}.\end{eqnarray*} \begin{remark} Ein zulässiger Punkt für dieses relaxierte Problem ist $x=0$, $z_{i}=\left|b_{i}\right|$ für $i\in\mathcal{G}$, $z_{i}=\max\left\{ 0,-b_{i}\right\} $ für $i\in\mathcal{U}$. Die Iterierte $\left(x_{k},z_{k}\right)$ liefert für das ursprüngliche QP einen zulässigen Punkt, gdw. $f\left(z_{k}\right)=f\left(0\right)=0$ erreicht wird. \end{remark} \chapter{Sequentielle Quadratische Programmierung (SQP)} \begin{equation} \min f\left(x\right)\qquad\textrm{s.d.}\quad c\left(x\right)=0\label{eq:SQPNLP}\end{equation} mit $f,g\in C^{2,1}$, d.h. 2-mal stetig differenzierbar mit Lipschitz-stetiger 2. Ableitung, $x\in\mathbb{R}^{n}$. \textbf{Idee}:\index{SQP} Modelliere das Problem am aktuellen Iterationspunkt mittels eines QP-Subproblems und benutze dessen Lösung als nächsten Iterationspunkt. Herausforderung: Die Wahl des QP-Subproblemes. \section{Basis Modell von SQP} Bezeichnungen: \begin{eqnarray*} L\left(x,\lambda\right) & = & f\left(x\right)+\lambda^{T}c\left(x\right)\\ B\left(x,\lambda\right) & = & \nabla_{xx}^{2}L\left(x,\lambda\right)\\ A\left(x\right)^{T} & = & \left(\nabla c_{1}\left(x\right),\nabla c_{2}\left(x\right),\ldots,\nabla c_{m}\left(x\right)\right)\end{eqnarray*} Für einen Iterationspunkt $\left(x_{k},\lambda_{k}\right)$,\begin{eqnarray*} B_{k} & := & B\left(x_{k},\lambda_{k}\right)\\ A_{k} & := & A\left(x_{k},\lambda_{k}\right).\end{eqnarray*} Sei $p_{k}$ eine Lösung von\index{SQP!Basismodell} \begin{eqnarray} & \min_{p} & \nabla f\left(x_{k}\right)p+\frac{1}{2}p^{T}B_{k}p\label{eq:SQPProblem}\\ & \textrm{s.d.} & A_{k}p+c\left(x_{k}\right)=0\nonumber \end{eqnarray} für den aktuellen Iterationspunkt $x_{k}$. Dann\begin{eqnarray*} x_{k+1} & = & x_{k}+p_{k}\\ \lambda_{k+1} & = & ?\end{eqnarray*} Durch die Wahl von $B_{k}$ entstehen verschiedene SQP-Verfahren. \newpage \textbf{Voraussetzungen}: \begin{enumerate} \item LICQ, d.h. $\nabla c_{i}\left(x_{k}\right)$, $i=1,\ldots,m$, sind linear unabhängig. Mit anderen Worten: $A_{k}$ hat vollen Zeilenrank. \item $d^{T}B_{k}d>0$ für alle $d\neq0$ mit $A_{k}^{T}d=0$. \end{enumerate} Das Problem (\ref{eq:SQPProblem}) hat eine eindeutige Lösung $\left(p_{k},\mu_{k}\right)$, welches die KKT-Bedingungen erfüllt:\begin{equation} \left[\begin{array}{cc} B_{k} & A_{k}^{T}\\ A_{k} & 0\end{array}\right]\left[\begin{array}{c} p_{k}\\ \mu_{k}\end{array}\right]=\left[\begin{array}{c} -\nabla f\left(x_{k}\right)\\ -c\left(x_{k}\right)\end{array}\right].\label{eq:SQPKKT}\end{equation} Andererseits: Das KKT von (\ref{eq:SQPNLP}) ist\begin{eqnarray*} F\left(x,\lambda\right):=\nabla f\left(x\right)+A^{T}\left(x\right)\lambda & = & 0\\ c\left(x\right) & = & 0\end{eqnarray*} Wenden wir auf $F\left(x,\lambda\right)=0$ das gewöhnliche Newton-Verfahren\index{Newton-Verfahren} im Punkt $\left(x_{k},\lambda_{k}\right)$ an, so erhalten wir\begin{eqnarray*} F\left(x_{k},\lambda_{k}\right)+\frac{\partial F\left(x_{k},\lambda_{k}\right)}{\partial\left(x,\lambda\right)}\left(\begin{array}{c} \Delta x_{k}\\ \Delta\lambda_{k}\end{array}\right) & = & 0,\end{eqnarray*} d.h.\begin{eqnarray} \left[\begin{array}{cc} B_{k} & A_{k}^{T}\\ A_{k} & 0\end{array}\right]\left[\begin{array}{c} \Delta x_{k}\\ \Delta\lambda_{k}\end{array}\right] & = & \left[\begin{array}{c} -\nabla f\left(x_{k}\right)-A_{k}^{T}\lambda_{k}\\ -c\left(x_{k}\right)\end{array}\right]\nonumber \\ \textrm{bzw. }\left[\begin{array}{cc} B_{k} & A_{k}^{T}\\ A_{k} & 0\end{array}\right]\left[\begin{array}{c} \Delta x_{k}\\ \lambda_{k}+\Delta\lambda_{k}\end{array}\right] & = & \left[\begin{array}{c} -\nabla f\left(x_{k}\right)\\ -c\left(x_{k}\right)\end{array}\right]\label{eq:SQPNewton}\end{eqnarray} Wenn $p_{k}=\Delta x_{k}$ und $\mu_{k}=\lambda_{k+1}=\lambda_{k}+\Delta\lambda_{k}$, so ist (\ref{eq:SQPKKT}) äquivalent zu (\ref{eq:SQPNewton}). Zeigen: (\ref{eq:SQPProblem}) ist äquivalent zu \begin{eqnarray*} & \min_{p} & \nabla_{x}L\left(x_{k},\lambda_{k}\right)p+\frac{1}{2}p^{T}B_{k}p\\ & \textrm{s.d.} & A_{k}p+c_{k}=0\end{eqnarray*} Es ist\[ \nabla_{x}L\left(x_{k},\lambda_{k}\right)p=\nabla f\left(x_{k}\right)^{T}p+\left(A_{k}^{T}\lambda_{k}\right)^{T}p=\nabla f\left(x_{k}\right)^{T}p-\underbrace{\lambda_{k}^{T}c_{k}}_{\textrm{unabhängig von }x_{k}}.\] \subsection{Ungleichungsrestriktionen} \begin{eqnarray*} & \min & f\left(x\right)\\ & \textrm{s.d.} & c_{i}\left(x\right)=0\quad i\in\mathcal{G}\\ & & c_{i}\left(x\right)\geq0\quad i\in\mathcal{U}\end{eqnarray*} \begin{enumerate} \item Möglichkeit: Approximation der Arbeitsmenge\[ W_{k}\approx A\left(x_{*}\right):=\left\{ i\in\mathcal{U}\left|c_{i}\left(x\right)=0\right.\right\} \cup\mathcal{G}.\] Also Zurückführung auf (\ref{eq:SQPNLP}). \item Möglichkeit:\begin{eqnarray*} & \min_{p} & \nabla f\left(x_{k}\right)p+\frac{1}{2}p^{T}B_{k}p\\ & & \nabla c_{i}\left(x_{k}\right)^{T}p+c_{i}\left(x_{k}\right)=0\qquad i\in\mathcal{G}\\ & & \nabla c_{i}\left(x_{k}\right)^{T}p+c_{i}\left(x_{k}\right)\geq0\qquad i\in\mathcal{U}\end{eqnarray*} \end{enumerate} \section{Die Wahl von $B_{k}$} \subsection{Quasi-Newton} \[ B_{k}\approx\nabla_{xx}^{2}L\left(x_{k},\lambda_{k}\right)\] mittels z.B. BFGS-Formel\index{BFGS-Aufdatierung}\index{Quasi-Newton}\begin{eqnarray*} B_{k+1} & = & B_{k}-\frac{B_{k}\Delta x_{k}\Delta x_{k}^{T}B_{k}}{\Delta x_{k}B_{k}\Delta x_{k}}+\frac{y_{k}y_{k}^{T}}{\Delta x_{k}^{T}y_{k}}\\ \textrm{mit }y_{k} & = & \nabla_{x}L\left(x_{k+1},\lambda_{k+1}\right)-\nabla_{x}L\left(x_{k},\lambda_{k}\right).\end{eqnarray*} Damit $B_{k+1}\succ0$, braucht man $\Delta x_{k}^{T}y_{k}>0$. Dies ist aber nicht immer erfüllt! Gedämpfte BFGS Aufdatierung\index{BFGS-Aufdatierung!gedämpft}: Man definiert \[ y_{k}=\nabla_{x}L\left(x_{k+1},\lambda_{k+1}\right)-\nabla_{x}L\left(x_{k},\lambda_{k}\right)\] wie zuvor und\begin{eqnarray*} r_{k} & = & \theta_{k}y_{k}+\left(1-\theta_{k}\right)B_{k}\Delta x_{k}\\ \theta_{k} & = & \begin{cases} 1 & \Delta x_{k}^{T}y_{k}\geq\frac{1}{5}\nabla x_{k}^{T}B_{k}\nabla x_{k},\\ \frac{\frac{4}{5}\Delta x_{k}^{T}B_{k}\Delta x_{k}}{\Delta x_{k}^{T}B_{k}\Delta x_{k}-\Delta x_{k}^{T}y_{k}} & \textrm{sonst}\end{cases}\in\left(0,1\right],\end{eqnarray*} also die größte Zahl $\leq1$, so daß $r_{k}^{T}\Delta x_{k}\geq\frac{1}{5}\Delta x_{k}^{T}B_{k}\Delta x_{k}$. \[ B_{k+1}=B_{k}-\frac{B_{k}\Delta x_{k}\Delta x_{k}^{T}B_{k}}{\Delta x_{k}^{T}B_{k}\Delta x_{k}}+\frac{r_{k}r_{k}^{T}}{\Delta x_{k}^{T}r_{k}}.\] Dann ist $\Delta x_{k}^{T}r_{k}\geq\frac{1}{5}\Delta x_{k}^{T}B_{k}\Delta x_{k}>0$, d.h. $B_{k+1}\succ0$. Ist $\theta_{k}=1$, so ist $r_{k}=y_{k}$ (gewöhnliche BFGS). Wäre $\theta_{k}=0$, so wäre $B_{k+1}=B_{k}$, d.h. keine Aufdatierung. Dieses Verfahren ist nicht effektiv, wenn $\nabla_{xx}^{2}L\left(x_{k},\lambda_{k}\right)$ nicht positiv definit ist! \subsection{Symmetrische Rank1 (SR1)-Aufdatierung} \[ B_{k+1}=B_{k}+\frac{\left(y_{k}-B_{k}\Delta x_{k}\right)\left(y_{k}-B_{k}\Delta x_{k}\right)^{T}}{\left(y_{k}-B_{k}\Delta x_{k}\right)^{T}\Delta x_{k}^{T}}.\] $B_{k+1}$ ist nicht unbedingt positiv definit.\index{SR1-Aufdatierung} Gute Wahl für Trust-Region SQP-Verfahren. \subsection{Erweiterte (Augmented) Lagrange Funktion} Für Line-Search SQP-Verfahren benötigt man die positive Definitheit von $B_{k}$. Ersetze die Lagrange-Funktion durch\index{Augmented Lagrange Funktion} \[ L_{A}\left(x,\lambda;\mu\right)=f\left(x\right)+\lambda^{T}c\left(x\right)+\frac{1}{\mu}\left\Vert c\left(x\right)\right\Vert ^{2}\] Dann existiert ein $\mu^{*}$, so daß $\nabla_{xx}^{2}L_{A}\left(x_{k},\lambda_{k};\mu\right)\succ0$ für alle $\mu<\mu^{*}$. Problem: Die Wahl von $\mu^{*}$. \section{Schrittweiten-Bestimmung} Dafür benötigen wir eine Hilfsfunktion, z.B.\[ \varphi_{1}\left(x,\mu,\eta\right)=f\left(x\right)+\frac{1}{\mu}\sum_{i\in\mathcal{G}}\left|c_{i}\left(x\right)\right|+\frac{1}{\eta}\sum_{i\in\mathcal{U}}c_{i}^{+}\left(x\right)\] für $\mu,\eta\in\mathbb{R}_{+}$, $c_{i}^{+}:=\max\left(0,c_{i}\right)$. \begin{remark} $\varphi_{1}$ ist nicht überall differenzierbar. Aber sie ist richtungsdifferenzierbar! \end{remark} \begin{lem} Wenn $\left(p_{k},\lambda_{k+1}\right)$ eine Lösung von (\ref{eq:SQPProblem}) ist, so gilt für die Richtungsableitung von $\varphi_{1}\left(x,\mu,\eta\right)$ in der Richtung $p_{k}$:\[ \nabla\left(\varphi_{1}\left(x_{k},\mu,\eta\right);p_{k}\right)\leq-p_{k}^{T}B_{k}p_{k}-\left(\mu^{-1}-\left\Vert \lambda_{k+1}\right\Vert _{\infty}\right)\left\Vert c\left(x_{k}\right)\right\Vert _{1}\] \end{lem} \begin{proof} \cite{NoWr} \end{proof} \begin{conclusion} $p_{k}$ ist eine Abstiegsrichtung für $\varphi_{1}\left(x_{k},\mu\right)$, wenn $B_{k}\succ0$ und $\mu$ hinreichend klein ist. \end{conclusion} Die Lösung des Problemes\[ \min_{\alpha_{k}\in\mathbb{R}_{+}}\varphi_{1}\left(x_{k}+\alpha_{k}p_{k},\mu,\eta\right)\] liefert dann eine Schrittweite $\alpha_{k}$. Eine andere Hilfsfunktion ist Fletcher's erweiterte Lagrange-Funktion\index{Lagrange Funktion!augmented!Fletcher}\[ \varphi_{F}=L_{A}\left(x,\lambda\right)=L\left(x,\lambda\right)+\frac{1}{2\mu}\left\Vert c_{\mathcal{G}}\left(x\right)\right\Vert ^{2}+\frac{1}{2\eta}\left\Vert c_{\mathcal{U}}^{+}\left(x\right)\right\Vert ^{2}.\] Die Hilfsfunktion ist eine weitere Variationsmöglichkeit, um verschiedene SQP-Verfahren zu erzeugen. \section{Line-Search SQP-Algorithmus} \begin{algorithm} [Line-Search SQP]\index{SQP!Line-Search}Man wähle $\eta\in\left(0,\frac{1}{2}\right)$ und $\tau\in\left(0,1\right)$. \begin{algorithmic}\REQUIRE Startpunkt $\left(x_{0},\lambda_{0}\right)$ und $B_{0}\approx\nabla_{xx}^{2}L\left(x_{0},\lambda_{0}\right)$. \COMPUTE \STATE Man berechne $f\left(x_{0}\right)$, $\nabla f\left(x_{0}\right)$, $c\left(x_{0}\right)$ und $A_{0}=\nabla c\left(x_{0}\right)$ \FOR{$k=0,1,2,\ldots$} \IF{Test} \STATE STOP \ENDIF \STATE Löse (\ref{eq:SQPProblem}) mit Voraussetzung $B_{k}\succ0$ und erhalte $\left(p_{k},\lambda_{k+1}\right)$. \STATE Finde $\mu_{k}$, so daß $p_{k}$ eine Abstiegsrichtung für $\varphi\left(x,\mu\right)$ an der Stelle $x_{k}$ ist. \STATE Setze $\alpha_{k}=1$ \WHILE{$\varphi\left(x_{k}+\alpha_{k}p_{k},\mu_{k},\eta\right)>\varphi\left(x_{k},\mu_{k}\right)+\eta\alpha_{k}D\left(\varphi\left(x_{k},\mu_{k},\eta\right);p_{k}\right)$} \STATE $\alpha_{k}:=\tau_{\alpha}\alpha_{k}$ für $\tau_{\alpha}\in\left(0,\tau\right)$ \ENDWHILE\STATE $x_{k+1}:=x_{k}+\alpha_{k}p_{k}$ \STATE Berechne $f\left(x_{k+1}\right)$, $\nabla f\left(x_{k+1}\right)$, $c\left(x_{k+1}\right)$, $A_{k+1}=\nabla c\left(x_{k+1}\right)$ \IF{$\lambda_{k+1}$ nicht bekannt} \STATE Berechne $\lambda_{k+1}$ aus $\min_{\lambda}\left\Vert \nabla f\left(x_{k+1}\right)+A_{k+1}^{T}\lambda\right\Vert _{2}^{2}$, d.h. $\lambda_{k+1}=-\left[A_{k+1}A_{k+1}^{T}\right]^{-1}A_{k+1}\nabla f\left(x_{k+1}\right)$ ({}``kleinste Quadrate Methode'') \ENDIF \STATE $\Delta x_{k}:=x_{k+1}-x_{k}=\alpha_{k}p_{k}$ \STATE $y_{k}:=\nabla_{x}L\left(x_{k+1},\lambda_{k+1}\right)-\nabla_{x}L\left(x_{k},\lambda_{k}\right)$ \STATE Berechne $B_{k+1}$ aus $B_{k}$, z.B. mittels BFGS-Aufdatierung \ENDFOR \end{algorithmic} \end{algorithm} \section{Trust-Region SQP-Verfahren} \index{SQP!Trust-Region}Betrachte das Problem\begin{eqnarray*} & \min & \nabla f\left(x_{k}\right)p+\frac{1}{2}p^{T}B_{k}p\\ & \textrm{s.d.} & A_{k}p+c\left(x_{k}\right)=0\\ & & \left\Vert p\right\Vert \leq\Delta_{k}\end{eqnarray*} \section{Alternativer Zugang zur Beschreibung von SQP-Verfahren} \begin{eqnarray*} & \min & f\left(x\right)\\ & \textrm{s.d.} & c\left(x\right)\leq0\\ & & \textrm{MFCQ erfüllt}\end{eqnarray*} Das KKT-System lautet:\begin{eqnarray*} \nabla f\left(x\right)+y^{T}A\left(x\right) & = & 0\\ y & \geq & 0\\ c\left(x\right) & \leq & 0\\ \left\langle y,c\left(x\right)\right\rangle & = & 0.\end{eqnarray*} Das (gestörte) Kojima-System\index{Kojima-System} ist \begin{eqnarray*} \nabla f\left(x\right)+\left(y^{+}\right)^{T}A\left(x\right) & = & 0\\ c\left(x\right) & = & y^{-}+\varepsilon y^{+}\end{eqnarray*} mit $y_{i}^{+}=\max\left\{ 0,y_{i}\right\} $, $y_{i}^{-}=\min\left\{ 0,y_{i}\right\} $, $\varepsilon y^{+}=\left(\varepsilon_{1}y_{1}^{+},\ldots,\varepsilon_{m}y_{m}^{+}\right)^{T}$ für ein kleines $\varepsilon>0$. Man wende ein {}``verallgemeinertes'' (d.h. mit verallgemeinerter Ableitung) Newton-Verfahren auf das gestörte Kojima-System an. Durch die Wahl von $\varepsilon$ bekommt man eine Menge von Verfahren, die {}``äquivalent'' zur Menge der SQP-Verfahren ist. (\cite{KuKl}) Zu einem KKT-Punkt $\left(x,y\right)$ gehört der Punkt $\left(x,y+c\left(x\right)\right)$ im Kojima-System, zu einem Kojima-Punkt $\left(x,y\right)$ der KKT-Punkt $\left(x,y^{+}\right)$. \chapter{Erweiterte Lagrange Funktionen und Methoden (Augmented Lagrangian)} \index{Lagrange Funktion!augmented}Literatur: \cite{Be} \section{Quadratische Straffunktion für Gleichheitsrestriktionen} \begin{eqnarray*} & \min & f\left(x\right)\\ & \textrm{s.d}. & h\left(x\right)=0\in\mathbb{R}^{m}\end{eqnarray*} Für beliebigen Strafparameter $\mu>0$ betrachte das Problem\index{Straffunktion, quadratisch}\[ \min_{x\in\mathbb{R}^{n}}Q\left(x;\mu\right)\equiv f\left(x\right)+\frac{1}{2}\mu\left\Vert h\left(x\right)\right\Vert _{2}^{2}\] \textbf{Standardvoraussetzung:} \[ f,h\in C^{1,1}\left(\mathbb{R}^{n}\right)\qquad\lim_{\left\Vert x\right\Vert \rightarrow\infty}f\left(x\right)=\infty.\] \begin{thm} Die Standardvoraussetzung impliziert für jede Folge $\mu_{k}\rightarrow\infty$ die Existenz einer Folge von globalen Minima \[ x_{k}\equiv\argmin\left\{ Q\left(x,\mu_{k}\right)\right\} \] Diese hat mindestens einen Häufungspunkt $x_{*}\in\mathbb{R}^{n}$ und es gilt\[ x_{*}\in\argmin\left\{ f\left(x\right)\left|x\in\mathbb{R}^{n},h\left(x\right)=0\right.\right\} .\] Ist an $x_{*}$ die LICQ Bedingung $\rank\left(\nabla h\left(x_{*}\right)\right)=m$ erfüllt und sind $\lambda_{*}\in\mathbb{R}^{m}$ die entsprechenden Lagrange-Multiplikatoren, so gilt\[ \lim_{k\rightarrow\infty}\mu_{k}h\left(x_{*}\right)=\lambda_{*}\in\mathbb{R}^{m}.\] Wobei o.B.d.A. angenommen wurde, daß nur ein Häufungspunkt existiert. \end{thm} \begin{proof} Wegen $Q\left(x,\mu\right)\geq f\left(x\right)$ gilt $\lim_{\left\Vert x\right\Vert \rightarrow\infty}Q\left(x,\mu\right)=\infty$, es existieren also immer globale Minima $x_{k}\in\argmin\left\{ Q\left(x,\mu_{k}\right)\right\} $. Die Folge ist beschränkt, hat also einen Häufungspunkt $x_{*}$ (o.B.d.A. eindeutig).\[ f_{*}=\inf_{x\in\mathbb{R}^{n}}f\left(x\right)\leq f\left(x_{k}\right)\leq Q\left(x_{k},\mu_{k}\right)\leq Q\left(\tilde{x}_{*},\mu_{k}\right)=f\left(\tilde{x}_{*}\right)<\infty,\] wobei $\tilde{x}_{*}$ ein globaler Minimalpunkt von $f$ ist. Daraus folgt im Grenzwert\[ \lim_{k\rightarrow\infty}\frac{1}{\mu_{k}}\left[Q\left(x_{k},\mu_{k}\right)-f\left(x_{k}\right)\right]\leq\lim_{k\rightarrow\infty}\frac{1}{\mu_{k}}\left(f\left(\tilde{x}_{*}\right)-f_{*}\right)=0\] Wegen \[ \lim_{k\rightarrow\infty}\frac{1}{\mu_{k}}\left[Q\left(x_{k},\mu_{k}\right)-f\left(x_{k}\right)\right]=\lim_{k\rightarrow\infty}\frac{1}{2}\left\Vert h\left(x_{k}\right)\right\Vert ^{2}=\frac{1}{2}\left\Vert h\left(x_{*}\right)\right\Vert ^{2}\] folgt \[ h\left(x_{*}\right)=0.\] Also ist der Häufungspunkt zulässig und wegen\[ \lim_{k\rightarrow\infty}f\left(x_{k}\right)=f\left(x_{*}\right)\leq f\left(\tilde{x}_{*}\right)\] auch global optimal. Zum Beweis der letzten Aussage nutze die Optimalitätsbedingung\[ \nabla_{x}Q\left(x_{k},\mu_{k}\right)=\nabla f\left(x_{k}\right)+\mu_{k}\nabla h\left(x_{k}\right)^{T}h\left(x_{k}\right)=0.\] Wegen LICQ hat $A_{k}:=\nabla h\left(x_{k}\right)\rightarrow\nabla h\left(x_{*}\right)$ für alle großen $k$ vollen Rank hat, so daß $A_{k}A_{k}^{T}\in\mathbb{R}^{n\times n}$ regulär ist.\begin{eqnarray*} A_{k}\nabla f\left(x_{k}\right)+A_{k}A_{k}^{T}\mu_{k}h\left(x_{k}\right) & = & 0\\ \Rightarrow\,\mu_{k}h\left(x_{k}\right) & = & -\left(A_{k}A_{k}^{T}\right)^{-1}A_{k}\nabla f\left(x_{k}\right)\\ \xrightarrow{k\rightarrow\infty}\quad\lambda_{*} & = & \left(A_{*}A_{*}^{T}\right)^{-1}A_{*}\nabla f\left(x_{*}\right)\end{eqnarray*} ist eindeutige Lsg. von KKT-Bedingung $\nabla f\left(x_{*}\right)+A\left(x_{*}\right)^{T}\lambda_{*}=0$ \end{proof} \textbf{Vorteil}: Einfache Anwendung eines unbeschränkten Optimierungsverfahrens auf ein beschränktes Problem. Die $Q\left(x,\mu\right)$ sind im Gegensatz zu $L_{1}$-Straffunktionen so oft diff.bar wie $f$ und $h$ selbst. \textbf{Nachteil}: Mit wachsenden $\mu$ werden die unbeschränkten Probleme $Q\left(x,\mu\right)$ immer schwerer lösbar, insbesondere geht die Konditionszahl der Hessematrix\[ \nabla_{x}^{2}Q\left(x_{*};\mu\right)=\nabla^{2}f\left(x_{*}\right)+\mu\nabla h\left(x_{*}\right)^{T}\nabla h\left(x_{*}\right)+\sum_{i=1}^{n}\lambda_{*}^{\left(i\right)}\nabla^{2}h\left(x_{*}\right)\] mit $\mu$ gegen unendlich. Es sei denn, es gilt $m=n$, was im allgemeinen nicht von Interesse ist. \begin{remark} Der Satz gilt im {}``wesentlichen'' schon dann, wenn die Teilprobleme so genau gelöst werden, daß\[ \lim_{k\rightarrow\infty}\left\Vert \nabla_{x}Q\left(x_{k},\mu_{k}\right)\right\Vert \equiv\lim_{k\rightarrow\infty}\left\Vert \nabla f\left(x_{k}\right)+\mu\nabla h\left(x_{k}\right)^{T}h\left(x_{k}\right)\right\Vert =0.\] \end{remark} \section{Lagrange Multiplikatoren Methode (1. Ordnung)} \textbf{Idee}: Vermeide $\mu_{k}\rightarrow\infty$ durch Addition des Lagrange-Termes $\lambda^{T}h\left(x\right)$. Dies ergibt die erweiterte Lagrange Funktion\index{Lagrange Funktion!erweitert}\[ E\left(x,\lambda,\mu\right)\equiv f\left(x\right)+\lambda^{T}h\left(x\right)+\frac{1}{2}\mu\left\Vert h\left(x\right)\right\Vert ^{2}.\] \begin{thm} \label{thm:ErwLagrangeFunktion}Voraussetzung: $f,h\in C^{2,1}\left(\mathbb{R}^{n}\right)$ und $x_{*}\in h^{-1}\left(0\right)$ ist lokales Minimum an dem die LICQ und die hinreichenden lokalen Optimalitätsbedingungen erster und zweiter Ordnung erfüllt sind. Behauptung: Für alle hinreichend großen $\mu$ ist $x_{*}$ dann unbeschränktes lokales Minimum von $E\left(x,\lambda_{*},\mu\right)$, wobei $\lambda_{*}\in\mathbb{R}^{m}$ die exakten Lagrangemultiplikatoren an $x_{*}$ sind. \end{thm} \begin{proof} Optimalitätsbedingung 1. Ordnung: \[ \nabla_{x}E\left(x_{*},\lambda_{*},\mu\right)=\nabla f\left(x_{*}\right)+\nabla h\left(x_{*}\right)^{T}\lambda_{*}+\mu\underbrace{\nabla h\left(x_{*}\right)^{T}h\left(x_{*}\right)}_{=0}=\nabla_{x}L\left(x_{*},\lambda_{*}\right)=0\] Optimalitätsbedingung 2. Ordnung: Es ist zu zeigen, daß für hinreichend große $\mu$\[ \nabla_{x}^{2}E\left(x_{*},\lambda_{*},\mu\right)=\nabla_{x}^{2}L\left(x_{*},\lambda_{*}\right)+\mu\nabla h\left(x_{*}\right)^{T}\nabla h\left(x_{*}\right)\succ0\] gilt. Mit $A=\nabla h\left(x_{*}\right)$ und $Z\in\mathbb{R}^{\left(n-m\right)\times n}$, so daß $AZ^{T}=0$ und $ZZ^{T}=I\in\mathbb{R}^{\left(n-m\right)\times\left(n-m\right)}$ läßt sich jeder Vektor $u$ in $\mathbb{R}^{n}$ darstellen als \[ u=Z^{T}v+A^{T}w.\] Multiplikation von $\nabla_{x}^{2}E\left(x_{*},\lambda_{*},\mu\right)$ von rechts und links mit $u$ ergibt\begin{eqnarray*} u^{T}\nabla_{x}^{2}E\left(x_{*},\lambda_{*},\mu\right)u & = & \left(v^{T}Z+w^{T}A\right)\left(\nabla_{x}^{2}L\left(x_{*},\lambda_{*}\right)+\mu A^{T}A\right)\left(Z^{T}v+A^{T}w\right)\\ & = & v^{T}Z\nabla_{x}^{2}L\left(x_{*},\lambda_{*}\right)Z^{T}v+\mu w^{T}\left(AA^{T}\right)^{2}w+2w^{T}A\nabla_{x}^{2}L\left(x_{*},\lambda_{*}\right)Z^{T}v\\ & & +w^{T}A\nabla_{x}^{2}L\left(x_{*},\lambda_{*}\right)A^{T}w\\ & \geq & \underbrace{\alpha\left\Vert v\right\Vert ^{2}}_{\substack{\textrm{mit $\alpha>0$ wegen SSC=}\\ \textrm{2. Order Sufficiency Condition}} }+\underbrace{\beta\mu\left\Vert w\right\Vert ^{2}}_{\substack{\textrm{mit $\beta>0$}\\ \textrm{wegen LICQ}} }-\underbrace{2\gamma\left\Vert w\right\Vert \left\Vert v\right\Vert }_{\textrm{mit }\gamma<\infty}-\underbrace{\delta\left\Vert w\right\Vert ^{2}}_{\textrm{mit }\delta<\infty}\\ & = & \alpha\left\Vert v\right\Vert ^{2}-2\gamma\left\Vert w\right\Vert \left\Vert v\right\Vert +\left(\beta\mu-\delta\right)\left\Vert w\right\Vert ^{2}\\ & = & \left(\left\Vert v\right\Vert ,\left\Vert w\right\Vert \right)^{T}\left(\begin{array}{cc} \alpha & -\gamma\\ -\gamma & \beta\mu-\delta\end{array}\right)\left(\begin{array}{c} \left\Vert v\right\Vert \\ \left\Vert w\right\Vert \end{array}\right).\end{eqnarray*} Die Matrix hat $\alpha$ als positives diagonales Element, ist also positiv definit, gdw. ihre Determinante $\alpha\left(\beta\mu-\delta\right)-\gamma^{2}$ positiv ist. Dies ist sicherlich für hinreichend große $\mu$ der Fall. \end{proof} \begin{remark} Entsprechend hat $E\left(x,\lambda,\mu\right)$ auch für alle $\lambda\approx\lambda_{*}$ und $\mu$ hinreichend groß ein lokales Minimum nahe $x_{*}$ (Im allgemeinen nicht in $h^{-1}\left(0\right)$, d.h. nicht zulässig bzgl. $h\left(x\right)=0$). \end{remark} \textbf{Frage:} Wie kann man die $\lambda_{k}$ justieren, so daß im Grenzwert $\lambda_{k}\rightarrow\lambda_{*}$ gilt? \textbf{Antwort:} Multiplikatorenregel\index{Multiplikatorenregel} 1. Ordnung: Wegen\[ x_{k+1}=\argmin\left\{ E\left(x,\lambda_{k},\mu_{k}\right)\right\} \] ist\begin{eqnarray*} 0 & = & \nabla f\left(x_{k+1}\right)+\nabla h\left(x_{k+1}\right)^{T}\lambda_{k}+\nabla h\left(x_{k+1}\right)^{T}h\left(x_{k+1}\right)\mu_{k}\\ & = & \nabla f\left(x_{k+1}\right)+\nabla h\left(x_{k+1}\right)^{T}\underbrace{\left(\lambda_{k}+h\left(x_{k+1}\right)\mu_{k}\right)}_{\equiv\lambda_{k+1}}.\end{eqnarray*} Dies ergibt die (natürliche) Aufdatierung\[ \lambda_{k+1}=\lambda_{k}+\mu_{k}h\left(x_{k+1}\right).\] Außerdem sollte $\mu_{k+1}\geq\mu_{k}$ möglichst so gewählt werden, daß $\nabla_{x}^{2}E\left(x_{k+1},\lambda_{k+1},\mu_{k+1}\right)$ positiv definit ist. \begin{thm} Unter den Voraussetzungen von Satz \ref{thm:ErwLagrangeFunktion} folgt aus $\mu$ hinreichend groß und $\left\Vert x_{0}-x_{*}\right\Vert $, $\left\Vert \lambda_{0}-\lambda_{*}\right\Vert $ hinreichend klein die lineare Konvergenz der Iteration\begin{eqnarray*} x_{k+1} & = & \argmin_{x}\left\{ E\left(x,\lambda_{k},\mu_{k}\right)\right\} ,\\ \lambda_{k+1} & = & \lambda_{k}+h\left(x_{k}\right)\mu\end{eqnarray*} zum beschränkten Minimum $x_{*}$ mit Multiplikatoren $\lambda_{*}$. \end{thm} \begin{proof} Benutze Standardergebnis für Iteration $z=G\left(y\right)$ mit $y_{*}=G\left(y_{*}\right)$ und $G\in C^{1,1}\left(\mathbb{R}^{n}\right)$: Falls alle Eigenwerte von $G_{y}\left|_{y_{*}}\right.\equiv\frac{dz}{dy}\left|_{y_{*}}\right.\in\mathbb{R}^{\left(n+m\right)\times\left(n+m\right)}$ betragsmäßig kleiner eins sind, folgt aus $\left\Vert y_{0}-y_{*}\right\Vert \approx0$, daß $y_{k+1}=G\left(y_{k}\right)\xrightarrow{k\rightarrow\infty}y_{*}$. Hier haben wir eine implizite Beziehnung zwischen $y=\left(x_{k},\lambda_{k}\right)$ und $z=\left(x_{k+1},\lambda_{k+1}\right)$, nämlich\[ 0\equiv F\left(y,z\right)=F\left(\left(x_{k},\lambda_{k}\right),\left(x_{k+1},\lambda_{k+1}\right)\right):=\left[\begin{array}{c} \nabla_{x}L\left(x_{k+1},\lambda_{k}\right)+\mu\nabla h\left(x_{k+1}\right)^{T}h\left(x_{k+1}\right)\\ \lambda_{k+1}-\lambda_{k}-\mu h\left(x_{k+1}\right)\end{array}\right].\] Daraus folgt\begin{eqnarray*} F_{y}=\frac{\partial F}{\partial\left(x_{k},\lambda_{k}\right)} & = & \left[\begin{array}{cc} 0 & \nabla h\left(x_{k+1}\right)^{T}\\ 0 & -I\end{array}\right]\\ F_{z}=\frac{\partial F}{\partial\left(x_{k+1},\lambda_{k+1}\right)} & = & \left[\begin{array}{cc} \nabla_{x}^{2}E\left(x_{k+1},\lambda_{k+1}\right) & 0\\ -\mu\nabla h\left(x_{k+1}\right) & I\end{array}\right].\end{eqnarray*} $\rho$ ist Eigenwert von $\frac{dz}{dy}=F_{z}^{-1}F_{y}$, gdw. \begin{eqnarray*} & & 0=\det\left[\rho I+F_{z}^{-1}F_{y}\right]\\ & \Leftrightarrow & 0=\det\left(\rho F_{z}+F_{y}\right)\qquad\textrm{da }\det F_{z}\neq0\\ & \Leftrightarrow & \rho=0\textrm{ oder }P\left(\rho\right):=\det\left(\left[\begin{array}{cc} \nabla_{x}^{2}E\left(x_{k+1},\lambda_{k+1}\right) & \frac{1}{\rho}\nabla h^{T}\left(x_{k+1}\right)\\ -\mu\nabla h\left(x_{k+1}\right) & \left(1-\frac{1}{\rho}\right)I\end{array}\right]\right)=0.\end{eqnarray*} Zu zeigen: Alle solche $\rho$ sind reell und liegen für hinreichend großes $\mu>0$ im Inneren des Intervalles $\left(-1,1\right)$. Ab jetzt werden alle Ableitungen direkt an der Lösung $\left(x_{*},\lambda_{*}\right)$ ausgewertet. Die hergeleiteten Ungleichheitsaussagen gelten dann aufgrund der Stetigkeit in einer Umgebung. Fall $\rho=1$: Dann ist\[ P\left(1\right)=\det\left(\left[\begin{array}{cc} \nabla_{x}^{2}L+\mu\nabla h^{T}\nabla h & \nabla h^{T}\\ \nabla h & 0\end{array}\right]\right)\left(-\mu\right)^{m}.\] Wegen der Optimalitätsbedingung 2. Ordnung gilt für alle $Z$ mit $\nabla h\, Z=0$\[ Z\nabla_{x}^{2}LZ^{T}=Z\left(\nabla_{x}^{2}L+\mu\nabla h^{T}\nabla h\right)Z^{T}\succ0.\] Daraus folgt wiederum $P\left(1\right)\neq0$. Fall $\rho\neq1$: Dann ist\begin{eqnarray*} P\left(\rho\right)\frac{\rho^{m}}{\left(\rho-1\right)^{m}} & = & \det\left(\left[\begin{array}{cc} \nabla_{x}^{2}E & \frac{1}{\rho-1}\nabla h^{T}\\ -\mu\nabla h & I\end{array}\right]\right)\\ & = & \det\left(\left[\begin{array}{cc} \nabla_{x}^{2}E+\frac{\mu}{\rho-1}\nabla h^{T}\nabla h & 0\\ 0 & I\end{array}\right]\right)\\ & = & \det\left(\nabla_{x}^{2}L+\left(\mu+\frac{\mu}{\rho-1}\right)\nabla h^{T}\nabla h\right)\\ & = & \det\left(\nabla_{x}^{2}L+\frac{\mu\rho}{\rho-1}\nabla h^{T}\nabla h\right).\end{eqnarray*} Da die Matrix symmetrisch ist, sind alle Wurzeln $\rho$ reell. Außerdem wurde im Beweis von Satz \ref{thm:ErwLagrangeFunktion} gezeigt, daß die Matrix positiv definit ist, falls $\frac{\mu\rho}{\rho-1}\geq\bar{\mu}$ für eine untere Schranke $\bar{\mu}$. Für $\mu\geq\bar{\mu}$ und $\rho>1$ oder $\mu\geq2\bar{\mu}$ und $\rho\leq-1$ ist dies erfüllt. \end{proof} \section{Verallgemeinerung auf Ungleichungen} $g_{j}\left(x\right)\leq0$ gdw. $g_{j}\left(x\right)+s_{j}^{2}=0$ für eine Schlupfvariable $s_{j}^{2}$. Dies ergibt in der erweiterten Lagrange-Funktion den Term\[ E_{j}=\lambda_{j}\left(g_{j}\left(x\right)+s_{j}^{2}\right)+\frac{\mu}{2}\left(g_{j}\left(x\right)+s_{j}^{2}\right)^{2},\] wobei dafür Sorge getragen wird, daß immer $\lambda_{j}\geq0$ gilt. $E_{j}$ hat bzgl. $g_{j}\left(x\right)+s_{j}^{2}$ das unbeschränkte Minimum\[ g_{j}\left(x\right)+s_{j}^{2}=-\frac{\lambda_{j}}{\mu}.\] Dies kann nur durch die Wahl von $s_{j}$ erreicht werden, falls $g_{j}\left(x\right)\leq-\frac{\lambda_{j}}{\mu}$. Dann ist $E_{j}\equiv-\frac{\lambda_{j}^{2}}{2\mu}$. Ist $g_{j}\left(x\right)>-\frac{\lambda_{j}}{\mu}$, so ist $E_{j}=-\lambda_{j}g_{j}\left(x\right)+\frac{\mu}{2}g_{j}^{2}\left(x\right)$. Für $\mu\rightarrow\infty$ erhält man $E_{j}\rightarrow\begin{cases} 0 & \textrm{für }g_{j}\leq0\\ \infty & \textrm{für }g_{j}>0\end{cases}$. $E_{j}$ wird auch {}``\emph{shifted quadratic penalty function}''\index{shifted quadratic penalty function} genannt. $E_{j}$ ist bzgl. $g_{j}$und damit $x$ zwar stetig differenzierbar, hat aber am {}``Schaltpunkt'' $g_{j}=-\frac{\lambda_{j}}{\mu}$ einen Sprung in der zweiten Ableitung. \begin{remark} Während die Slackvariable $s_{j}^{2}$ explizit eliminiert wurde, bleibt der Lagrange Multiplikator $\lambda_{j}$ erhalten. Die Aufdatierungsformel 1. Ordnung ist einfach\[ \lambda_{j}^{\left(n+1\right)}=\max\left\{ 0,\lambda_{j}^{\left(n\right)}+\mu^{\left(n\right)}g_{j}\left(x^{\left(n\right)}\right)\right\} .\] \end{remark} \appendix \newpage \chapter{Interior Point Methods} Vorlesung von Prof. Florian A. Potra, University of Maryland Baltimore Country. Dantzig, 1947: Simplex Algorithm\index{LP} \[ \left(P\right)\qquad\begin{array}{rl} \min & c^{T}x\\ \textrm{s.t.} & Ax=b\\ & x\geq0\end{array}\qquad\qquad\left(D\right)\begin{array}{rl} \max & b^{T}y\\ \textrm{s.t.} & A^{T}y+s=c\\ & s\geq0\end{array}\] \begin{eqnarray*} \mathcal{P} & := & \left\{ x\in\mathbb{R}^{n}\left|Ax=b,x\geq0\right.\right\} \\ \mathcal{D} & := & \left\{ \left(s,y\right)\in\mathbb{R}^{m}\left|A^{T}y+s=c,s\geq0\right.\right\} \end{eqnarray*} $x$ primal variables, $y$ dual variables, $s$ slack variables We observe, that\[ c^{T}x-b^{T}y=x^{T}c-y^{T}b=x^{T}\left(A^{T}y+s\right)-y^{T}Ax=x^{T}A^{T}y+x^{T}s-y^{T}Ax=x^{T}s\geq0\] So, for all $x\in\mathcal{P}$ and all $\left(s,y\right)\in\mathcal{D}$\index{duality gap}\[ 0\leq c^{T}x+b^{T}y=x^{T}s=:\textrm{duality gap}\] \begin{thm} $\exists x^{*}\in\mathcal{P}$, $\left(s^{*},y^{*}\right)\in\mathcal{D}$, s.t. $x^{*T}s^{*}=0$ $\Rightarrow$ $c^{T}x^{*}=\min_{x\in\mathcal{P}}c^{T}x$, $b^{T}y^{*}=\max_{\left(y,s\right)\in\mathcal{D}}b^{T}y$ \end{thm} \begin{proof} ~\begin{eqnarray*} c^{T}x & \geq & b^{T}y\qquad\forall x\in\mathcal{P},\left(s,y\right)\in\mathcal{D}\end{eqnarray*} So, \[ c^{T}x\geq b^{T}y^{*}=c^{T}x^{*}\qquad\forall x\in\mathcal{P}\] on the one hand, and\[ b^{T}y\leq c^{T}x^{*}=b^{T}y^{*}\qquad\forall\left(s,y\right)\in\mathcal{D}.\qedhere\] \end{proof} And we see,\begin{eqnarray*} 0 & \leq & c^{T}x-c^{T}x^{*}\leq x^{T}s\qquad x\geq0\\ 0 & \leq & b^{T}y^{*}-b^{T}y\leq x^{T}s\qquad s\geq0\end{eqnarray*} $x^{T}s=0$ $\Leftrightarrow$ $x_{i}s_{i}=0$ $\forall i\in\left\{ 1,\ldots,n\right\} $. (complementarity) \[ \left(PD\right)\qquad\begin{array}{rcl} xs & = & 0\\ Ax-b & = & 0\\ A^{T}y+s-c & = & 0\\ x & \geq & 0\\ s & \geq & 0\end{array}\] Denote by\begin{eqnarray*} z & = & \left(\begin{array}{c} x\\ s\\ y\end{array}\right)\\ F\left(z\right) & = & \left(\begin{array}{c} xs\\ Ax-b\\ A^{T}y+s-c\end{array}\right)\qquad F:\mathbb{R}^{2n+m}\rightarrow\mathbb{R}^{2n+m}.\end{eqnarray*} We need to find a point $z$ with $F\left(z\right)=0$. Using Newton's Method, we obtain \begin{eqnarray*} z^{+} & = & z-F'\left(z\right)^{-1}F\left(z\right)\\ F'\left(z\right) & = & \left(\begin{array}{ccc} S & X & 0\\ A & 0 & 0\\ 0 & I & A^{T}\end{array}\right)\qquad S=\left(\begin{array}{ccc} s_{1} & & 0\\ & \ddots\\ 0 & & s_{n}\end{array}\right)\end{eqnarray*} Let $A$ have full rank. \textbf{Annahme:} $\exists x^{0}\in\mathcal{P},\left(s^{0},y^{0}\right)\in\mathcal{D}$, $x^{0}>0$, $s^{0}>0$ (\emph{interior point condition})\index{interior point condition} Now, denote by\begin{eqnarray*} F_{\tau}\left(z\right) & := & \left(\begin{array}{c} xs-\tau e\\ Ax-b\\ A^{T}y+s-c\end{array}\right),\qquad e=\left(\begin{array}{c} 1\\ \vdots\\ 1\end{array}\right)\in\mathbb{R}^{n},\end{eqnarray*} so $F_{0}\left(z\right)=F\left(z\right)$. The idea of path following interior point method\index{interior point method, path following} is: Given a point \[ z\in\mathcal{F}^{0}=\left\{ \left.z=\left(\begin{array}{c} x\\ s\\ y\end{array}\right)\right|Ax=b,A^{T}y+s=c,x>0,y>0\right\} \] compute $\mu=\frac{x^{T}s}{n}$ and apply a newton step to $F_{\bar{\mu}}\left(t\right)=0$ where $\bar{\mu}=\left(1-\gamma\right)\mu$ for some $\gamma\in\left(0,1\right)$.\begin{eqnarray} & & \Delta z=\left(\begin{array}{c} u\\ v\\ w\end{array}\right)=-F_{\tau}'\left(z\right)^{-1}F_{\tau}\left(z\right)\nonumber \\ \Leftrightarrow & & F_{\tau}'\left(z\right)\Delta z=-F_{\tau}\left(z\right)\nonumber \\ \Leftrightarrow & & \begin{array}{rcl} Su+Xv & = & \tau e-Xs\\ Au & = & b-Ax=0\\ v+A^{T}w & = & c-s-A^{T}y=0\end{array}\label{eq:IPnewtonstep}\end{eqnarray} \begin{lem} If $z\in\mathcal{F}^{0}$ then the solution of the system (\ref{eq:IPnewtonstep}) satisfies \begin{enumerate} \item $u^{T}v=0$ \item $\left\Vert Du\right\Vert ^{2}+\left\Vert D^{-1}v\right\Vert ^{2}=\left\Vert \left(XS\right)^{-\frac{1}{2}}\left(\tau e-xs\right)\right\Vert ^{2}=:\left\Vert a\right\Vert ^{2}$, $\sqrt{x}=\left(\sqrt{x_{1}},\ldots,\sqrt{x_{n}}\right)^{T}$ \item $\left\Vert uv\right\Vert \leq\frac{1}{\sqrt{8}}\left\Vert a\right\Vert ^{2}$ \item $\left\Vert uv\right\Vert _{\infty}\leq\frac{1}{4}$$\left\Vert a\right\Vert ^{2}$ \item $\left\Vert uv\right\Vert _{1}\leq\frac{1}{2}\left\Vert a\right\Vert ^{2}$ \end{enumerate} where $D=X^{-\frac{1}{2}}S^{\frac{1}{2}}$. \end{lem} \begin{proof} ~ \begin{enumerate} \item $Au=0$ $\Leftrightarrow$ $u\in\ker\left(A\right)$ and $v+A^{T}w=0$ $\Leftrightarrow$ $v\in\textrm{Range}$$\left(A^{T}\right)$ So, $u^{T}v=0$, because $u^{T}v=\left\langle u,v\right\rangle =\left\langle u,-A^{T}w\right\rangle =\left\langle Au,-w\right\rangle =0$. \item $Su+Xv=\tau e-xs$ $\bar{u}=Du$, $\bar{v}=D^{-1}v$, $a=\bar{u}+\bar{v}$ $Du+D^{-1}v=a$, \[ \left\Vert a\right\Vert ^{2}=\left\Vert Du+D^{-1}v\right\Vert ^{2}=\left\Vert Du\right\Vert ^{2}+\left\Vert D^{-1}v\right\Vert ^{2}+2\left(Du\right)^{T}\left(Dv\right)^{-1}=2u^{T}v=0\qedhere\] \end{enumerate} \end{proof} \begin{thm} Assume $z\in\mathcal{F}^{0}$, $\left\Vert xs-\mu e\right\Vert \leq\alpha\mu$, $\alpha\in\left(0,1\right)$, $\mu=\frac{x^{T}s}{n}$. Then by taking $\gamma=\ldots$ and \begin{eqnarray*} \bar{\mu} & = & \tau=\left(1-\gamma\right)\mu,\\ \bar{z} & = & z-\left(F_{\tau}'\left(z\right)\right)^{-1}F\left(z\right),\end{eqnarray*} we have $\bar{z}\in\mathcal{F}^{0}$ and\[ \left\Vert \bar{x}^{T}\bar{s}-\bar{\mu}e\right\Vert \leq\alpha\bar{\mu}.\] \end{thm} \begin{proof} $z\left(\gamma\right)=\bar{z}$, \begin{eqnarray*} \bar{x}\bar{s} & = & \left(x+u\right)\left(s+v\right)=xs+su+xv+uv\\ & = & xs+\tau e-xs+uv\qquad\textrm{by (\ref{eq:IPnewtonstep})}\\ & = & \tau e+uv\\ & = & \left(1-\gamma\right)\mu e+uv\\ \Rightarrow\,\bar{\mu} & = & \frac{\bar{x}^{T}\bar{s}}{n}=\left(1-\gamma\right)\mu,\qquad\textrm{because }\sum\bar{x}_{i}\bar{s}_{i}=n\left(1-\gamma\right)\mu+\underbrace{\sum u_{i}v_{i}}_{=0}\\ \Rightarrow\,\bar{x}\bar{s}-\bar{\mu}e & = & uv\end{eqnarray*} so $\left\Vert \bar{x}\bar{s}-\bar{\mu}e\right\Vert \leq\alpha\bar{\mu}$ $\Leftrightarrow$ $\left\Vert uv\right\Vert \leq\alpha\bar{\mu}$. From the Lemma: $\left\Vert uv\right\Vert \leq\frac{1}{\sqrt{8}}\left\Vert \left(XS\right)^{-\frac{1}{2}}\left(\left(1-\gamma\right)\mu e-xs\right)\right\Vert ^{2}$, so we have to determine $\gamma$ such that $\left\Vert \left(XS\right)^{-\frac{1}{2}}\left(\left(1-\gamma\right)\mu e-xs\right)\right\Vert \leq\alpha\bar{\mu}=\alpha\left(1-\gamma\right)\mu$. We have\begin{eqnarray*} \left\Vert xs-\mu e\right\Vert ^{2} & = & \sum\left(x_{i}s_{i}-\mu\right)^{2}\leq\alpha^{2}\mu^{2}\\ \Rightarrow\,\left|x_{i}s_{i}-\mu\right| & \leq & \alpha\mu\\ \Rightarrow\, x_{i}s_{i} & \geq & \left(1-\alpha\right)\mu\\ \Rightarrow\,\left\Vert \left(XS\right)^{-\frac{1}{2}}\left(\left(1-\gamma\right)\mu e-xs\right)\right\Vert ^{2} & \leq & \frac{1}{\min x_{i}s_{i}}\left\Vert \left(1-\gamma\right)\mu e-xs\right\Vert ^{2}\\ & \leq & \frac{\left\Vert \mu e-xs\right\Vert ^{2}+\left\Vert \left(1-\gamma\right)\mu e-\mu e\right\Vert ^{2}}{\left(1-\alpha\right)\mu}\\ & \leq & \frac{1}{\left(1-\alpha\right)\mu}\left(\alpha^{2}\mu^{2}+\mu^{2}\gamma^{2}n\right)\\ & = & \frac{\alpha^{2}+\gamma^{2}n}{\left(1-\alpha\right)}\mu\\ \Rightarrow\,\left\Vert \bar{x}\bar{s}-\bar{\mu}e\right\Vert & \leq & \frac{1}{\sqrt{8}}\cdot\frac{\alpha^{2}+\gamma^{2}n}{\left(1-\alpha\right)}\mu\\ & \stackrel{?}{\leq} & \alpha\left(1-\gamma\right)\mu\end{eqnarray*} Assume the stronger condition $\alpha\leq\frac{\sqrt{8}}{1+\sqrt{8}}$. So $\alpha^{2}\left(1+\sqrt{8}\right)<\sqrt{8}\alpha$ implies $\alpha^{2}<\sqrt{8}\alpha\left(1-\alpha\right)$. Then, $\frac{\alpha^{2}+\gamma^{2}n}{\left(1-\alpha\right)\sqrt{8}}\mu\leq\alpha\left(1-\gamma\right)\mu$ is satisfied, iff\begin{eqnarray*} \alpha^{2}+\gamma^{2}n & \leq & \sqrt{8}\alpha\left(1-\alpha\right)\left(1-\gamma\right)\\ \Leftrightarrow\,\gamma^{2}n+\underbrace{\sqrt{8}\alpha\left(1-\alpha\right)}_{2\lambda}\gamma-\underbrace{\left(\sqrt{8}\alpha\left(1-\alpha\right)-\alpha^{2}\right)}_{\sigma>0} & \leq & 0\\ \gamma^{*} & = & \frac{-\lambda+\sqrt{\lambda^{2}+\sigma n}}{n}\\ & = & \frac{\sigma n}{n\left(\sqrt{\lambda^{2}+\sigma n}+\lambda\right)}\\ & \xrightarrow{n\rightarrow\infty} & \frac{\sqrt{\sigma}}{\sqrt{n}}.\end{eqnarray*} Want an elegant lower bound for $\gamma^{*}$. Set $\gamma^{*}=\frac{\delta}{\sqrt{n}}$,\[ \delta=\frac{\sigma}{\sqrt{\frac{\lambda^{2}}{n}+\sigma}+\frac{\lambda}{\sqrt{n}}}\geq\frac{\sigma}{\sqrt{\frac{\lambda^{2}}{2}+\sigma}+\frac{\lambda}{\sqrt{2}}}=:\hat{\delta}\] If we take $\hat{\gamma}=\frac{\hat{\delta}}{\sqrt{n}}$ we have $\left\Vert \bar{x}\bar{s}-\bar{\mu}e\right\Vert \leq\alpha\bar{\mu}$ provided $\alpha<\frac{\sqrt{8}}{1+\sqrt{8}}$. \end{proof} \begin{algorithm} ~ \begin{algorithmic}\REQUIRE $z^{0}\in\mathcal{F}^{0}$, $z^{0}\in N_{\alpha}$ where $N_{\alpha}=\left\{ z\left|\left\Vert xs-\mu e\right\Vert \leq\alpha\right.\right\} $ \COMPUTE \FOR{$k=0,1,\ldots$} \STATE $z\leftarrow z^{k}$ \STATE determine $\gamma$ such that $\bar{z}\in N_{\alpha}$ \STATE $z^{k+1}\leftarrow\bar{z}$, $\gamma_{k}\leftarrow\gamma$ \ENDFOR \end{algorithmic} Properties: $z^{k}\in N_{\alpha}$, $\mu_{k+1}\leq\left(1-\gamma_{k}\right)\mu_{k}$, $k=01,2,\ldots$ \end{algorithm} \begin{cor} $\mu_{k}\leq\left(1-\frac{\hat{\delta}}{\sqrt{n}}\right)^{k}\mu_{0}$ \end{cor} $\frac{\mu_{k}}{\mu_{0}}\leq\varepsilon$ for $\left(1-\frac{\hat{\delta}}{\sqrt{n}}\right)^{k}\leq\varepsilon$, i.e. $k\log\left(1-\frac{\hat{\delta}}{\sqrt{n}}\right)\leq\log\varepsilon$. \begin{remark} If $A,b,c$ are integers and $L$ ist the binary length of data and $\varepsilon<2^{-L}$ then $z$ can be rounded to an exact solution of LP in $O\left(n^{3}\right)$ operations. An exact solution of (PD) is obtained in $O\left(\sqrt{n}L\right)$ iterations + $O\left(n^{3}\right)$ arithmetic operations. \end{remark} Smart algorithms determine $\gamma_{k}$ such that $\gamma_{k}\nearrow1$. From $su+xv=\left(1-\gamma\right)\mu e-xs$ obtain the \begin{itemize} \item affine scaling direction\index{scaling direction!affine} from \begin{equation} su^{a}+xv^{a}=-xs\label{eq:affscalingdir}\end{equation} Newtonstep: $-F'\left(z\right)^{-1}F_{0}\left(z\right)$ \item centering scaling direction\index{scaling direction!centering} from \begin{equation} su^{c}+xv^{c}=\mu e-xs\label{eq:centscalingdir}\end{equation} Newtonstep: $-F'\left(z\right)^{-1}F_{\mu}\left(z\right)$ \end{itemize} With $u=\gamma u^{a}+\left(1-\gamma\right)u^{c}$, $v=\gamma v^{a}+\left(1-\gamma\right)v^{c}$, it holds\begin{eqnarray*} su+xv & = & \gamma\left(su^{a}+xv^{a}\right)+\left(1-\gamma\right)\left(su^{c}+xv^{c}\right)\\ & = & \gamma\left(-xs\right)+\left(1-\gamma\right)\left(\mu e-xs\right)\\ & = & \left(1-\gamma\right)\mu e-xs\end{eqnarray*} $\left\Vert uv\right\Vert \leq\alpha\left(1-\gamma\right)\mu$ holds, if\begin{equation} \left\Vert \gamma^{2}u^{a}v^{a}+\gamma\left(1-\gamma\right)\left(u^{a}v^{c}+u^{c}v^{a}\right)+\left(1-\gamma\right)^{2}u^{c}v^{c}\right\Vert \leq\alpha\left(1-\gamma\right)\mu\label{eq:affcentscalingdirgamma}\end{equation} Solve (\ref{eq:affcentscalingdirgamma}) for $\gamma$. \begin{algorithm} [Gonzaga, Largest Step]Solve (\ref{eq:affscalingdir}) and (\ref{eq:centscalingdir}) and then get the longest $\gamma$, that satisfies (\ref{eq:affcentscalingdirgamma}). \end{algorithm} \textbf{Note}: (\ref{eq:affscalingdir}) and (\ref{eq:centscalingdir}) have the same Jacobian. \begin{eqnarray*} N_{\alpha} & = & \left\{ z\left|\left\Vert \frac{xs}{\mu}-e\right\Vert _{2}\leq\alpha\right.\right\} \\ N_{\alpha}^{-} & = & \left\{ z\left|\left\Vert \left(\frac{xs}{\mu}-e\right)^{-}\right\Vert _{\infty}\leq\alpha\right.\right\} \\ \mathcal{D}_{\beta} & = & N_{\left(1-\beta\right)}^{-}=\left\{ z\left|xs\geq\beta\mu e\right.\right\} \\ \lim_{\alpha\nearrow1}N_{\alpha}^{-} & = & \mathcal{F}\end{eqnarray*} \begin{thebibliography}{NoWr} \bibitem[JaSt]{JaSt}Jarre, Stoer, Optimierung (gründlich) \bibitem[NoWr]{NoWr}Nocedal, Wright, Numerical Optimization (lesbar) \bibitem[KiSc]{KiSc}Kielbasinski, Schwetlick, Numerische lineare Algebra \bibitem[KuKl]{KuKl}Kummer, Klatte, Nonsmooth Equations in Optimization \bibitem[Be]{Be}D. Bertsehas, Constraint Optimization and Lagrange Multiplier Methods, 1982 \end{thebibliography} \printindex{} \end{document}