11 Implementierungen, Anwendungsbeispiele, Backpropagation

11.1 Exkurs – Gradienten und Repräsentation

Die Berechnung von Gradienten ist ebenso wesentlich wie schwierig in numerischen Algorithmen. Viele erfolgreiche Algorithmen basieren auf effizient berechenbaren Darstellungen oder Approximationen des Gradienten.

Ein erstes Beispiel ist der stochastische Gradientenabstieg, der auf einen Schätzer statt des eigentlichen Gradienten baut.

In der Optimierung mit (partiellen) Differentialgleichungen wie \[\begin{equation*} \mathcal J(u) = \frac 12 \int_0^T\|x(s)-x^*\|^2 + \|u(s)\|^2\,\mathsf{d}s\to \min_{u} \quad \text{s.t. }\dot x(t) = Ax(t) + Bu(t) \end{equation*}\] macht es einen entscheidenden Unterschied, dass der Gradient \[\partial_u \mathcal J(u) = u + B^Tp \] (Jan beachte, dass \(u\) eine Funktion ist, also ein \(\infty\)-dimensionales Objekt) über die Lösung der adjungierten Gleichung \[\begin{equation*} -\dot p(t) = A^T p(t) + x(t)-x^*, \quad p(T) = 0 \end{equation*}\] definiert ist (siehe beispielsweise (Kurdila and Zabarankin 2005, Theorem 5.5.1) für ein abstraktes Resultat und bspw. (Tröltzsch 2009, Abschnitt 5.9.1) für eine Umsetzung in einem Gradientenabstiegsverfahren mit PDEs).

Während diese Resultate in der Theorie die Wohlgestelltheit der Optimierungsprobleme helfen zu analysieren, werden Sie in der Praxis gerne verwendet weil die Lösung einer Differentialgleichung einfacher umzusetzen ist (und im Zweifel auch effizienter) als die Berechnung von sehr abstrakten Gradienten.

Um abstraktere Gradienten zu charakterisieren hilft oftmals die Definition der totalen Ableitung einer Funktion \(f\colon X\to Y\) bei einem \(x\in X\) als die lineare Abbildung \(L(x)\colon X\to Y\) derart, dass \[\begin{equation*} f(x+h) - f(x) - L(x)[h] = o(\| h\|) \end{equation*}\] für \(h\in X\) gilt.

Seien beispielsweise die Räume als \(X=Y=\mathbb R^{n\times n}\) gegeben und \[\begin{equation} f(S) = A^TS + SA -SRS + Q \tag{11.1} \end{equation}\] dann hat der Riccati Operator \(f\) mit Koeffizienten \(A\), \(R\), \(Q \in \mathbb R^{n\times n}\) demnach die Realisierung des Gradienten (an der Stelle \(S_0\)) gegeben als \[\begin{equation*} L(S_0)[h] = (A^T-S_0R)h + h(A-RS_0). \end{equation*}\]

Für neuronale Netze ist der Gradient von \(f(x; \tilde A, b)=\tilde Ax+b\) bezüglich der Gewichte \(\tilde A\in \mathbb R^{m\times n-1}\) und \(b\in \mathbb R^{m}\) interessant. Zunächst mal verstehen wir die affine lineare Abbildung im projektiven Raum vermöge \[\begin{equation*} y = Ax + b \longleftrightarrow \begin{bmatrix} y \\ 1 \end{bmatrix} = \begin{bmatrix}A & b \\ 0 & 1 \end{bmatrix}\begin{bmatrix} x \\ 1 \end{bmatrix} \end{equation*}\] und betrachten einfach \(f(x; A)=Ax\) mit \(A\in \mathbb R^{m\times n}\). Gemäß der Ableitungsformel gilt \[\begin{equation*} f(x;A+h) - f(x;A) - hx = (A+h)x - Ax - hx = 0 = o(\|h\|) \end{equation*}\] also im Prinzip ist \(L(A)\leftrightarrow x\) die Ableitung wenn auch nicht als Vektor sondern als \[\begin{equation*} L_{[x]}\colon \mathbb R^{n\times m} \to \mathbb R^{n}\colon h\to L_{[x]} h = hx. \end{equation*}\]

Weil hier an der Stelle \(A\) abgeleitet wird, sollte hier \(L(A)\) stehen. Da bei linearen Abbildung die Ableitung überall gleich ist, gilt \(L(A)\equiv L\) unabhängig von \(A\). Für später wird es interessant sein, dass die Ableitung direkt mit \(x\) zusammenhängt, sodass wir das \(x\) in die Definition mit aufnehmen.

Diese Abbildung \(L_{[x]}\) ist linear und könnte als Matrix (bzw. Tensor im Sinne einer höherdimensionalen Datenstruktur) im \(\mathbb R^{n \times (n\times m)}\) realisiert werden.

11.2 Backpropagation

Die besondere Architektur neuronaler Netze, wie sie im Machine Learning verwendet werden, impliziert eine besondere Struktur in der Gradientenberechnung.

Wir betrachten ein 2-layer Netzwerk \[\begin{equation*} s(x) = \sigma(A_2 \sigma(A_1x)) \end{equation*}\] mit den bias-Vektoren \(b_1\) und \(b_2\) wie oben erläutert bereits in den Matrizen \(A_1\in \mathbb R^{n_0\times n_1}\) und \(A_2 \in \mathbb R^{n_2\times n_1}\) integriert. Für das learning wird zu Datenpunkten \((x, y) \in \mathbb R^{n_0}\times \mathbb R^{n_2}\) die Ableitung der Funktion \[\begin{equation*} (A_1, A_2) \to l((x,y); A_1, A_2) = \|y-s(x)\|_2^2 \end{equation*}\] bezüglich des Parameter“vektors” \(p=(A_2, A_1)\) gesucht.

Eine erste Anwendung der Kettenregel ergibt \[\begin{equation*} \partial_p l((x,y); p) = -2\bigl \langle y-s(x;p), \, \partial_p s(x;p)\bigr\rangle \end{equation*}\] sodass wir uns erstmal um die Ableitung von \(s\) bezüglich \(p\) kümmern können.

Bezüglich \(A_2\) ergibt die Kettenregel direkt \[\begin{equation*} \partial_{A_2} = \sigma'(A_2x_1)L_{[x_1]} \end{equation*}\] wobei

• \(x_1:=\sigma(A_1x)\) der Zwischenwert nach der vorletzten (hier der ersten) layer ist,
• \(\sigma'(A_2x_1)\in \mathbb R^{n_2'\times n_2}\) die (diagonale) Jacobimatrix der komponentenweise definierten Aktivierungsfunktion \(\sigma\)
• und \(L_{[x_1]}\in \mathbb R^{n_2\times (n_2\times n_1)}\) die Ableitung von \(A_2x_1\) nach \(A_2\) (vgl. oben).

Bezüglich \(A_1\) ergibt die Kettenregel direkt \[\begin{equation*} \partial_{A_1} = \sigma'(A_2x_1)A_2\sigma'(A_1x)L_{[x]} \end{equation*}\]

11.2.1 Rekursion

Wir bemerken, dass die \(A_i\) Komponenten des Gradienten von \(s\) über \[\begin{equation*} \delta_N := \sigma'(A_Nx_{N-1}), \quad \partial_{A_N} s(x) = \delta_N L_{[x_{N-1}]} \end{equation*}\] für die letzte layer und dann rekursiv über

\[\begin{equation*} \delta_{n-1} := \delta_nA_n\sigma'(A_{n-1}x_{n-1}), \quad \partial_{A_{n-1}} s(x) = \delta_{n-1} L_{[x_{n-2}]} \end{equation*}\] von der letzten bis zur ersten layer berechnet werden können.

11.3 Praktische Berechnung des Gradienten

Im schönsten machine learning-Sprech sagt Jan hier immer Gradient obwohl eher die Ableitung oder die Jacobi Matrix gemeint ist.

Mit \(\hat y = -2(y-s(x;p)) \in \mathbb R^{n_2\times 1}\) bekommen wir aus obigen Formeln, dass \[\begin{equation*} \begin{split} \partial_{A_2}l((x,y); p) = \bigl \langle \hat y, \, \partial_{A_2}s(x;p)\bigr\rangle &= \hat y^T \sigma'(A_2x_1)L_{[x_1]} \\ &= \tilde y^TL_{[x_1]}, \end{split} \end{equation*}\] mit \(\tilde y=\sigma'(A_2x_1)\hat y\), was wegen der Diagonalität von \(\sigma'\) nur eine Skalierung der einzelnen Elemente von \(\hat y\) darstellt.

Eine Betrachtung der Dimensionen \(\hat y \in \mathbb R^{n_2 \times 1}\), \(\sigma'(A_2x_1)\in \mathbb R^{n_2\times n_2}\) und \(L_{[x_1]}\in \mathbb R^{n_2\times(n_2\times n_1)}\) ergibt für das Produkt, dass \[\begin{equation*} \tilde y^T L_{[x_1]} \in \mathbb R^{1\times (n_2 \times n_1)} \end{equation*}\] sodass wir schon fast den Gradienten zum Parameter \(A_2 \in \mathbb R^{n_2 \times n_1}\) addieren können.

Wir berechnen die Einträge \(g_{1ij}\) des Gradienten \(\tilde y^T L_{[x_1]}\) als die \(1\)-te (und einzige – da der Bildraum 1-dimensional ist) Komponente (bzgl. die kanonischen Basis) auf die der \(ij\) kanonische Basisvektor aus dem Urbildraum abgebildet wird. Hier bietet sich als kanonische Basis \(\{e^{(ij)}\}\subset \mathbb R^{n_2\times n_1}\) die Menge der Matrizen an, die \(0\) sind bis auf eine \(1\) in der \(i\)-ten Zeile an der \(j\)-ten Stelle. In der Tat gilt dann \[\begin{equation*} \mathbb R^{n_2\times n_1}\ni A = [a_{ij}] \leftrightarrow A = \sum_{i,j}a_{ij}e^{(ij)} \end{equation*}\] und weiter \[\begin{equation*} g_{1ij}=\tilde y^TL_{[x_1]}e^{(ij)} = \tilde y^T e^{(ij)}x_1 = \tilde y_i (x_{1})_j \end{equation*}\] sodass wir als \(1\)-te und einzige (matrixwertige) Komponente des Gradienten das äußere Produkt von \(\tilde y\) und \(x_1\) erhalten \[\begin{equation*} (\tilde y^T L_{[x_1]})_1 = \tilde y x_1^T \in \mathbb R^{n_2 \times n_1}. \end{equation*}\]

Für die Ableitung bezüglich \(A_k\) in einem \(N\)-layer Netzwerk bekommen wir mit obiger Rekursionsformel \[\begin{equation*} \partial_{A_k} l(x;p) = \hat y^T \delta_k L_{[x_{k-1}]} = \tilde{\tilde y} x_{k-1}^T \end{equation*}\] mit \[\begin{equation*} {\tilde {\tilde y}}^T = \hat y^T \sigma'(A_Nx_N)A_N\sigma'(A_{N-1})A_{N-1} \dotsm \sigma'(A_kx_{k-1}). \end{equation*}\]

11.4 Implementierungen und Beispiele

Jan könnte als grobe Einschätzung sagen:

Die Eleganz und Effizienz von AD für Optimierungsanwendungen ist unbestritten. Dennoch bedeutet der Einsatz eine Abwägung von Extra-Implementierungsaufwand und Selbstbeschränkung in Struktur und Funktionsumfang des Codes (das sind zwei Seiten der gleichen Medaille – auch nonstandard but clever Funktionen können mit entsprechendem Aufwand A-differenziert werden). Deswegen ist AD in der Forschung (hier wird der Mehraufwand gescheut) und in der industriellen Anwendung (viel legacy code mit spezifisch cleveren Funktionen) eine Randerscheinung wenn auch mit leuchtenden Erfolgsgeschichten.

Sicherlich ist es am besten im Bezug auf Aufwand und Funtionalität, den code sofort mit Augenmerk auf AD zu entwickeln.

Soll der AD overhead nicht selbst mitentwickelt werden, empfiehlt es sich auf Programmumgebungen zurückzugreifen, die einen AD framework mitbringen. Hier greift obige Abwägung – die AD verursacht kaum Mehraufwand allerdings muss der eigene Code an die vorgegebene Struktur und Funktionalität angepasst sein.

Für bestehenden Code gibt es zwei Herangehensweisen (neben der “Um”entwicklung).

1. precompiling – “AD ready” Code wird automatisiert aus bestehendem Code generiert.

2. overloading – Es werden AD Variablentypen und Klassen definiert und die Grundoperationen überladen sodass mit minimalen Änderungen ein bestehender Code die AD Funktionalität bekommt.

11.4.1 Anwendungsbeispiele

• Algorithmic differentiation of an industrial airfoil design tool coupled with the adjoint CFD method – ein Beispiel in dem ein kompletter Simulationscode automatisch differenziert wurde um Tragflügel zu optimieren.

11.4.2 Implementierungen

• CASADI – eine AD aware (für python und Matlab/Octave) Umgebung zur Optimierung (auch mit Differentialgleichungen)
• PyTorch – eine machine learning toolbox mit AD integriert zur Gradientenberechnung
• ADOL-C – ein overloading framework für C und C++
• Tapenade – ein precompiler der aus Code AD Code generiert

11.5 Aufgaben

11.5.1 Ableitung und Newton der Riccati Gleichung

Weisen Sie nach, dass der Gradient der (matrixwertigen) Riccati Abbildung (11.1) tatsächlich die gegebene Form hat und formulieren Sie damit ein Newton Verfahren zur Lösung der Riccati Gleichung \(f(S)=0\).

11.5.2 Ableitung von \(Ax\) nach der Matrix

Sei \(A=[a_{ij}]\in \mathbb R^{m\times n}\). Bestimmen Sie partielle Ableitung der Funktion \(f_x(A)=A(x)\) nach einer Komponente \[\begin{equation*} \frac{\partial {f_x}}{{\partial a_{ij}}}(A) \end{equation*}\] und geben Sie davon ausgehend eine Darstellung der totalen Ableitung \(\partial f_x(A) \colon \mathbb R^{m\times n}\to \mathbb R^{m}\) an

11.5.3 Backward Propagation

Leiten Sie aus \(\mathbf A = \begin{bmatrix} A & b \\ 0 & 1 \end{bmatrix}\) und \(\mathbf x = (x, 1)\) sowie \(\mathbf y = (y, 1)\) an der Stelle \([A\,b]=[A_0\, b_0]\) die Formel für den Gradienten \[\begin{equation*} \partial_{[A\,b]} (\frac 12 \|\mathbf y-\sigma(\mathbf A\mathbf x)\|_2^2)\Bigr|_{[A_0\, b_0]}=\Bigl (\sigma'(A_0x+b_0)\bigr (y-\sigma(A_0x+b_0)\bigr )\Bigr)\,\mathbf x^T \end{equation*}\] Bestätigen Sie sie numerisch, beispielsweise mit Hilfe der Code snippets aus Beispiel Implementierung, für ein Beispiel Netzwerk \(s\colon \mathbb R^{3}\to \mathbb R^{2}\colon \xi \to \sigma(\tilde A\xi + b)\) mit \(\sigma = \tanh\) an der Stelle \(A_0=\begin{bmatrix}2 & 1 & 1 \\ 1 & 2& 0 \end{bmatrix}\) und \(b_0=\begin{bmatrix}0\\1 \end{bmatrix}\) beim Datenpunkt \((\tilde x, y) = ((1,1,1), (1, 1))\).

11.5.4 AD and Autograd in pytorch

Arbeiten Sie sie sich durch die pytorch tutorials

1. Introduction to torch.autograd
2. AD with torch.autograd

und erweitern Sie die Codes um den Gradienten (in Bezug auf \(A\in \mathbb R^{2\times 2}\)) von \(x\to x^TAx\) für \(x\in \mathbb R^{2}\) an der Stelle \(A_0=\begin{bmatrix}2 & 1 \\ 1 & 2 \end{bmatrix}\) zu bestimmen.

Referenzen

Kurdila, A.J., Zabarankin, M.: Convex functional analysis. Birkhäuser (2005)

Tröltzsch, F.: Optimale steuerung partieller differentialgleichungen. Vieweg+Teubner, Wiesbaden, Germany (2009)