Die hier vorgestellten Bibliotheken erlauben, dass eine Variable als Zahlenwert mehr als nur einen einzelnen Skalar beinhaltet. Sie bieten eine hocheffiziente Verwaltung und Verarbeitung eines Arrays von Werten gleichen (homogoenen) Typs. Damit wird dem Schritt zu Vektoren, Matrizen und Tensoren Tür und Tor geöffnet, um numerische Berechnungen effizient ausdrücken zu können.
Die Programmbibliothek NumPy
liefert mit numpy.ndarray
einen Rang-n Tensor Typ,
der effizientes Rechnen mit Vektoren und Matritzen erlaubt.
All diese Berechnungen auf Basis dieses np.ndarray
ist gemeinsam,
dass sie dank der zugrundeliegenden Vektorisierung sehr effizienter ablaufen.
Für performanten (und häufig auch besser lesbaren) Code ist es daher immer notwendig,
auf diese Operationen zurückzugreifen -- auch wenn sie anfangs etwas verwirrend sein mögen :-)
Lit.: Stefan Van Der Walt, S. Chris Colbert, Gaël Varoquaux: The NumPy array: a structure for efficient numerical computation, Computing in Science and Engineering 13, 2 (2011) 22-3, PDF
Die Bibliothek numpy
wird üblicherweise immer mit dem Kurznamen Alias np
importiert:
import numpy as np
Beispiel: Die Variablen $v$ und $w$ sollen einen Vektor aus $\mathbb{R}^3$ repräsentieren.
Die Addition und Multiplikation werden elementweise durchgeführt,
während .dot()
das Skalarprodukt (engl. dot product) durchführt.
Jetzt definieren wir $v$ und $w$, wobei $v = (4,\,3,\,-1)$ sein soll und $w$ ein zufällig gewählter Vektor. Dann führen wir die Berechnungen durch.
v = np.array([4,3,-1]) # Konstruktion aus einer Python "list"
w = np.random.rand(3) # 3 zufällige uniform verteilte Zahlen
v
array([ 4, 3, -1])
w
array([ 0.70613368, 0.61589773, 0.31350427])
v + w
array([ 4.70613368, 3.61589773, -0.68649573])
v - w
array([ 3.29386632, 2.38410227, -1.31350427])
v * w
array([ 2.82453471, 1.84769319, -0.31350427])
Das (generalisierte) Skalarprodukt von $v$ und $w$
v.dot(w)
4.358723634005492
Der Typ von v
ist das schon erwähnte ndarray
type(v)
numpy.ndarray
Die Typen der Zahlen in den jeweiligen ndarray
Objekten ist homogen und kann mit .dtype
abgerufen werden.
Standardmäßig sind dies int64
oder float64
.
Insbesondere bedeutet dies, dass es sich nicht um native Python-Integer Werte handelt und es daher einen Überlauf geben kann. siehe dtype.
Der große Vorteil ist, dass dadurch vektorisierte Berechnungen möglich sind und mit einer viel höheren Geschwindigkeit verarbeitet werden können.
v.dtype
dtype('int64')
w.dtype
dtype('float64')
Bei der Instanzierung eines ndarray
kann der dtype
explizit angegeben werden -- es kann hier z.B. der Python object
-Typ gewält werden um mit int
-Werten zu rechnen:
python_ints = np.array([11111111111111111131111111,
14441111111111111111111111],
dtype=object)
python_ints.dtype
dtype('O')
python_ints.sum()
25552222222222222242222222
Neben diesen Grundoperationen, werden alle NumPy-spezifischen Funktionen automatisch über tensor-wertige Argumente ausgeführt. Das heißt beispielsweise, dass die $\sin$-Funktion über einen Vektor $x := (1,2,3)$ ergibt wieder einen Vektor, dessen Einträge die numerischen Werte von $(\sin(1),\,\sin(2),\,\sin(3))$ beinhalten.
x = np.array([1,2,3])
np.sin(x)
array([ 0.84147098, 0.90929743, 0.14112001])
$(e^1,\,e^2,\,e^3)$
np.exp(x)
array([ 2.71828183, 7.3890561 , 20.08553692])
Verschachtelte Listen werden zu einer n x m Matrix instanziert.
m = np.array([[ 1,-2, 3, 8],
[ 4, 6, 9,-9],
[-1, 0, 5, 2],
[11,13,15,16]])
m
array([[ 1, -2, 3, 8], [ 4, 6, 9, -9], [-1, 0, 5, 2], [11, 13, 15, 16]])
Berechnung des Cosinus über alle Einträge der Matrix m
:
np.cos(m)
array([[ 0.54030231, -0.41614684, -0.9899925 , -0.14550003], [-0.65364362, 0.96017029, -0.91113026, -0.91113026], [ 0.54030231, 1. , 0.28366219, -0.41614684], [ 0.0044257 , 0.90744678, -0.75968791, -0.95765948]])
Der Doppelpunkt :
steht hierbei für die gesamte Spalte, n:m
für ein Intervall, wobei n
oder m
weggelassen werden können. Negative Werte werden vom Ende weg gezählt.
Dies ist eine Verallgemeinerung der Listen-Indizierung.
m
array([[ 1, -2, 3, 8], [ 4, 6, 9, -9], [-1, 0, 5, 2], [11, 13, 15, 16]])
Element an Position (2, 2)
m[2, 2]
5
Setzen des Wertes an einer Position:
m[2, 2] = 987
m
array([[ 1, -2, 3, 8], [ 4, 6, 9, -9], [ -1, 0, 987, 2], [ 11, 13, 15, 16]])
Spalte 0
m[:, 0]
array([ 1, 4, -1, 11])
Zeile 1
m[1, :]
array([ 4, 6, 9, -9])
Letzte Spalte
m[:, -1]
array([ 8, -9, 2, 16])
Auswahl von zwei Elementen: Zeile 2, Spalte 0 und 2
m[2, [0, 2]]
array([ -1, 987])
m[-1:]
array([[11, 13, 15, 16]])
"Teilmatrix, vom Anfang bis zur vorletzten Zeile und vor-vorletzten Spalte"
m[:-1, :-2]
array([[ 1, -2], [ 4, 6], [-1, 0]])
Die .diagonal()
Funktion gibt die Diagonalelemente zurück.
m.diagonal()
array([ 1, 6, 987, 16])
Ein doppelter Doppelpunkt ::
gibt für ein n:m
Intervall zusätzlich eine Schrittweite an.
m[::2, :]
array([[ 1, -2, 3, 8], [ -1, 0, 987, 2]])
Die Schrittweite kann auch negativ sein, z.B. zum Umdrehen der Richtung:
m[:, ::-1]
array([[ 8, 3, -2, 1], [ -9, 9, 6, 4], [ 2, 987, 0, -1], [ 16, 15, 13, 11]])
Eine wichtige Eigenschaft dieser Extraktionen ist, dass sie keine Kopie der Daten liefern. Dies ist im Allgemeinen schlecht, da es Zeit kostet und unnötigen Speicherplatz verbraucht.
Daher wirken sich Änderungen an diesen "Ansichten" auf die Originaldaten aus -- was das folgende Beispiel verdeutlicht.
Möchte man tatsächlich eine Kopie haben, so gibt es die .copy()
Methode.
orig = np.array([[1, 2], [3, 4]])
orig
array([[1, 2], [3, 4]])
# Kopie
copy = orig.copy()
# View der zweiten Spalte, als 1-dim Array
view = orig[:, 1]
view
array([2, 4])
Zweites Element
view[1]
4
Zuweisen des Wertes 99
an das zweite Element des Vektors view
.
view[1] = 99
view
array([ 2, 99])
orig
hat ebenfalls (!) die 99
in den Daten.
orig
array([[ 1, 2], [ 3, 99]])
... copy aber nicht:
copy
array([[1, 2], [3, 4]])
m.T
array([[ 1, 4, -1, 11], [ -2, 6, 0, 13], [ 3, 9, 987, 15], [ 8, -9, 2, 16]])
m[:, 0].dot(m[:, 1])
165
.ravel()
gibt alle der Matrix oder dem Tensor zugrundeliegende Werte zurück
(je nach dem in welcher Reihenfolge die Werte gespeichert werden, also spaltenweise oder zeilenweise)
m.ravel()
array([ 1, -2, 3, 8, 4, 6, 9, -9, -1, 0, 987, 2, 11, 13, 15, 16])
Darüber hinaus gibt es einige hilfreiche Funktionen, um einen Raster für das Auswerten einer Funktion an Punkten zu ermöglichen.
np.arange
ist ähnlich wie Python's range
Funktion,
liefert jedoch gleich np.ndarray
Vektoren zurück.
np.arange(1, 5, .5)
array([ 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5])
10 gleichmäßig verteilte Punkte von 0 to 2
np.linspace(0, 2, 10)
array([ 0. , 0.22222222, 0.44444444, 0.66666667, 0.88888889, 1.11111111, 1.33333333, 1.55555556, 1.77777778, 2. ])
10 logarithmisch verteilte Punkte von $\text{base}^{-1}$ nach $\text{base}^{2}$ mit base = 10
np.logspace(-1, 2, 10, base=10)
array([ 0.1 , 0.21544347, 0.46415888, 1. , 2.15443469, 4.64158883, 10. , 21.5443469 , 46.41588834, 100. ])
Zur Auswertung von mehrstelligen Funktionen (also über mehrere Achsen), kann ein entsprechendes Grid von Punkten aus zwei oder mehr Vektoren konstruiert werden.
xx, yy = np.meshgrid([1, 2, 3], [0, 10, 20])
xx
listet alle x-Werte auf, die jeweils mit den ...
xx
array([[1, 2, 3], [1, 2, 3], [1, 2, 3]])
... entsprechenden y-Werten in yy
korrespondieren:
yy
array([[ 0, 0, 0], [10, 10, 10], [20, 20, 20]])
Also: (1, 0), (2, 0), (3, 0), (1, 10), (2, 10), ...
Dies ist zur Auswertung von mehrstellige Funktionen sehr praktisch:
1 + xx + 2 * yy
array([[ 2, 3, 4], [22, 23, 24], [42, 43, 44]])
Eine Indexmatrix ist ähnlich, liefert aber unmittelbar die (x, y)-Index Einträge der Positionen.
xx, yy = np.indices((3, 2))
xx
array([[0, 0], [1, 1], [2, 2]])
yy
array([[0, 1], [0, 1], [0, 1]])
Die Dimensionierung eines Vektors oder Matrix lässt sich jederzeit ändern.
Das liegt daran, dass ein np.ndarray
aus einem Vektor aller Werte und einer Information,
wie lang die Zeilen bzw. Spalten sind besteht.
Diese Shape-Information kann jederzeit geändert werden und damit ändert sich auch die Dimensionierung des Arrays -- es muss jedoch die Anzahl der Elemente gleich bleiben!
1d-Vektor der Länge 25
c = np.arange(25)
c.shape
(25,)
Ändern der Rang-1 shape (25,)
auf Rang-2 (5, 5)
, also eine 5x5-Matrix mit 25 Elementen:
c.shape = (5, 5)
c
array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14], [15, 16, 17, 18, 19], [20, 21, 22, 23, 24]])
oder <array>.reshape((tupel))
, wobei -1
verwendet werden kann um die Dimensionierung automatisch anzupassen.
# Vektor mit 8 Elementen
c = np.arange(8)
# Umwandlung zu einer Matrix mit 4 Spalten
c = c.reshape((-1, 4))
c
array([[0, 1, 2, 3], [4, 5, 6, 7]])
c.shape
(2, 4)
c.ndim
2
Rang 3, mit 2x2x2 Elementen:
c.shape = (2, 2, 2)
c
array([[[0, 1], [2, 3]], [[4, 5], [6, 7]]])
c.ndim
3
Tipp: Instanzierung und reshaping kann leicht kombiniert werden:
np.linspace(-2, 2, 9).reshape((-1, 3))
array([[-2. , -1.5, -1. ], [-0.5, 0. , 0.5], [ 1. , 1.5, 2. ]])
Manchmal bietet sich aus rechnerischen Gründen an, einen Vektor zu einer n x 1
-Matrix zu verwandeln.
Dies -- und andere Redimensionierungen, die eine neue Achse einfügen -- können mit np.newaxis
erledigt werden:
f = np.linspace(0, 1, 5)
f.shape
(5,)
f2 = f[:, np.newaxis]
f2.shape
(5, 1)
f2
array([[ 0. ], [ 0.25], [ 0.5 ], [ 0.75], [ 1. ]])
f2.shape
(5, 1)
Bestehenden Vektoren und Matritzen können mittels r_
und c_
Operatoren zusammengefügt werden.
a = np.arange(10)
b = np.arange(100, 110)
Zeilenweise mit r_
(für "row")
np.r_[a, b]
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109])
Spaltenweise mit c_
(für "column")
np.c_[a, b, b]
array([[ 0, 100, 100], [ 1, 101, 101], [ 2, 102, 102], [ 3, 103, 103], [ 4, 104, 104], [ 5, 105, 105], [ 6, 106, 106], [ 7, 107, 107], [ 8, 108, 108], [ 9, 109, 109]])
Eine der wichtigsten (und komplexeren) Techniken ist "Broadcasting". Hierbei werden Skalare, Vektoren und Matritzen in Rechenoperationen gemischt. Das funktioniert so, dass im ersten Rang, der einem der Argumente fehlt, alle Werte auf die Art wiederholt werden, sodass die Dimensionierung in diesem Rang passt.
Vereinfachtes Beispiel: Eine Addition des Vektors $[1,\,2,\,5]$ mit dem Skalar $10$ führt dazu, dass der Skalar intern wie ein $[10,\,10,\,10]$ Vektor behandelt wird und daher das Ergebnis $[11,\,12,\,15]$ ist.
np.array([1,2,5]) + 10
array([11, 12, 15])
oder Matrix + Vektor
a = np.arange(16).reshape((4, 4))
a
array([[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11], [12, 13, 14, 15]])
b = np.array([-1, 0, 10, 100])
b
array([ -1, 0, 10, 100])
a + b
array([[ -1, 1, 12, 103], [ 3, 5, 16, 107], [ 7, 9, 20, 111], [ 11, 13, 24, 115]])
Setzen einer ganzen Zeile:
m[-1, :] = 777
m
array([[ 1, -2, 3, 8], [ 4, 6, 9, -9], [ -1, 0, 987, 2], [777, 777, 777, 777]])
Bei der vorhergehenden Addition wird der "horizontale" b
-Vektor entlang der ersten Dimension der Matrix a
(also zeilenweise) vervielfacht.
In Kombination mit np.newaxis
lassen sich "äußere" Operationen durchführen:
Hier wird ein Vektor zu einer 4x1-Matrix verändert und dann addiert.
a = np.arange(4)
a
array([0, 1, 2, 3])
b = a[:, np.newaxis]
b
array([[0], [1], [2], [3]])
a + b
array([[0, 1, 2, 3], [1, 2, 3, 4], [2, 3, 4, 5], [3, 4, 5, 6]])
Hier werden also, damit es sich "ausgeht", beide Vektoren a
und b
entlang der jeweils anderen Achse vervielfacht.
Multiplikation eines Tensors 3-ter Stufe mit einem Vektor
a = np.arange(8)
a.shape = (2,2,2)
a
array([[[0, 1], [2, 3]], [[4, 5], [6, 7]]])
b = np.array([1,11])
b
array([ 1, 11])
a * b
array([[[ 0, 11], [ 2, 33]], [[ 4, 55], [ 6, 77]]])
... und deren (generalisiertes) Skalarprodukt ergibt eine 2x2-Matrix:
a.dot(b)
array([[11, 35], [59, 83]])
In der Praxis werden häufig Funktionen auf Vektoren oder Matrizen angewendet. Siehe auch Numpy Routines
m = np.array([[5, 6, -1, 5],
[1, 1, 1, 1],
[0,-5, 10, 2]])
m
array([[ 5, 6, -1, 5], [ 1, 1, 1, 1], [ 0, -5, 10, 2]])
m.sum()
26
m.sum(axis=0)
array([ 6, 2, 10, 8])
m.sum(axis=1)
array([15, 4, 7])
m.min()
-5
m.ptp
steht für "peak-to-peak", gibt also den Wertebereich an.
Hier wird außerdem mit axis=0
nur entlang der Zeilen gearbeitet.
m.ptp(axis=0)
array([ 5, 11, 11, 4])
Die Submodule numpy.linalg und scipy.linalg beinhalten diverse Routinen für lineare Algebra.
import numpy.linalg as LA
Definition der Matrix $m$
m = np.array([[ 4, -1.1],
[-2, 8 ]])
Eigenwerte & Eigenvektoren vom $m$
eigenvalues, eigenvectors = LA.eig(m)
print(eigenvalues)
print(eigenvectors)
[ 3.51002008 8.48997992] [[-0.91347498 0.23795303] [-0.40689491 -0.97127666]]
Testen, ob der 0
-te Eigenwert/Vektor passt: beide Vektoren müssen bis auf Rundungsfehler gleich sein (das kontrolliert np.allclose
)
np.allclose(
m.dot(eigenvectors[:,0]),
eigenvalues[0] * eigenvectors[:,0]
)
True
Determinante von $m$
LA.det(m)
29.799999999999997
Invertieren der Matrix $m$
LA.inv(m)
array([[ 0.26845638, 0.03691275], [ 0.06711409, 0.13422819]])
LA.inv(m).dot(m)
array([[ 1., 0.], [ 0., 1.]])
Ein wichtiger Zweig der Numerischen Mathematik ist Optimierung. Diese besteht aus einer Sammlung standardisierter Problemstellungen, welche unter eventuell gegebenen Nebenbedingungen einen Punkt $x \in \mathbb{R}^n$ suchen, wo eine gegebene Zielfunktion $f: \mathbb{R}^n \rightarrow \mathbb{R}$ einen minimalen oder maximalen Wert annimmt.
Problemstellungen dieser Art finden sich in mannigfaltigen Anwendungsgebieten und bringen diverse Varianten dieser Fragestellung hervor. Die folgenden Beispiele zeigen exemplarisch einige Varianten vor.
Link: scipy.optimize
from scipy.optimize import minimize
$f_1(x): \mathbb{R}^3 \rightarrow \mathbb{R}$ wird optimiert.
f1 = lambda x: (4 * x[0] - x[1])**2 + (x[1] - x[2])**2 * x[0] - x[2]
$x_0$ ist der Startwert in $\mathbb{R}^3$
x0 = [1., 2., 0.]
bounds
ist die Box der Randbedingungen an $x$:
bounds = [(-10, 10), (-10, 10), (0, 1)]
opti1 = minimize(f1, x0, bounds = bounds)
opti1
fun: -0.99999999999996092 hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64> jac: array([ 1.52100554e-06, -4.32986980e-07, -9.99999905e-01]) message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' nfev: 64 nit: 13 status: 0 success: True x: array([ 0.24999999, 0.9999998 , 1. ])
x_opti = opti1["x"]
x_opti
array([ 0.24999999, 0.9999998 , 1. ])
Eine höherdimensionale Funktion kann mittels der Standardalgorithmen "BFGS" bzw. "L-BFGS-B" optimiert werden. SciPy liefert hier eine moderne und "offizielle" Fortranimplementation dieser Verfahren aus.
Hier eine Formulierung für ein gesuchtes $x$, welches in jeder Dimension beschränkt ist (d.h. "in einer Box liegt").
$$\begin{align} \min & \quad f(x) && \quad f:\mathbb{R}^n \rightarrow \mathbb{R}\\ \text{s.t.} & \quad x_{lb} \leq x \leq x_{ub} && \quad x,\,x_{lb},\,x_{ub} \in \mathbb{R}^n \end{align}$$Das Ziel ist das $x$ in dieser "Box" zu finden, welches $f(x)$ minimiert.
from scipy.optimize import fmin_l_bfgs_b
Mit diesen (wichtigsten) Argumenten:
func : callable f(x): Function to minimise.
x0 : ndarray : Initial guess.
approx_grad : bool : Whether to approximate the gradient numerically.
bounds : list : ``(min, max)`` pairs for each element in ``x``
defining the bounds on that parameter.
Beispiel:
f1
, welches mit dem Variablenvektor aufgerufen wird und den Wert zurückgeben soll.[0, 0]
wenn in $\mathbb{R}^2$approx_grad = True
, wenn die Zielfunktion keinen Gradienten zurückgibtbounds = [(-5, 5), (-5, 5)]
f1 = lambda x : x[0]**3 + (1 + x[1]**2) - (2 - x[0]) * x[1]
x_opt, f_x_opt, info = fmin_l_bfgs_b(f1,
x0=[0, 0],
bounds=[(-5, 5), (-5, 5)],
approx_grad=True)
"Minimum: f(%s) = %.6f" % (x_opt, f_x_opt)
'Minimum: f([-5. 3.49999914]) = -136.250000'
from scipy.optimize import fsolve
f2 = lambda x : 3 * x**3 - x**2 + 1
fsolve(f2, 1.)
array([-0.5981935])
Vektorwertige Funktion $f_3(x): \mathbb{R}^2 \rightarrow \mathbb{R}^2$
from scipy.optimize import root
from numpy import sqrt
f3 = lambda x : [(x[0] + x[1])**2 - 4, x[0] + sqrt(x[1])]
f3_root = root(f3, [2., 1.])
f3_root
fjac: array([[-0.97680164, -0.21414609], [ 0.21414609, -0.97680164]]) fun: array([ 2.07300843e-12, 1.27009514e-13]) message: 'The solution converged.' nfev: 12 qtf: array([ -8.12325364e-09, 1.26538374e-09]) r: array([-4.8930465 , -3.96638873, 0.61326347]) status: 1 success: True x: array([-2., 4.])
f3_root.x
array([-2., 4.])