Die Maxwellgleichungen anschaulich verstehen

Die Maxwellgleichungen beschreiben die grundlegenden physikalischen Gesetze im Zusammenhang mit Elektrizität und Magnetismus. 

James Clerk Maxwell

Die mathematische Formulierung  erscheint zunächst schwierig, die Gleichungen lassen sich aber anschaulich verstehen.

\[\begin{flalign} & (1) \qquad \operatorname{div}\mathbf{E} = \frac{\rho}{\varepsilon_0} \\ & (2) \qquad \operatorname{div} \mathbf{B} = 0 \\ & (3) \qquad \operatorname{rot} \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} \\ & (4) \qquad \operatorname{rot} \mathbf{B} = \mu_0 \left(\mathbf{J} + \varepsilon_0\frac{\partial \mathbf{E}}{\partial t} \right) \end{flalign}\]

Bertrachten wir zunächst die erste Gleichung, das sogenannte Gaußsche Gesetz.

  \[\operatorname{div} \mathbf{E} = \frac{\rho}{\varepsilon_0}\]

Die Divergenz \(\operatorname{div}\) eines Vektorfeldes \(\mathbf{E}\) ist dabei als Skalarprodukt des Nabla Operators \( \nabla = \left( \frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z}\right)\) mit dem Feld \(\mathbf{E}\), also als \(\nabla \cdot \mathbf{E}\) definiert.

Entfernen wir uns im zwei dimensionalen Raum einen kleinen Schritt von der Position \(\left(x_0,y_0\right)\), so zeigt das Vektorfeld in der Regel in eine andere Richtung. Der rote Vektor im Bild oben ist der Differenzvektor.

Betrachten wir alle Differenzvektoren in der Umgebung von \(\left(x_0,y_0\right)\), und mitteln über die Skalarprodukte aus Abstandsvektor und Differenzvektor, so ergibt sich im Grenzübergang zu sehr kleinen Abständen \(\left( \partial / \partial x, \partial / \partial y\right) \cdot \left(E_x, E_y\right) \). Da sich das Skalarprodukt aus dem Produkt der Beträge der Vektoren mal dem Cosinus des Zwischenwinkels errechnet, wird es maximal wenn Abstandsvektor und Differenzvektor in die gleiche Richtung zeigen. Der Ausdruck \(\nabla\cdot \mathbf{E}\) beschreibt also wieviel 'Feld' aus dem Punkt \(\left(x_0,y_0\right)\) gemittelt über alle Richtungen 'herausfließt'. Das Gaußsche Gesetzt sagt nun, dass die 'Menge des herausfließenden elektrischen Feldes' proportional zur Punktladungsdichte ist. Punktladungen sind also die Quellen und Senken des elektrischen Feldes.

Die zweite Gleichung 

\[\operatorname{div} \mathbf{B} = 0\]

wird häufig auch als Gaußsches Gesetz des Magnetismus bezeichnet. Die Divergenz der magnetischen Flussdichte \(\mathbf{B}\) ist Null. Es gibt keine Quellen oder Senken, d.h. keine magnetischen Monopole, die analog zu elektrostatischen Ladungen wirken würden. Das magnetische Feld verhält sich wie eine nicht komprimierbare Flüssigkeit. 

Die dritte Gleichung

   \[\operatorname{rot} \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}\]

wird als Faradaysches Induktionsgesetzt bezeichnet.

Als Rotation \(\operatorname{rot}\) des Vektorfeldes \(\mathbf{E}\) bezeichnet man dabei das Vektorprodukt des Nabla Operators \( \nabla = \left( \frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z}\right)\) mit dem Feld \(\mathbf{E}\), also \(\nabla \times \mathbf{E}\).

Im Gegensatz zum Skalarprodukt, ist das Vektorprodukt ein Maß dafür, ob zwei Vektoren zu einander senkrecht stehen. Sein Betrag errechnet sich als Produkt der Beträge, mal dem Sinus des Zwischenwinkels. Der Betrag des Vektorprodukts wird also für einen Zwischenwinkel von 90° maximal. Betrachten wir wieder die Differenzvektoren in einigem Abstand von \(\left(x_0,y_0\right)\), so entspricht die Rotation \(\operatorname{rot}\) dem Grenzwert des Mittelwertes aus Abstandsvektoren \(\times\) Differenzvektoren für kleine Abstände im Übergang zum Punkt \(\left(x_0,y_0\right)\). Das Induktionsgesetz besagt nun, dass eine zeitliche Änderung des Magnetfeldes bzw. der magnetischen Flussdichte zu einer Rotation des elektrischen Feldes führt. Anders ausgedrückt resultiert also die zeitliche Änderung des durch eine Leiterschleife eingeschlossenen magnetischen Flusses in einer Spannung an den Enden der Leiterschleife. Die induzierte Spannung wirkt der Änderung entgegen. Das Bild unten von Michael Lenz veranschaulicht dies.

Die vierte Gleichung

\[ \operatorname{rot} \mathbf{B} = \mu_0 \left(\mathbf{J} + \varepsilon_0\frac{\partial \mathbf{E}}{\partial t} \right) \]

wird als Ampèresches Gesetz mit Maxwells Ergänzung bezeichnet. Der erste Term auf der rechten Seite der Gleichung besagt, dass die Rotation des Magnetfeldes bzw. die magnetische Flußdichte \( \mathbf{B}\) proportional zum elektrischen Stromfluss bzw. der Stromdichte \( \mathbf{J}\) ist. Ist das elektrische Feld \(\mathbf{E}\) zeitlich konstant, können wir den zweiten Term vernachlässigen. Für einen von Gleichstrom durchflossenen Leiter gilt also das untenstehende Bild.

Ändert sich das elektrische Feld \(\mathbf{E}\), ist die Rotation zusätzlich proportional zur zeitlichen Änderung des elektrischen Feldes. Diese Änderung wird häufig auch Verschiebungsstrom genannt, und stellt die sogenannte Maxwell Ergänzung dar.

Aus der Maxwell Ergänzung im Wechselstromfall

\[ \operatorname{rot} \mathbf{B} = \mu_0 \left(\varepsilon_0\frac{\partial \mathbf{E}}{\partial t} \right) \]

lässt sich durch erneutes Anwenden des Rotationoperators recht einfach die Wellengleichung

$$\Delta \mathbf{E} = \mu_0 \varepsilon_0 \frac{\partial^2 \mathbf{E}}{\partial t^2} $$

herleiten, \(\Delta\) bezeichnet dabei den sogenannten Laplace Operator. Eine Lösung dieser Differentialgleichung sind ebene Wellen der Form

$$\mathbf{E} = E_0 f(k \cdot \mathbf{x} - ct)$$

Sie beschreiben die Ausbreitung einer elektromagnetischen Welle im Raum.

Viel Spaß beim selber Rechnen und Verstehen...

Ein neuronales Netzwerk von Grund auf programmieren

Es stehen mittlerweile viele Werkzeuge wie Keras, als Applikationsschnittstelle für Bibliotheken wie tensorflow, zur Verfügung um neuronale Netzwerke zu programmieren. Um ein besseres Verständnis der Zusammenhänge zu erlangen, möchte ich jedoch ein neuronales Netzwerk von Grund auf, d.h. im Wesentlichen nur unter Verwendung von Phyton und numpy programmieren.

Ein neuronales Netzwerk, genauer gesagt ein Feedforward Neural Network (FNN), besteht aus einer Eingabeschicht, einer oder mehreren verborgener Schichten und einer Ausgabeschicht.

Die einzelnen Schichten \(A\) beinhalten sogenannte Neuronen, die über lineare Abbildungen 

$$ Z^{l+1}=W^{l+1}\cdot A^l + B^{l+1}$$

miteinander verknüpft sind. Dabei bezeichnet \(W^{l+1}\) eine Matrix, deren Dimension der Zahl der Neuronen der Schicht \(l+1\) mal der Zahl der Neuronen der Schicht \(l\) entspricht. \(B^{l+1}\) bezeichnet einen Vektor der addiert wird. Sollen Bilder verarbeitet werden, stellen die Neuronen z.B. die einzelnen Pixel des Bildes dar. Die Matrizen \(W\) und \(B\) werden im Rahmen eines Trainings so angepasst, dass die Ausgabeschicht das gewünschte Ergebnis möglichst gut erreicht. Da durch lineare Abbildungen nur lineare Zusammenhänge 'angefittet' werden können, werden die errechneten \(Z \) Werte jeder Schicht durch nicht lineare Funktionen 'aktiviert', und dadurch in Werte von \(A\) abgebildet. Die einfachste, und sehr häufig verwende Aktivierungsfunktionen ist die sogenannte Rectified Linear Unit (ReLU) Funktion. Die ReLu Funktion liefert 0 wenn \(z<0\) und \(z\) wenn \(z\geq 0\).

Sollen z.B. Bilder klassifiziert werden, dient diese Funktion im Netzwerk dazu die Ausgabe in eine Wahrscheinlichkeitsverteilung über die vorherzusagenden Ausgabeklassen umzuwandeln. Für Mehrklassenprobleme wird im Allgemeinen die Softmax-Funktion verwendet.

\[ \text{softmax}(z_j)= s_j = \frac{e^{z_j}}{\sum_{i=1}^{k} e^{z_i}} \]

Dabei ist \(z\) der Eingabevektor mit \(k\) Komponenten, \(z_j\) die \(j\)-te Komponente von \(z\), und die Summe im Nenner läuft über alle \(k\) Komponenten. Grafisch lässt sich die Funktion in Abhängigkeit einer Komponente wie folgt darstellen:

Die Funktion kann als normalisierte Exponentialfunktion betrachtet werden, und ist eine Verallgemeinerung der logistischen Funktion auf mehrere Dimensionen. Bei welchem \(z\) Wert der 'Übergang' von 0 nach 1 einer Komponente stattfindet, hängt von den Werten der anderen Komponenten ab. Durch die Summe im Nenner ist jedoch stets gewährleistet das die Summe aller Einzelwahrscheinlichkeiten 1 beträgt.

Das Training des neuronalen Netzes läuft nun wie folgt ab. Zu Beginn des ersten Durchlaufs werden die Matrizen \(W\) und \(B\) zufällig initialisiert. Anschließend erfolgt die sogenannte 'forward propagation' mit aufbereiteten Trainingsdaten. Die errechneten Wahrscheinlichkeiten erreichen in der Regel nicht 100% für eine Kategorie. Wie stark die Ergebnisse abweichen, kommt in der Kostenfunktion zum Ausdruck. Soll mehr als eine  Kategorie bestimmt werden, empfiehlt es sich die sogenannte Kreuzentropie zweier Wahrscheinlichkeitsvektoren als Kostenfunktion zu verwenden.

$$ L(y,s) = - \sum_{i=1}^{C} y_i \cdot \log(s_i) $$

Dabei bezeichnet \( C \) die Anzahl der Klassen, \(s_i\) die Ausgabewerte der letzten Schicht, und \(y_i\) den Ergebnisvektor in 'one hot' oder '1 aus n' Kodierung.  Je kleiner die Kreuzentropie, desto ähnlicher sind die beiden Wahrscheinlichkeitsverteilungen, denn \(-log(s_i)\) wird umso kleiner je mehr sich \(s_i\) der 1 annähert.

Im weiteren Verlauf des Trainings werden die Parameter \(W\) und \(B\) nach der sogenannten 'gradient descent' Methode optimiert. Dazu wird die Änderung der Kostenfunktion für jeden Layer und jeden Parameter bestimmt. Anschließend wird das neuronale Netzwerk in der 'back propagation' durchlaufen, und die Parameter für den nächsten Trainingsdurchlauf angepasst. Der letzte nicht aktivierte Layer \(ll\) soll sich also iterativ wie folgt ändern:

$$ Z^{ll}:=\nabla_{Z^{ll}}L = Z^{ll}-\alpha\frac{dL}{dZ^{ll}}$$

Dabei bezeichnet \(\alpha\) die sogenannte Lernrate. Um die Ableitung der Kostenfunktion \(L\) nach \(Z\) bestimmen zu können, müssen wir zunächst die Softmaxfunktion ableiten. Das erfordert etwas Aufwand, aber wenn wir die partielle Ableitung des Logarithmus bilden, können wir die Rechnung vereinfachen.

$$ \frac{\partial}{\partial z_j} \log(s_i) = \frac{1}{s_i} \cdot \frac{\partial s_i}{\partial z_j}, $$

Der Ausdruck auf der rechten Seite folgt direkt aus der Kettenregel. Als Nächstes ordnen wir die obere Formel um.

$$\frac{\partial s_i}{\partial z_j} = s_i \cdot  \frac{\partial}{\partial z_j} \log(s_i)$$

Die linke Seite ist genau die partielle Ableitung, die wir suchen. Wir betrachten zuerst den Logarithmus von s.

$$\begin{align*}\log s_i & = \log \left( e^{z_i} / \sum_{i=1}^n e^{z_i} \right) = z_i - \log \left( \sum_{i=1}^n e^{z_i} \right)\end{align*}$$

Die partielle Ableitung des resultierenden Ausdrucks ist:

$$\frac{\partial}{\partial z_j}\log s_i = \frac{\partial z_i}{\partial z_j}-\frac{\partial}{\partial z_j}\log \left( \sum^n_{i=1}e^{z_i}\right)$$

Der erste Term auf der rechten Seite

$$ \frac{\partial z_j}{\partial z_j}= \begin{cases}1, & i = j \\ 0, & \text{sonst} \end{cases}$$

kann prägnant mit der Indikatorfunktion \(1\{i=j\}\) geschrieben werden. Die Indikatorfunktion nimmt den Wert 1 an, wenn \(i=j\) ist, und sonst 0. Der zweite Term auf der rechten Seite kann durch Anwendung der Kettenregel ausgewertet werden.

$$\frac{\partial}{\partial z_j} \log s_i = 1\{i = j\} - \frac{1}{\sum_{l=1}^n e^{z_l}} \left( \frac{\partial}{\partial z_j} \sum_{l=1}^n e^{z_l}\right)$$

Die partielle Ableitung der Summe zu erhalten, ist trivial.

$$\frac{\partial}{\partial z_j} \sum_{i=1}^n e^{z_i} = \frac{\partial}{\partial z_j} [e^{z_1} + e^{z_2} + ... + e^{z_j} + ... + e^{z_n}] = \frac{\partial}{\partial z_j} [e^{z_j}] = e^{z_j}$$

Eingesetzt in die obige Formel führt dies zu:

$$\frac{\partial}{\partial z_j} \log s_s = 1\{i=j\} - \frac{e^{z_i}}{\sum_{i=1}^n e^{z_i}} = 1\{i=j\}-s_j$$

Schließlich müssen wir den obigen Ausdruck mit \(s_i\) multiplizieren wie am Anfang des Abschnitts gezeigt:

$$\frac{\delta s_i}{\delta z_j}=s_i \frac{\partial}{\partial z_j}\log(s_i) =s_i \left( 1\{i=j\}-s_j \right)$$

Damit ist unsere Ableitung abgeschlossen. Wir haben eine Formel für alle Elemente der Jacobi-Matrix \(J\) erhalten. Für den Spezialfall von n = 4 gilt zB.:

$$J_{\text{softmax}} = \begin{bmatrix}s_1 \cdot (1 - s_1) & -s_1 \cdot s_2 & -s_1 \cdot s_3 & -s_1 \cdot s_4 \\ -s_2 \cdot s_1 & s_2 \cdot (1 - s_2) & -s_2 \cdot s3 & -s_2 \cdot s_4 \\ -s_3 \cdot s_1 & -s_3 \cdot s_2 & s_3 \cdot (1 - s_3)  &  -s_3 \cdot s4 \\ -s_4 \cdot s_1 &  -s_4 \cdot s_2  &  -s_4 \cdot s_3 & s_4 \cdot (1 - s_4) \end{bmatrix}$$ 

Als nächstes wollen wir nun die Ableitung der Kostenfunktion in Bezug auf die Eingabe \(z\) der letzten Aktivierungsschicht berechnen.

$$\frac{\partial L}{\partial z_j} = - \frac{\partial}{\partial z_j}\sum_{i=1}^{c} y_i \cdot \log(s_i) = - \sum_{i=1}^{c} y_i \cdot \frac{\partial}{\partial z_j}\log(s_i) = -\sum_{i=1}^{c} \frac{y_i}{s_i} \cdot \frac{\partial s_i}{\partial z_j}$$

Fügen wir die Ableitungen ein, die wir oben erhalten haben,

$$\frac{\partial L}{\partial z_j} = -\sum_{i=1}^{c} \frac{y_i}{s_i} \cdot s_i \cdot (1\{i=j\}-s_j) = -\sum_{i=1}^{c} y_i \cdot (1\{i=j\}-s_j)$$

und multiplizieren aus

$$\frac{\partial L}{\partial z_j}=  \sum_{i=1 }^c y_i \cdot s_j-\sum_{i=1 }^c y_i \cdot 1\{i=j\})$$

dann ergibt sich

$$\frac{\partial L}{\partial z_j} = \sum_{i=1}^{c} y_i \cdot s_j - y_j$$

Als nächstes ziehen wir \(s_i\) aus der Summe, da es nicht vom Index \(i\) abhängt

$$\frac{\partial L}{\partial z_j} = s_j \cdot \sum_{i=1}^c y_i -y_j = s_j - y_j$$.

In prägnanter Vektorschreibweise erhalten wir:

$$\nabla_{Z^{ll}}L = S -Y$$

Als nächstes wollen wir bestimmen, wie sich \(W\) nach einem Trainingsdurchlauf ändern soll. Für die Änderung von \(W\) soll analog zur Änderung von \(Z\) gelten:

$$W^l:=W^l-\alpha '\frac{dL}{dW^l}=W^l-\alpha '\frac{dL}{dZ^l}\cdot\frac{dZ^l}{dW^l}$$

Der Ausdruck \(dZ/dW\) lässt sich mathematisch ähnlich schwer formulieren wie der bereits berechnete Ausdruck \(dL/dZ\). Es gilt:

$$\frac{\partial z_{ik}^l}{\partial w_{mn}^l}=\frac{\partial}{\partial w_{mn}^l}\sum_{j=1}^{\text{dim(layer l)}}w_{ij}^l a_{jk}^{l-1} + b_{ik}^l=1\{m=i,n=j\}a_{nk}^{l-1}$$

Dabei bezeichnet \(a_{nk}\) das Element einer Matrix der Dimension \(dim(\text{layer l}) * N\), wobei \(N\) die Zahl der Trainingsdatensätze darstellt. Durch die Verwendung des Skalarprodukt \(dL/dZ \cdot dZ/dW\) möchten wir die Änderung in \(W\) über alle Datensätze mitteln. Wir normieren daher auf \(1/N\) und transponieren \(A^{l-1}\). Es gilt:

$$\lambda ' \left( \frac{dL}{dZ^l}\cdot \frac{dZ^l}{dW^l} \right)_{mn}=\frac{\lambda}{N}\sum_{k=1}^N \left( \frac{dL}{dz^l}\right)_{mk}\cdot a_{kn}^{l-1}$$ 

oder übersichtlicher

$$\nabla_{W^l}L = \frac{1}{N} \nabla_{Z^l}L\cdot {A^{l-1}}^T$$

Auch für die iterative Änderung von \(B\) gilt:

$$B^l:=B^l-\alpha '\nabla_{B^l}L=B^l-\alpha '\frac{dL}{dB^l}= B^l-\alpha '\frac{dL}{dZ^l}\cdot\frac{dZ^l}{dB^l}=B^l-\alpha '\frac{d L}{dZ^l}^l\cdot 1$$

Dabei bezeichnet \(1\) eine Einheitsmatrix der Dimension \(dim(\text{layer l}) * N\). In Vektorschreibweise gilt wieder:

$$\nabla_{B^l}L=\frac{1}{N}\nabla_{Z^l}L\cdot 1^T$$

Wir multiplizieren erneut mit der transponierten (Einheits)matrix um über alle Datensätze zu mitteln.

Zu guter Letzt gilt für das update der inneren Schichten \(Z\):

$$Z^l:=Z^l-\alpha '\nabla_{Z^l}L = Z^l-\alpha '\frac{dL}{dZ^l}=Z^l-\alpha '\frac{dZ^{l+1}}{dA^l}\cdot\frac{dL}{dZ^{l+1}}\cdot\frac{dA^l}{dZ^l}$$

und in Vektorschreibweise bedeutet das:

$$\nabla_{Z^l}L={W^{l+1}}^T\cdot\nabla_{Z^{l+1}}L\odot\text{ReLu}'(Z^l)$$

Dabei steht \(\odot\) für die elementweise Multiplikation bzw. das Hadamard Produkt und \(ReLu '\) für die Ableitung der ReLu Aktivierungsfunktion. Diese Ableitung ist trivialerweise durch

$$ReLu'(z) = \begin{cases} 1 & \text{wenn } z > 0 \\0 & \text{wenn } z \leq 0 \end{cases}$$

gegeben. Damit sind die mathematischen Grundlagen unseres neuronalen Netzwerkes im Wesentlichen erarbeitet, und wir können uns an die Programmierung wagen.

Unser neuronales Netzwerk soll als Trainingsdaten den berühmten MNIST Datensatz verwenden. Dieser besteht aus 42000 handgeschriebenen Ziffern in 8 bit Graustufen mit der Auflösung 28 x 28 Pixel. Die Eingabeschicht \(A^0\) besteht daher aus 784 Einheiten, die den 784 Pixeln in jedem 28 x 28-Eingabebild entsprechen. Der Programmcode ist so gestaltet, dass beliebig viele innere Schichten eingefügt werden können. Bereits mit nur einer inneren Schicht lassen sich jedoch sehr gute Ergebnisse erzielen. Diese eine versteckte Schicht \(A^1\) besteht aus 10 Einheiten mit ReLU-Aktivierung, und die Ausgabeschicht \(A^2\) besteht ebenfalls aus 10 Einheiten, die den zehn Ziffernklassen mit Softmax-Aktivierung entsprechen.

Der gesamte Code steht als Kaggle notebook zum Ausprobieren zur Verfügung. Ich habe versucht die einzelnen Schritte möglichst gut im Programmcode zu kommentieren.

 

import numpy as np
import pandas as pd
import time
from datetime import datetime
from matplotlib import pyplot as plt


# set structure of neural network
layer_dim = [784,10,10] # number of neurons per layer (including input and output layer)
L = len(layer_dim) # number of layers
LL = L-1 # index of Last Layer

# set training parameter
test_size = 1000 # set number of test samples, note: test dataset  will be removed from training dataset
training_rate = 0.10 # set training rate
iterations = 200 # set number of iterations

# set global variables 
params = {} # create dictionary for parameter W and b
cache = {} # create dictionary for A and Z
grads = {} # create dictionary for parameter gradients
cost = np.zeros(iterations) # create array for cost value per iteration
accuracy = np.zeros(iterations) # create array for achieved training accuracy per iteration

# load and format datasets

data = pd.read_csv('/kaggle/input/training-set/train.csv') # import *.csv using pandas
data = np.array(data) # convert data into numpy array
m, n = data.shape # m = number of samples, n= number of input neurons + label

np.random.shuffle(data) # shuffle lines before splitting into test and training sets

# create test dataset and reformat
data_test = data[0:test_size].T # take 1st test_size rows and transpose 
Y_test = data_test[0] # 1st row contains labels
X_test = data_test[1:n] / 255. # remaining rows contain normalized gray values

# create training dataset and reformat
data_train = data[test_size:m].T # treat training set in the same way
Y_train = data_train[0]
X_train = data_train[1:n] / 255.

# define one hot representation of label vector
def one_hot(Y):
    one_hot_Y = np.zeros((Y.max() + 1, Y.size)) # create zero matrix w/ shape (number of categories) x (number of samples)
    one_hot_Y[Y, np.arange(Y.size)] = 1 # set 1 in row Y for every column (from 0 to number of samples -1)
    return one_hot_Y

Oh_Y_train = one_hot(Y_train) # calculate one_hot_y and set as a global variable
Oh_Y_test = one_hot(Y_test)


#define the neural network

# define initial model parameter
def init_params():
    for l in range (1,L):
        params['W'+str(l)] = np.random.normal(size=(layer_dim[l], layer_dim[l-1])) * np.sqrt(1./layer_dim[l-1]) # create array of dim(layer_n+1) x dim(layer_n) with values around 0 and standard deviation of 1 scaled to sqrt(1/# neurons layer_n) 
        params['b'+str(l)] = np.random.normal(size=(layer_dim[l], 1)) * np.sqrt(1./layer_dim[l]) # create array of dim(layer_n+1) x 1 with values around 0 and standard deviation of 1 scaled to sqrt(1/# neurons layer_n+1) 
    return params

# define non linear activation functions
# define Rectified Linear Unit function for activation of inner layers 
def ReLU(Z):
    return np.maximum(Z, 0)
# define softmax function to calculate propability for last layer (logistic regression)
def softmax(Z):
    A = np.exp(Z) / sum(np.exp(Z))
    return A

# define forward propagation
def forward_prop(X, params):
    cache['A0'] = X
    for l in range (1,LL): # calculate inner layers
        cache['Z'+str(l)] = params['W'+str(l)].dot(cache['A'+str(l-1)]) + params['b'+str(l)] # Z = W x A + b
        cache['A'+str(l)] = ReLU(cache['Z'+str(l)]) # A = ReLu(Z)
    cache['Z'+str(LL)] = params['W'+str(LL)].dot(cache['A'+str(LL-1)]) + params['b'+str(LL)] # calculate outer layer
    cache['A'+str(LL)] = softmax(cache['Z'+str(LL)]) # A_Last_Layer = softmax (Z_Last_Layer) 
    return  cache['A'+str(LL)], cache

# define derivative of Rectified Linear Unit function
def ReLU_deriv(Z):
    return Z > 0 # return all Z values > 0 , else return 0

# define backward propagation
def backward_prop(params, cache, Oh_Y):   
    grads['Z'+str(LL)] = cache['A'+str(LL)] - Oh_Y # calculate change of Z_LastLayer, Note: repeat {Z = Z -alpha * dL/dZ} for gradient descent, and dL/dA * dA/dZ = A-Y when using L = -(Y * ln(A) + (1-Y) * ln(1-A)) as a cost function, and softmax to activate A
    grads['W'+str(LL)] = 1 / m * grads['Z'+str(LL)].dot(cache['A'+str(LL-1)].T) # calculate change of W_LL , dW_l = 1/m * dZ_l * A_(l-1).T 
    grads['b'+str(LL)] = 1 / m * np.sum(grads['Z'+str(LL)]) # calculate change of b_LL, dB_l = 1/m * sum(dZ_l)
    for l in range (LL-1,0,-1):  # calculate change of inner layers    
        grads['Z'+str(l)] = params['W'+str(l+1)].T.dot(grads['Z'+str(l+1)]) * ReLU_deriv(cache['Z'+str(l)]) # calculate change of Z_l = W_(l+1).T x dZ_(l+1) * dReLu(Z_l)
        grads['W'+str(l)] = 1 / m * grads['Z'+str(l)].dot(cache['A'+str(l-1)].T) # calculate change of W_l
        grads['b'+str(l)] = 1 / m * np.sum(grads['Z'+str(l)]) # calculate change of b_l
    return grads

# define update of parameters
def update_params(params, grads, alpha):
    for l in range (1,L):
        params['W'+str(l)] = params['W'+str(l)] - alpha * grads['W'+str(l)] # gradient descent W_l, W_l = W_l - alpha * dW_l
        params['b'+str(l)] = params['b'+str(l)] - alpha * grads['b'+str(l)] # gradient descent b_l, b_l = b_l - alpha * dB_l
    return params


# define evaluation of model
# define prediction function
def get_predictions(AL):
    return np.argmax(AL, 0) # position of biggest value in Activated_LastLayer is predicted number

# define calculation of accuracy
def get_accuracy(predictions, Oh_Y):
    return np.sum(predictions == np.argmax(Oh_Y, 0)) / Oh_Y.shape[1] # sum up matches and normalize to samplesize

# define calculation of cost function
def compute_cost(AL, Oh_Y):
    cost = (-1 / m) * np.sum(np.multiply(Oh_Y, np.log(AL)))
    return (cost)

# define cost and accuracy plot   
def plot_graph(x_values, cost, accuracy):   
    fig, axs = plt.subplots(1, 2, figsize=(9, 3))
    axs[0].set_xlabel('iteration')
    axs[0].set_ylabel('cost [a.u.]')
    axs[0].plot(x_values,cost)
    axs[1].set_xlabel('iteration')
    axs[1].set_ylabel('accuracy [%]')
    axs[1].set_ylim(0,100)
    axs[1].plot(x_values,accuracy)
    plt.show()

Nachdem alle notwendigen Funktionen definiert sind, starten wir das Training.

#start the training
start = time.time() 
params = init_params() # initialize W and b
for i in range(iterations):
    AL, cache = forward_prop(X_train, params) # perform forward propagation
    cost[i] = compute_cost(AL,Oh_Y_train) # calculate cost
    accuracy[i] = get_accuracy(get_predictions(AL), Oh_Y_train) # get average accuracy for all samples
    grads = backward_prop(params, cache, Oh_Y_train) # perform backward propagation
    params = update_params(params, grads, training_rate) # update model
    # print out cost and accurracy
    if i % 20 == 0:
        print("Iteration: "+str(i))
        print("Cost :"+str(cost[i])+" Accuracy: "+str(accuracy[i]))
end = time.time()
duration = end - start # training duration in seconds

plot_graph(list(range(0,iterations)), cost, accuracy*100) # plot cost and accuracy vs iterations

Bereits nach 200 Iterationen sättigt die Genauigkeit merklich. Wir erreichen eine Genauigkeit von ~85%.

Wir sichern das Model,

# Save model
now = datetime.now()
date = now.strftime('%Y-%m-%dT%H-%M-%S') + ('-%02d' % (now.microsecond / 10000))
name = ""
for i in range (len(layer_dim)):
    name += str(layer_dim[i])+"_"
name = date+"_dim_"+name+"samples_"+str(m)+"_rate_"+str(training_rate).replace('.', '')+"_iterations_"+str(iterations)
np.save(name, params) # save to ./kaggle/working

und können es wieder laden.

# Load model
model = np.load(name +".npy",allow_pickle='TRUE').item()

# define how to predict label for a certain sample
def make_prediction(X, model):
    AL ,_ = forward_prop(X, model)
    prediction = get_predictions(AL)
    return prediction

# define how to test prediction vs label and image
def test_prediction(index, model):
    current_image = X_test[:, index, None].reshape((28, 28)) * 255
    prediction = make_prediction(X_test[:, index, None], model)
    label = Y_test[index]
    print("Prediction: ", prediction)
    print("Label: ", label)
    plt.gray()
    plt.imshow(current_image, interpolation='nearest')
    plt.show()
    return

Dann berechnen wir die für den Testdatensatz erzielte Genauigkeit. Alternativ könnten wir uns auch einige Testbilder und die aus dem Model bestimmte Kategorie anzeigen lassen. 

# calculate overall accuracy of model for test dataset
test_accurracy = get_accuracy(make_prediction(X_test, model), Oh_Y_test)
print(test_accurracy)

Für den obigen Durchlauf erzielen wir auch für den nicht trainierten Testdatensatz eine Genauigkeit von fast 89%. Um die Performance der untersuchten Konfiguration zu dokumentieren, können wir abschließend noch ein log file speichern.

# write log file for training session
log = ["date: "+date+"\n",
       "training samples: "+str(m)+"\n",
       "training duration: "+str(round(duration/60, 1))+" minutes\n",
       "training accuracy: "+str(round(accuracy[iterations-1]*100, 1))+" %\n",
       "test samples: "+str(test_size)+"\n",
       "test accurracy: "+str(round(test_accurracy*100, 1))+" %\n"]
with open("log_"+str(name)+".txt", 'w') as f:
  f.writelines(log)

Viel Spaß beim selber programmieren...

Die Signalstärkeanzeige eines Funkgerätes kalibrieren

Wie viele auf dem Markt verfügbaren Funkgeräte liefert mein Kenwood TH-79E Handfunkgerät nur eine sehr rudimentäre Aussage über die Stärke eines empfangenen Signal. Unten ein Auszug aus dem Handbuch.

Weder im Handbuch noch im Servicemanual findet sich eine Aussage über den quantitativen Wert in µV oder dBm.

Im Funkverkehr wird eine Funkverbindung üblicherweise mit einem Rapport bestätigt. Aus der Morsewelt hat man das R(eadability)-S(ignal strength)-T(one) System übernommen um die die Qualität der empfangene Aussendung zu beschreiben. Für die R-Werte gilt:

und für die S-Werte ab 144 MHz:

Der T-Wert wird nur für Telegrafiesignale jedoch nicht für Telefoniesendungen verwendet.

Für eine korrekte Einschätzung des S-Wertes ist also ein zusätzliches Messgerät, ein sogenanntes S-Meter notwendig. Ich möchte im Folgenden die 'Kalibrierung' der Anzeige meines Handfunkgeräts mit einen NanoVNA durchführen.

Ich habe zunächst die Ausgangsleistung des NanoVNA-H mit dem Oszilloskop an einem 50 Ohm Abschlusswiderstand bestimmt.

Bei 1 MHz zeigt sich eine Peak zu Peak Spannung von Vpp = 162mV. Die effektive Spannung U_eff des Rechtecksignal beträgt Vpp/2 = 81mV. Der NanoVNA liefert an 50 Ohm also eine Leistung von P = U^2/R = 128µW oder ca. -9 dBm. Ich konnte diese Signalstärke auch bei 100kHz, 10MHz und 100MHz bestätigen.

Um einen S-Wert von z.B. 9 zu erzielen sind also zusätzliche Dämpfungsglieder notwendig. Ich habe mir daher analoge Dämpfungsglieder im Verbund von 0, -10, -20, und -30 dB, 

sowie eines digitales Abschwächglied das von 0 bis -31.5 dB in Stufen von 0,5 dB einstellbar ist, zugelegt.

Zur Verbesserung des Übersprechverhaltens habe ich die SMA Buchsen beider Produkte auf der Rückseite zusätzlich verlötet.

In einem ersten Schritt habe ich die einzelnen Dämpfungsglieder in einem Frequenzbereich von 1 MHz bis 500 MHz charakterisiert. Bei einer nominellen Dämpfung von 10dB

habe ich Werte zwischen -9.79 dB und -8,98 dB gemessen. Bei nominell -20dB

Werte zwischen -19,54 dB und -18,65dB, und bei nominell -30dB

Werte zwischen -29,49 dB und -28,25dB. Alles in Allem stimmen die Werte gut mit dem Aufdruck überein, und sogar bei unprofessioneller Verschaltung aller Dämpfungsglieder mit Drahtbrücken 

lässt dich eine Gesamtdämpfung von 60dB mit einer Genauigkeit von +/- 1dB erzielen.

Das Digitale Dämpfungsglied habe ich bei einer Versorgungsspannung von 3,1V mit offenen Digitaleingängen auf 0 dB abgeglichen.

Auch für dieses Dämpfungsglied zeigt sich zwischen 1 MHz und 500 MHz mit +/- 1 dB  eine recht gute Übereinstimmung mit der nominellen Dämpfungskaskade.

Für die weitere Messung mit dem Handfunkgerät, hat sich der NanoVNA-H selbst jedoch auf Grund einer fehlenden Abschirmung als deutliche Störquelle herausgestellt. Um die Störung hörbar zumachen habe ich im Video unten einen Sweep von 145,995 kHz bis 146,005 kHz auf dem NanoVNA durchgeführt, und die Empfangsfrequenz am  Funkgerät auf 146 MHz eingestellt. 

Auch in mehreren Metern Entfernung ist die Abstrahlung des NanoVNA deutlich zu hören. Erst die Schirmung in einer Metallschachtel schafft Abhilfe.

Ich habe zunächst versucht den NanoVNA in die Schachtel zu bauen, 

bin dann aber dazu übergegangen die gesamte Elektronik in die Schachtel zu packen. Ich habe lediglich ein kurzes Kabel zum Funkgerät geführt, und die Löcher im Metallgehäuse mit Alufolie verschlossen. 

Für verschiedene Dämpfungskombinationen konnte ich dann folgende Skalenteile am Funkgerät abgelesen.

In der Tabelle oben habe ich zusätzlich die rechnerische Signalstärke und die 'Peak to Peak' Spannung eines entsprechenden Rechtecksignals aufgeführt. Grafisch aufgetragen ergibt sich folgendes Bild.

Die Zahl der Balken auf dem Display entspricht recht genau den S-Werten. Ein ziemlich erstaunliches Ergebnis wie ich finde, wenn man sich vor Augen führt, dass Spannungen im sub-µV-Bereich gemessen werden müssen. Im Video unten ist die Messung bei einer Dämpfung von 120 dB zu sehen. Die Signalanzeige des Funkgerätes ändert sich beim Aufsetzen des Metalldeckels von Vollausschlag auf 2 Balken. Dies entspricht einem einem S-Wert von 2-3.

Viel Spaß beim selber experimentieren...

Die Fourier-Transformation und die Unschärferelation in der Quantenmechanik

Die Fourier-Transformation hat vielerlei Anwendungen in der Mathematik, den  Ingenieurswissenschaften und in der Physik. Insbesondere in der Quantenmechanik vermittelt sie zwischen dem sogenannten Ortsraum eines Systems und dem Impulsraum. In ihren Eigenschaften kommt unter Anderem die Unschärfe zwischen Ort und Impuls eines quantenmechanischen Teilchens zum Ausdruck.

Beginnen möchte ich meine Überlegungen mit der zeitunabhängigen Schrödingergleichung.

$$\left(-\frac{\hbar^2}{2m}\Delta_r + V(\vec{r})\right)\psi (\vec{r}) = E \psi (\vec{r})$$

Dabei bezeichnet V(r) das Potential in dem sich das Teilchen bewegt, \(E\) die Energie des Systems, und \( \psi (\vec{r})\) die Wellenfunktion des Teilchens im Ortsraum. Des Weiteren bezeichnet \(\hbar\) das Wirkungsquantum, \(m\) die Masse des Teilchens, und \(\Delta_r = \nabla^2\) den sogenannten Laplaceoperator. Um den Zusammenhang zwischen Impuls p und Wellenlänge \(\lambda\) bzw.  Wellenzahl k zu verstehen, ist es hilfreich den einfachen Fall einer freien Bewegung, d.h. mit \( V(\vec{r}) = 0\), in einer Dimension zu betrachten. Die Schrödingergleichung vereinfacht sich zu

$$-\frac{\hbar^2}{2m}\left(\frac{\partial}{\partial x}\right)^2\psi=E\psi=\frac{1}{2}mv^2\psi=\frac{p^2}{2m}\psi .$$

Aus dieser Darstellung wird unmittelbar ersichtlich, dass der Operator für den Impuls sinnvollerweise als 

$$ \hat{p}:=\frac{\hbar}{i}\frac{\partial}{\partial x} \quad \mathsf{oder} \quad \hat{p}:=\frac{\hbar}{i}\nabla_r$$

definiert werden kann. Wenden wir den Operator z.B. auf eine gebundene Kreisbewegung in der komplexen Zahlenebene an, dann gilt

$$\begin{eqnarray}\hat{p}\, \psi (r)& = & \frac{\hbar}{i}\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y}\right)\left[\cos\left(\frac{2\pi}{\lambda}x\right)+i\sin\left(\frac{2\pi}{\lambda}x\right)\right] \\ & = &\frac{2\pi}{\lambda}\frac{\hbar}{i}\left[-\sin\left(\frac{2\pi}{\lambda}x\right)+i\cos\left(\frac{2\pi}{\lambda}x\right)\right] \\ & = & \frac{h}{\lambda} \, \psi(r),\end{eqnarray}$$

und daraus folgt \( p=\frac{h}{\lambda}=\hbar k \). Der Impuls hängt also über das Planksche Wirkungsquantum mit der Wellenlänge der Schwingung zusammen.

Allgemein lässt sich jede Wellenfunktion \(\psi\), die periodisch in \(k=2\pi / \lambda\) ist, als Fourierreihe schreiben. Dabei gilt \(n\lambda = 2\pi R\) mit \(n\) aus \( \mathbb{N}_0\) und \(R\) als  gedachtem Kreisradius auf dem sich die Schwingung ausbildet.

Die Fourierreihe lautet:

$$\begin{eqnarray} \psi(x) & = & \sum_{n = 0}^\infty A_n \cos \left( \frac{n}{R} x\right) + B_n \sin \left( \frac{n}{R}x \right) \\ & = & \sum_{n = -\infty}^\infty \psi_n e^{inx/R} \\ & = &\sum_{n = -\infty}^\infty \psi_n e^{ikx}\end{eqnarray}$$

Dabei haben wir uns der Eulerschen Formel \( e^{ikx} = \cos(kx) + i \sin(kx) \) bedient, und die Koeffizienten \(\psi_n\) als \(\psi_n =A_n + i B_n\) definiert. Um die Koeffizienten \(\psi_n\) zu bestimmen multiplizieren wir beide Seiten mit \(e^{-im x/R}\) . Dadurch ergibt sich:

$$e^{-i m x/R}\psi(x) = \sum_{n = -\infty}^\infty \psi_n e^{i(n-m)x/R}$$

Als Nächstes  integrieren wir beide Seiten dieser Gleichung über den Umfang eines Kreises.

$$\int\limits_{-\pi R}^{\pi R} \mathrm{d}x~e^{-i m x/R}\psi(x) = \sum_{n = -\infty}^\infty \psi_n \int\limits_{-\pi R}^{\pi R} \mathrm{d}x~ e^{i(n-m)x/R}$$

Der kompliziert erscheinende Ausdruck lässt sich deutlich vereinfachen, denn es gilt

$$\int\limits_{-\pi R}^{\pi R} \mathrm{d}x~ e^{i(n-m)x/R}=\begin{cases}\displaystyle 2\pi R \quad ;n = m\\ \displaystyle 0 \qquad \ ;n\neq m\end{cases}$$

und daraus folgt

$$\int\limits_{-\pi R}^{\pi R} \mathrm{d}x~e^{-i m x/R}\psi(x)  = 2\pi R \psi_m.$$

Die Koeffizienten \(\psi_m\) ergeben sich nach Umstellen zu

$$\psi_n = \frac{1}{2\pi R}\int\limits_{-\pi R}^{\pi R} \mathrm{d}x~e^{-i n x/R}\psi(x)\quad (*)$$

Für \(R\rightarrow\infty\) lässt sich dieser Zusammenhang auf nicht periodische \(\psi(x)\) verallgemeinern. Zunächst führen wir \(\Delta k=\frac{1}{R}\) und eine 'Dichtefunktion' \(\hat{\psi}(k) = \sqrt{2\pi}R\psi_n\) ein. Dann gilt

$$\begin{eqnarray}\psi(x) & = & \sum_{n=-\infty}^\infty \psi_n e^{ikx}\\ & = &\sum_{k=n / R} \Delta k \cdot \frac{1}{\sqrt{2\pi}}\hat{\psi}(k)e^{ikx}\end{eqnarray}$$

und im Grenzfall \(R\rightarrow\infty\) ergibt sich

$$\psi(x) = \frac{1}{\sqrt{2\pi}} \int\limits_{-\infty}^\infty\mathrm{d}k~ e^{ikx} \hat{\psi}(k) \quad (1)$$

Die 'Koeffizientenfunktion' \(\hat{\psi}(k)\) ergibt sich unter Zuhilfenahme von (*) im Grenzfall  \(R\rightarrow\infty\) zu

$$\begin{eqnarray}\hat{\psi}(k) & = & \sqrt{2\pi}R \; \psi_n \\ & = & \sqrt{2\pi}R\frac{1}{2\pi R}\int\limits_{-\infty}^\infty \mathrm{d}x~ e^{-ikx}\psi(x)\ \\ & = & \frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^\infty \mathrm{d}x~ e^{-ikx}\psi(x) \quad (2) \end{eqnarray}$$ 

Die Gleichungen (1) und (2) bilden die weithin bekannte Fouriertransformation. Durch die Symmetrie der Gleichungen wird die Äquivalenz der Darstellung einer Wellenfunktion \(\psi (x)\) im Ortsraum und der Darstellung \(\psi (k)\) im Impulsraum deutlich. 

Inwiefern darin auch die Unschärfebeziehung zwischen Impuls und Ort zum Ausdruck kommt, möchte ich anhand eines Beispiels zeigen. Die Wellenfunktion \(\psi(x)\) eines Teilchens sei wie folgt von -a bis +a um einen Nullpunkt lokalisiert.

Da in der Quantenmechanik die Aufenthaltswahrscheinlichkeit eines Teilchens durch \(\psi^2\) gegeben ist, soll die Amplitude der Wellenfunktion \(1/\sqrt{2a}\) betragen. Dies gewährleistet, dass das Quadrat der Wellenfunktion eine Fläche von 1 zwischen  -a und +a einschließt.

Setzt man \(\psi(x)\) in Gleichung (1) ein, dann ergibt sich:

$$\hat{\psi}(k) = \frac{1}{\sqrt{2\pi}} \frac{1}{\sqrt{2a}}\int\limits_{-a}^a \mathrm{d}x~e^{-ikx} = \frac{1}{\sqrt{\pi a}} \frac{1}{-ik}(e^{-ika} - e^{ika} ) =  \frac{1}{\sqrt{\pi a} }\frac{\sin (ka)}{k}$$

Grafisch lässt sich das als gedämpfte Schwingung mit einem Maximum am Ort 0 darstellen.

Lokalisieren wir das Teilchen im Ortsraum weniger stark und vergrößern a, dann verkleinert sich die Periode der Schwingung, und es bildet sich im k-Raum ein starker Peak um 0 aus. Ein sehr kleines Intervall von k Werten kann also mit hoher Wahrscheinlichkeit beobachtet werden.

Wenn wir das Teilchen im Ortsraum sehr stark lokalisieren und a sehr klein machen, wird \(\psi(k)\) weitgehend unabhängig von k. Für die Sinusfunktion gilt die sogenannte  Kleinwinkelnäherung. Es gilt:

$$\hat{\psi}(k)=\frac{1}{\sqrt{\pi a}}\frac{ka}{k}=\sqrt{a}{\pi}$$

Die Wellenfunktion ist über den gesamten k-Raum verschmiert. Es lässt sich kein einzelner k-Wert finden, der sich durch ein großes \(\psi(k)\) auszeichnet. Der Impuls ist unbestimmbar.

Der Vollständigkeit halber sei zum Abschluss noch erwähnt, dass das Prinzip der Unschärfe natürlich für alle physikalischen Größen und deren Fouriertransformierte gilt. Die Unschärfebeziehung ist also keine Eigenschaft der Quantenmechanik sondern findet sich auch in anderen Bereichen wie z.B der Signalverarbeitung wieder. Bedingt durch die Unschärfebeziehung ist es z.B. nicht möglich das genaue Frequenzspektrum eines Signals zu einem bestimmten Zeitpunkt zu messen.

Die Formeln in diesem Beitrag wurden übrigens mit MathJax erstellt. Es genügt ein kleines 'Code-Snippet' als 'Gadget' ins Layout des blogs einzufügen, und schon werden die TeX Befehle, die man in den Editor schreibt, in der Webansicht wohl gesetzt dargestellt.

Viel Spaß beim selber rechnen, schreiben und nachdenken...