Eine Variable kann als Wert mehr als nur einen Skalar beinhalten.
Der Schritt zu Vektoren, Matrizen und Tensoren öffnet Tür und Tor um numerische Berechnungen ausdrücken zu können.
Die Programmbibliothek NumPy liefert mit numpy.ndarray
einen Rang-n Tensor,
der effizientes rechnen mit Vektoren und Matritzen erlaubt.
All diesen Berechnungen auf Basis dieses np.ndarray
ist gemeinsam,
dass sie dank der Vektorisierung viel 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 :-)
Die Bibliothek numpy
wird üblicherweise immer mit dem Kurznamen Alias np
importiert:
import numpy as np
Beispiel: Die Variablen $v$ und $w$ soll hier 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])
w = np.random.rand(3)
v
array([ 4, 3, -1])
w
array([ 0.54487339, 0.68192372, 0.11945565])
v + w
array([ 4.54487339, 3.68192372, -0.88054435])
v - w
array([ 3.45512661, 2.31807628, -1.11945565])
v * w
array([ 2.17949355, 2.04577117, -0.11945565])
Das (generalisierte) Skalarprodukt von $v$ und $w$
v.dot(w)
4.1058090694632829
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])
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]])
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 :
steht hierbei für die gesamte Spalte, n:m
für ein Intervall, wobei n
oder m
weggelassen werden kann. Negative Werte werden vom Ende weg gezählt.
m
array([[ 1, -2, 3, 8], [ 4, 6, 9, -9], [-1, 0, 5, 2], [11, 13, 15, 16]])
m[:, 0]
array([ 1, 4, -1, 11])
m[1, :]
array([ 4, 6, 9, -9])
m[2, 2]
5
m[:, -1]
array([ 8, -9, 2, 16])
m[2, [0, 2]]
array([-1, 5])
m[-1:]
array([[11, 13, 15, 16]])
m[:-1, :-2]
array([[ 1, -2], [ 4, 6], [-1, 0]])
Die .diagonal()
Funktion gibt die Diagonalelemente zurück.
m.diagonal()
array([ 1, 6, 5, 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, 5, 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, 5, 0, -1], [16, 15, 13, 11]])
m.T
array([[ 1, 4, -1, 11], [-2, 6, 0, 13], [ 3, 9, 5, 15], [ 8, -9, 2, 16]])
m[:, 0].dot(m[:, 1])
165
m.diagonal()
array([ 1, 6, 5, 16])
Darüber hinaus gibt es einige hilfreiche Funktionen, die das erzeugen von Punkten und eine Auswertung über einen Raster 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 equally spaced points from 0 to 1
np.linspace(0, 1, 10)
array([ 0. , 0.11111111, 0.22222222, 0.33333333, 0.44444444, 0.55555556, 0.66666667, 0.77777778, 0.88888889, 1. ])
10 equally spaced points on a logarithmic scale from base^1 to base^2 with base = 10
np.logspace(1, 2, 10, base=10)
array([ 10. , 12.91549665, 16.68100537, 21.5443469 , 27.82559402, 35.93813664, 46.41588834, 59.94842503, 77.42636827, 100. ])
Zur Auswertung von Funktionen über mehrere Achsen, welche aus zwei oder mehr Vektoren konstruiert werden.
xx, yy = np.meshgrid([1, 2, 3], [0, 10, 20])
xx
array([[1, 2, 3], [1, 2, 3], [1, 2, 3]])
yy
array([[ 0, 0, 0], [10, 10, 10], [20, 20, 20]])
Dies ist zur Auswertung von mehrstellige Funktionen sehr praktisch:
1 + xx + 2 * yy
array([[ 2, 3, 4], [22, 23, 24], [42, 43, 44]])
Indexmatrix
xx, yy = np.indices((3, 2))
xx
array([[0, 0], [1, 1], [2, 2]])
yy
array([[0, 1], [0, 1], [0, 1]])
Die Dimensionierung lässt sich jederzeit Ändern, es muss nur die Anzahl der Elemente gleich bleiben.
c = np.arange(25)
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.
c = np.arange(8)
c = c.reshape((-1, 4))
c
array([[0, 1, 2, 3], [4, 5, 6, 7]])
c.shape
(2, 4)
c.shape = (2, 2, 2)
c
array([[[0, 1], [2, 3]], [[4, 5], [6, 7]]])
c.ndim
3
c = np.arange(8).reshape((2, -1))
c = c[:, np.newaxis]
c
array([[[0, 1, 2, 3]], [[4, 5, 6, 7]]])
c.shape
(2, 1, 4)
c.ndim
3
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, die einem der Argumente fehlt, die entsprechenden Werte mehrfach wiederholt werden, sodass die Dimensionierung in dem 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]])
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)
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]])
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]]])
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(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
status: 0 success: True nfev: 64 fun: -0.99999999999996247 x: array([ 0.25 , 0.99999982, 1. ]) message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' jac: array([ 1.53210777e-06, -4.10782519e-07, -9.99999905e-01]) nit: 13
x_opti = opti1["x"]
x_opti
array([ 0.25 , 0.99999982, 1. ])
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
status: 1 success: True qtf: array([ -8.12325364e-09, 1.26538374e-09]) nfev: 12 r: array([-4.8930465 , -3.96638873, 0.61326347]) fun: array([ 2.07300843e-12, 1.27009514e-13]) x: array([-2., 4.]) message: 'The solution converged.' fjac: array([[-0.97680164, -0.21414609], [ 0.21414609, -0.97680164]])
f3_root.x
array([-2., 4.])