NumPy

Was ist NumPy?

Python ist in technisch-wissenschaftlichen Anwendungen eine der beliebtesten Programmiersprachen und wird von vielen Nicht-Informatikern eingesetzt. Diese Beliebtheit liegt neben der einfachen Sprachsyntax auch an den leistungsfähigen Bibliotheken, die es für mathematische Anwendungen in Python gibt.

Eine wichtige Bibliothek ist NumPy. Sie stellt viele Funktionen und Klassen bereit, die man für Berechnungen in den Natur- und Ingenieurswissenschaften benötigt.

NumPy ist Teil einer Sammlung von Python-Modulen für das wissenschaftliche Rechnen, die SciPy heißt.

Ein großer Vorteil von NumPy liegt darin, dass es die eigentlichen Berechnungen nicht in Python durchführt, sondern hinter den Kulissen auf sehr schnelle C/C++-Funktionen zurückgreift, die alle Beschleunigungsmöglichkeiten moderner Prozessoren voll ausnutzen. Deswegen laufen Berechnungen mit NumPy auch bei großen Datenmengen schnell ab.

NumPy

• Paket für wissenschaftliche Berechnungen mit Python
• n-dimensionale Arrays (Vektoren, Matrizen)
• Lineare Algebra
• Fourier Transformationen
• Zufallszahlen
• Statistik
• Open Source
• …

Sie lernen in dieser Vorlesung NumPy kennen, weil es Ihnen bei der täglichen Arbeit mit den mathematischen Herausforderungen einer Ingenieurswissenschaft helfen kann.

NumPy importieren

Das zentrale Element von NumPy ist das ndarray, ein mehrdimensionales Array von Daten, das anders als die Liste in Python eine feste Größe hat und nur Daten eines Datentyps aufnimmt. Die Daten im Array können verändert werden, nicht aber die Größe und die Dimensionen.

Bevor man NumPy nutzen kann, muss man es installieren (siehe hierzu die Dokumentation) und dann in das Python-Programm importieren.

• NumPy muss installiert werden
• NumPy muss mit import numpy in eigene Skripte importiert werden

import numpy as np

# Array mit Ganzzahlen anlegen
a = np.array([[1, 2, 3], [4, 5, 6]])
print(a)
# [[1 2 3]
# [4 5 6]]

# Array mit Fließkommazahlen anlegen
b = np.array([[1, 2, 3], [4, 5, 6]], float)
print(b)
# [[ 1.  2.  3.]
# [ 4.  5.  6.]]

NumPy mit dem Alias np zu installieren ist eine Konvention, die Sie häufig bei den Verwendern von NumPy sehen werden. Verpflichtend ist das natürlich nicht und Sie können jeden Alias verwenden, der Ihnen gefällt.

Ein Array in NumPy wird mit der Funktion array angelegt. Dieser können Sie normale Python-Listen übergeben, um diese in ein NumPy-Array umzuwandeln. Als optionalen Parameter können Sie den Datentyp des Array bestimmen; NumPy konvertiert dann alle Elemente in diesen Datentyp.

l = [ 1, 2.3, 3 ]
a = np.array(l, int)
print(a) # -> array([1, 2, 3])

b = np.array(l, float)
print(b) # -> array([1., 2.3, 3.])

c = np.array(l, str)
print(c) # -> array(['1', '2.3', '3'], dtype='<U1')

Das Array mit den string-Werten ist natürlich für Berechnungen nur begrenzt sinnvoll und soll hier nur den grundlegenden Mechanismus verdeutlichen.

NumPy Beispiele

NumPy ist so umfangreich, dass es hier nicht komplett behandelt werden kann. Trotzdem sollen einige Beispiele die Bereite der Möglichkeiten zeigen.

In weiteren Verlauf dieser Vorlesung werden wir nur auf die Funktionen zur linearen Algebra in NumPy eingehen.

import numpy as np

# Matrix nur mit Nullen
a = np.zeros((2, 3))
# array([[ 0.,  0.,  0.],
#        [ 0.,  0.,  0.]])

# Form der Matrix
a.shape # -> (2, 3)

# Diagonal-Matrix
b = np.diag([1, 2, 3])
# array([[1, 0, 0],
#        [0, 2, 0],
#        [0, 0, 3]])

Warum steht die (2, 3) im Beispiel in Klammern? Die Funktion zeros erwartet als ersten Parameter die Dimensionen des Arrays. Da beliebig viele Dimensionen unterstützt werden, muss man diese in Form eines Tupels übergeben.

np.zeros((2, 2, 2))
# array([[[0., 0.],
#        [0., 0.]],
#
#       [[0., 0.],
#        [0., 0.]]])

Man kann in NumPy auch sehr einfach Zufallszahlen erzeugen oder Daten speichern und laden.

# Matrix mit Zufallszahlen
a = np.random.random((2,3))
# array([[ 0.1834851 ,  0.49019724,  0.1798613 ],
#        [ 0.96952566,  0.884116  ,  0.40070798]])

# Speichern in eine Datei
np.savetxt("a_out.txt", a)

# Laden aus Datei
b = np.loadtxt("a_out.txt")

Matrizen sind veränderlich

Das ndarray ist die zentrale Datenstruktur von NumPy. Es hat nach der Erzeugung eine feste Größe, erlaubt aber die Veränderung der darin enthaltenen Daten. Damit steht es zwischen den Listen und Tupeln in Python. NumPy-Arrays stellen Matrizen im mathematischen Sinne dar.

Matrizen

• haben feste, unveränderliche Größe (→ Tupel)
• enthalten Daten, die geändert werden können (→ Liste)

import numpy as np

# 2x2-Matrix nur mit Nullen anlegen
a = np.zeros((2, 2))
print(a)
# [[ 0.  0.]
#  [ 0.  0.]]

# Element an der Position (0, 0) ändern
a[0, 0] = 7
print(a)
# [[ 7.  0.]
#  [ 0.  0.]]

Das Beispiel zeigt, wie man eine Null-Matrix anlegt, indem man die gewünschte Dimension als Tupel übergibt. Die Elemente der Matrix kann man dann über deren Index adressieren. Der Index hat genau so viele Dimensionen, wie die Matrix auch.

v = np.zeros(3)
print(v) # -> [0. 0. 0.]

v[2] = 42
print(v) # -> [ 0.  0. 42.]

v[0, 0] = 7.0 # Fehler
# IndexError: too many indices for array

Attribute von Matrizen

Jedes ndarray in NumPy kann bezüglich der Eigenschaften der dargestellten Matrize befragen. Hierzu bietet es eine Reihe von Attributen an.

Attribute von Matrizen

• ndim: Dimensionen
• shape: Form
• size: Anzahl der Elemente
• T: Transponierte Form
• dtype: Datentyp

a = np.array([[1, 2, 3], [4, 5, 6]])

print(a.ndim) # -> 2
print(a.shape) # -> (2, 3)
print(a.size) # -> 6
print(a.dtype) # -> int64
print(a.T)
# [[1 4]
# [2 5]
# [3 6]]

Die transponierte Form vertauscht Zeilen und Spalten:

a = np.array([[1, 2, 3], [4, 5, 6]])
print(a)
# [[1 2 3]
#  [4 5 6]]

print(a.T)
# [[1 4]
#  [2 5]
#  [3 6]]

Operatoren für Arrays

Der erste große Unterschied zwischen NumPy-Arrays und den Listen in Python ist, dass die arithmetischen Operatoren elementweise angewendet werden. Addiert man z. B. zwei Arrays, dann werden die Elemente der Arrays addiert. Multipliziert man ein Array mit einem Skalar, dann werden alle Elemente mit dem Skalar multipliziert.

Operatoren werden elementweise angewandt (+, -, *, /)

• Skalar: auf jedes Element des Arrays
• Array: auf die passenden Elemente der Arrays

a = np.ones((3, 3))

a + 1
# [[ 2.,  2.,  2.],
#  [ 2.,  2.,  2.],
#  [ 2.,  2.,  2.]]

a * 3
# [ 3.,  3.,  3.],
# [ 3.,  3.,  3.],
# [ 3.,  3.,  3.]]

Das Beispiel zeigt, wie die Skalare einfach auf alle Elemente der Matrize angewandt werden.

Sind beide Operanden NumPy-Arrays, erfolgt eine elementweise Anwendung der Operatoren. Die Elemente mit demselben Index in beiden Arrays werden verknüpft und das Ergebnis wieder an diesen Index im Ziel-Array geschrieben.

a = np.arange(4)
# array([0, 1, 2, 3])

b = np.array([2, 3, 2, 4])
# array([2, 3, 2, 4])

b * a # array([0, 3, 4, 12])
b - a # array([2, 2, 0, 1])

c = [2, 3, 4, 5]
a * c
# array([0, 3, 8,15])

In diesem Beispiel wird mit arange(4) ein Array a mit vier Elementen erzeugt, die den Bereich 0 bis 3 abdecken. Das zweite Array b wird über eine Python-Liste erzeugt.

Man sieht, wie die Operatoren auf die einzelnen Elemente der Arrays angewendet werden, um das Ergebnis-Array zu erzeugen.

Was aber passiert, wenn die Arrays unterschiedliche Größen haben? Dann fehlen ja bei einem Operanden die Elemente für die Anwendung des Operators? In diesem Fall kommt das Broadcasting zum Einsatz.

Broadcasting

Das Broadcasting kümmert sich darum, wie die Operationen ablaufen sollen, wenn die Arrays nicht dieselbe Dimension haben.

• Array broadcasting: Wenn NumPy auf zwei Arrays arbeiten soll, werden die Dimensionen verglichen
• die Dimensionen sind kompatibel, wenn
• sie dieselbe Größe haben oder
• eine von ihnen 1 ist

a = np.arange(3)
b = np.arange(4)

a + b
# ValueError: operands could not be broadcast together with shapes (3,) (4,)

In diesem Beispiel kann die Addition nicht durchgeführt werden, weil das Array a ein Element zu wenig hat. Ein Broadcasting ist nicht möglich, weil die Arrays nur eine Dimension haben und diese nicht identisch und auch nicht 1 ist.

a = np.array([ 0, 1, 2 ])
b = np.array([3])

c = a + b
print(c) # -> [3 4 5]

In diesem Fall funktioniert das Broadcasting, weil das Array b die Dimension 1 hat und somit sein Element (3) mit allen Elementen des Arrays a verknüpft werden kann.

Broadcasting führt also keine "magischen" Modifikationen an den Arrays durch, sorgt aber dafür, dass in bestimmten Fällen Python-Code vermieden wird.

Ohne Broadcasting müsste im Beispiel das Element aus dem Array b herausgeholt werden und als Skalar verwendet:

a = np.array([ 0, 1, 2 ])
b = np.array([3])

c = a + b[0]
print(c) # -> [3 4 5]

Was hier bei einem eindimensionalen Array noch relativ wenig Aufwand macht, kann bei mehrdimensionalen Arrays erhebliche Schreibarbeit erfordern – diese wird durch das Broadcasting minimiert.

Vektor-Operationen

NumPy bietet die üblichen Vektor-Operationen an. Diese sind allerdings nicht als Methoden auf dem NumPy-Array realisiert, sondern sind eigene Funktionen, denen man die Vektoren übergeben kann. Der Grund für diese Lösung liegt darin, dass Vektoren kein Sonderfall von Matrizen sind.

Vektor-Operatonen

• inner: Inneres Produkt
• outer: Äußeres Produkt
• dot: Matrix-Multiplikation
• cross: Kreuzprodukt (Vektorprodukt)

# Listen werden automatisch in Vektoren konvertiert
u = [1, 2, 3]
v = [1, 1, 1]

np.inner(u, v) # -> 6
np.dot(u, v) # -> 6
np.cross(u, v) # -> [-1,  2, -1]

np.outer(u, v) # -> [[1, 1, 1], [2, 2, 2], [3, 3, 3]]

Man muss keine Python-Listen verwenden, sondern kann auch eindimensionale NumPy-Arrays verwenden. Die Möglichkeit, direkt Listen einzusetzen ist als Vereinfachung gedacht, um die Konvertierung in Arrays zu vermeiden.

# Listen werden automatisch in Vektoren konvertiert
u = np.array([1, 2, 3])
v = np.array([1, 1, 1])

np.inner(u, v) # -> 6
np.dot(u, v) # -> 6
np.cross(u, v) # -> [-1,  2, -1]

Slicing

NumPy erweitert die Syntax für das Slicen von Liste auf die NumPy-Arrays, sodass man auch diese komfortable in Teilstücke zerlegen kann.

a = np.array([[1, 2, 3], [4, 5, 6]])

# 2 Zeile, alle Spalten
a[1, :] # -> [4, 5, 6]

# 1. und 2. Zeile, alle Spalten
a[0:2]
# [[1, 2, 3],
# [4, 5, 6]

# Alle Zeilen, Spalten 3 und 4
a[:, 2:4]
# [[3],
#  [6]]

Iterieren über Matrizen

So wie man über Listen in Python iterieren kann – sich also alle Elemente nacheinander ansehen – kann man auch über die Arrays von NumPy laufen. Diese sind so konstruiert, dass man die normale for-Schleife aus Python für die Iteration einsetzen kann.

a = np.array([[1, 2, 3], [4, 5, 6]])

# Zeilen und Spaltenweise iterieren
for row in a:
    for column in row:
        print(column)

# Über alle Elemente iterieren
for element in a.flat:
    print(element)

In beiden Fällen ist die Ausgabe "1 2 3 4 5 6". Allerdings weiß man im zweiten Fall nicht mehr, wann der Wechsel in die nächste Zeile erfolgt, da diese Informationen durch das .flat verloren geht.

Matrix aus Datei lesen

Größere Matrizen möchte man möglicherweise nicht im Python-Programm selbst pflegen, sondern in einer externen Datei ablegen. Oder die Daten, die man auswerten möchte stammen aus einer anderen Software, z. B. aus einem Laborsystem, mit dem man Messungen aufgezeichnet hat.

Für diese Fälle bietet NumPy die Möglichkeit, die Daten eines Arrays aus einer externen Datei zu laden.

• loadtxt erlaubt es ein Array aus einer Datei zu lesen

1 2 3.5
4 5 6
7.0 8 9.3

import numpy as np

m = np.loadtxt("matrix.txt")

print(m)
# [[ 1. ,  2. ,  3.5],
#  [ 4. ,  5. ,  6. ],
#  [ 7. ,  8. ,  9.3]])

Die Datei muss hierbei das Format haben, dass die Zeilen des Arrays auch Zeilen in der Datei sind. Jede Zeile wird durch ein Newline-Zeichen (\n) abgeschlossen.

Lineare Regression

Die folgenden Beispiele zur linearen Regression verwenden die Zeichenfunktionen aus matplotlib, die erst im nächsten Kapitel besprochen werden. Da die lineare Regression aber so besser darzustellen ist, sei dieser Vorgriff entschuldigt.

import numpy as np
import matplotlib.pyplot as plt

x = [  0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0 ]
y = [ 18.7, 20.0, 20.3, 21.5, 22.0, 23.0, 23.0, 25.5, 24.0 ]

fit = np.polyfit(x, y, 1)
fit_fn = np.poly1d(fit)

print(fit_fn) # => 0.7433 x + 19.03

plt.plot(x, y, 'o', x, fit_fn(x), 'k')
plt.xlim(0, 9)
plt.ylim(18, 25)
plt.show()

Die Magie liegt hier in der Funktion polyfit, die ausgehend von den x- und y-Werten ein Polynom bestimmt, dass die Daten mit den geringsten Fehlerquadraten annähert. Der letzte Parameter von polyfit bestimmt den Grad des Polynoms, der hier 1 ist. Damit wird eine lineare Funktion bestimmt, welche die Bedingung erfüllt und somit eine lineare Regression durchgeführt.

Mit np.poly1d(fit) wird aus der Beschreibung des Polynoms eine Funktion erzeugt, die dann benutzter werden kann, um die Gerade der linearen Regression zu zeichnen.

Im folgenden Beispiel sieht man das deutlicher:

x = [  0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0 ]
y = [ 18.7, 20.0, 20.3, 21.5, 22.0, 23.0, 23.0, 25.5, 24.0 ]

fit = np.polyfit(x, y, 1)
fit_fn = np.poly1d(fit)

print(fit_fn(0)) # -> 19.026666666666667
print(fit_fn(1)) # -> 19.77
print(fit_fn(2)) # -> 20.513333333333335
print(fit_fn(3)) # -> 21.256666666666668
print(fit_fn(8)) # -> 24.973333333333336

Regression höherer Ordnung

Nachdem die Funktionsweise von polyfit klar geworden ist, kann man leicht auch Regressionen höherer Ordnung bestimmen. Hier im Beispiel ein Polynom mit Grad 2.

import numpy as np
import matplotlib.pyplot as plt

x = [ 0.0, 1.0, 2.0, 3.0, 4.0,  5.0,  6.0,  7.0,  8.0,  9.0, 10.0 ]
y = [ 0.0, 5.1, 21.6, 40.1, 80.4, 120.6, 185.5, 240.3, 313.9, 397.3, 490.5 ]

fit = np.polyfit(x, y, 2)
fit_fn = np.poly1d(fit)

print(fit_fn)

plt.plot(x, y, 'o', x, fit_fn(x), 'k')
plt.xlim(0, 11)
plt.ylim(0, 510)
plt.show()