Hoofdstuk 9 — Lineaire algebra, random number generation en numerieke stabiliteit

Let op: De tekst is in het Nederlands; alle Python codevoorbeelden zijn in het Engels.

Doel van dit hoofdstuk:

Krachtig werken met lineaire algebra in NumPy: vectoren, matrices, products (@, matmul, einsum).
Lineaire systemen oplossen, least squares, rangen, condities en de pseudo-inverse.
Overzicht van decomposities (SVD, eigenwaarden/-vectoren) en hun toepassingen.
Moderne random number generation (RNG) met numpy.random.Generator voor reproduceerbare experimenten.
Numerieke stabiliteit: conditionering, schaalproblemen, stabiele patronen en valkuilen.

Aanbevolen voorkennis: Hoofdstuk 6–8 (arrays, indexeren/broadcasting, ufuncs/aggregaties), basis lineaire algebra.

9.1 Vectors en matrices in NumPy

NumPy vertegenwoordigt vectoren en matrices als ndarray’s. Een 1D-array is een vector; 2D is een matrix.

import numpy as np

x = np.array([1.0, 2.0, 3.0])          # vector (3,)
A = np.array([[1.0, 2.0],
              [3.0, 4.0],
              [5.0, 6.0]])             # matrix (3, 2)

print(x.shape)  # (3,)
print(A.shape)  # (3, 2)

Kolom- en rijvectoren expliciteren met None/np.newaxis (zie H6–H7):

col = x[:, None]  # (3, 1)
row = x[None, :]  # (1, 3)

9.2 Matrix- en vectorproducten: `@`, `matmul`, `einsum`

@ en np.matmul voeren matrixvermenigvuldiging uit (geen elementgewijs product).
np.dot is voor 1D gelijk aan inwendig product; voor 2D komt het overeen met matmul.
np.einsum geeft expliciete controle over indexcontracties (leesbaar bij complexe patronen).

u = np.array([1.0, 2.0, 3.0])      # (3,)
B = np.array([[1.0, 0.0, 2.0],
              [0.0, 1.0, 3.0]])    # (2, 3)

# matrix-vector
y = B @ u            # shape (2,)

# matrix-matrix
C = np.array([[1.0, 2.0], [3.0, 4.0]])   # (2,2)
D = B @ C             # (2,3) @ (2,2) -> ValueError (incompatible)
# Correct: (3,2) @ (2,2)
E = B.T @ C           # shape (3, 2) @ (2, 2) -> (3, 2)

# einsum voorbeelden
# inwendig product van u met zichzelf
ip = np.einsum('i,i->', u, u)      # scalar
# batch matrix-vector: (n,m) x (m,) voor meerdere matrices/vektoren
M = np.arange(12.0).reshape(3, 4)
w = np.arange(4.0)
res = M @ w                         # (3,)

Tip: controleer shapes expliciet om dimensionele fouten te voorkomen.

9.3 Normen, lengtes en projecties

Vectornormen: np.linalg.norm(x, ord=2) (Euclidisch); ord=1, ord=np.inf voor L1/L∞.
Matrixnormen: ook via np.linalg.norm; let op definitie (induced norms).
Normaliseren (unit vector) en projecties.

x = np.array([3.0, 4.0])
length = np.linalg.norm(x)          # 5.0
unit = x / length                   # [0.6, 0.8]

# Projecteer vector v op richting u (niet-nul)
u = np.array([1.0, 2.0, 2.0])
v = np.array([2.0, 1.0, 0.0])
proj_v_on_u = (v @ u) / (u @ u) * u

Orthogonaliteit: u @ v == 0 (binnen numerieke tolerantie) betekent loodrecht.

9.4 Lineaire systemen: `solve`, `lstsq`, `pinv`

Los systemen Ax = b op met stabiele routines; vermijd handmatig inverteren.

A = np.array([[3.0, 2.0],
              [1.0, 2.0]])
b = np.array([5.0, 5.0])

# Exact vierkant en niet-singulier
x = np.linalg.solve(A, b)

# Over- of onderbepaald: least squares
M = np.array([[1.0, 1.0],
              [1.0, 2.0],
              [1.0, 3.0]])
y = np.array([1.0, 2.0, 2.5])
coef, resid, rank, s = np.linalg.lstsq(M, y, rcond=None)

# Pseudo-inverse (Moore-Penrose)
M_pinv = np.linalg.pinv(M)     # gebaseerd op SVD
coef2 = M_pinv @ y

Belangrijk:

Gebruik solve/lstsq in plaats van inv(A) @ b. Inversen zijn duurder en numeriek minder stabiel.
Controleer rank en s (singuliere waarden) om conditionering te beoordelen.

9.5 Rang, determinant, condities, condition number

Rang (np.linalg.matrix_rank(A)) geeft het aantal lineair onafhankelijke rijen/kolommen.
Determinant (np.linalg.det(A)) vertelt over singulariteit maar is gevoelig voor schaal.
Condition number (np.linalg.cond(A)) meet gevoeligheid; hogere waarden impliceren instabiliteit.

A = np.array([[1.0, 2.0], [2.0, 4.000001]])
print(np.linalg.matrix_rank(A))
print(np.linalg.det(A))
print(np.linalg.cond(A))

Regels:

Slecht geconditioneerde systemen vergroten afrondingsfouten; her-schaal kolommen/rijen of gebruik regularisatie.
Inspecteer cond(A) en pas model/algoritme aan bij grote waarden.

9.6 Eigenwaarden/-vectoren en SVD (overzicht)

Eigen-decompositie: A v = λ v voor vierkante matrices.
SVD: A = U Σ Vᵀ werkt voor rechthoekige matrices en is robuust; nuttig voor PCA, pseudo-inversen en rangen.

A = np.array([[2.0, 0.0], [0.0, 1.0]])
# Eigenwaarden/-vectoren
w, V = np.linalg.eig(A)
# SVD
U, S, VT = np.linalg.svd(A, full_matrices=False)

# Reconstructie
A_rec = (U * S) @ VT

Opmerkingen:

Voor symmetrische (Hermitische) matrices is np.linalg.eigh numeriek stabieler dan eig.
SVD levert singuliere waarden in aflopende volgorde; kleine waarden wijzen op bijna-singulariteit.

9.7 Geen inverses bouwen: `solve`/`cho_solve`

Als A symmetrisch positief-definiet (SPD) is, gebruik Cholesky-factorisatie.

# SPD voorbeeld
S = np.array([[4.0, 2.0], [2.0, 3.0]])
L = np.linalg.cholesky(S)   # S = L @ L.T
# los Sx=b op via forward/back substitution
b = np.array([1.0, 2.0])
# Gebruik solve met triangular matrices via scipy.linalg (indien beschikbaar) of terugvallen op np.linalg.solve
x = np.linalg.solve(S, b)

Richtlijnen:

Vermijd inv(A). Gebruik specifieke oplossers.
Hergebruik factorisaties waar mogelijk.

9.8 Moderne RNG: `numpy.random.Generator`

Gebruik de nieuwe API voor reproduceerbare en onafhankelijke streams.

rng = np.random.default_rng(seed=42)   # PCG64 als default

# Distributies
u = rng.random(5)                      # Uniform [0,1)
z = rng.standard_normal(5)             # N(0,1)
k = rng.integers(low=0, high=10, size=5)

# Multivariate normaal
mean = np.array([0.0, 0.0])
cov = np.array([[1.0, 0.8], [0.8, 2.0]])  # symmetrisch PSD
samples = rng.multivariate_normal(mean, cov, size=1000, check_valid='warn')

# Permutaties
arr = np.arange(10)
rng.shuffle(arr)        # in-place
p = rng.permutation(5)  # nieuw array

Best practices:

Gebruik één Generator per task/worker; vermijd globale np.random state.
Zaai (seed) expliciet voor reproduceerbaarheid; gebruik SeedSequence voor hiërarchische seeds in parallelle settings.

9.9 Numerieke stabiliteit in de praktijk

Schaal feature-matrices (kolomgewijs centreren/standaardiseren) om conditionering te verbeteren.
Gebruik stabiele vormen: logsumexp in plaats van naïeve np.exp + sum + log.
Kies geschikte dtypes (float64 standaard; float32 voor snelheid/geheugen, maar controleer fouten).
Vermijd subtractie van bijna gelijke getallen (catastrofale annulatie).

# Softmax (stabiel per rij)
X = np.array([[1.0, 2.0, 3.0], [1.0, 2.0, 4.0]])
Xmax = X.max(axis=1, keepdims=True)
exp = np.exp(X - Xmax)
softmax = exp / exp.sum(axis=1, keepdims=True)

# Log-sum-exp truc
logsumexp = Xmax[:, 0] + np.log(exp.sum(axis=1))

Conditionering illustratie:

A = np.array([[1.0, 1.0], [1.0, 1.000001]])
b = np.array([2.0, 2.000001])
sol = np.linalg.solve(A, b)
print('cond(A)=', np.linalg.cond(A), ' solution=', sol)

9.10 Prestaties en geheugen

NumPy gebruikt BLAS/LAPACK (via OpenBLAS/MKL) voor @, svd, eig, etc.
Zorg voor contigue arrays met juiste dtype (float64 of float32) voor maximale performance.
Voorkom enorme tussenresultaten: gebruik @/einsum i.p.v. broadcast+sum wanneer het eigenlijk een matrixproduct is.
Batched operaties: @ en matmul ondersteunen broadcasting over leading dims (3D+).

# Batched matmul: (batch,n,m) @ (batch,m,p) -> (batch,n,p)
A = np.random.default_rng(0).random((10, 3, 4))
B = np.random.default_rng(1).random((10, 4, 5))
C = A @ B
assert C.shape == (10, 3, 5)

9.11 Veelvoorkomende patronen en valkuilen

Z-standaardisatie: Z = (X - X.mean(axis=0)) / X.std(axis=0, ddof=0).
Regressie met lstsq i.p.v. inv(X.T@X) @ X.T @ y.
Projectiematrix met orthonormale kolommen Q: P = Q @ Q.T (symmetrisch idempotent).
Controleer symmetrie (np.allclose(A, A.T)) voordat je eigh gebruikt.

9.12 Korte oefeningen

Geef A (shape (3,3)) en b (shape (3,)). Los Ax=b op met np.linalg.solve en verifieer door terug te vermenigvuldigen.
Maak een matrix X (shape (n,d)) en bereken de covariantiematrix S = X_centered.T @ X_centered / (n-1). Gebruik np.linalg.eig of np.linalg.svd om de hoofdrichtingen te vinden.
Demonstreer het verschil tussen inv(A) @ b en solve(A, b) qua snelheid/nauwkeurigheid op een willekeurige SPD matrix.
Gebruik de moderne RNG om 1000 samples uit N(0, I) in 2D te trekken en schat de empirische covariantie; vergelijk met de theoretische.
Laat met een slecht geconditioneerde matrix zien dat kleine verstoringen in b leiden tot grote verschillen in x.

Werkboek en uitwerkingen:

Oefeningen bij dit hoofdstuk vind je in ../workbook/hoofdstuk-09/oefeningen.md.
Uitwerkingen staan in ../solutions/hoofdstuk-09/uitwerkingen.md.

Verder lezen en vervolg:

Hoofdstuk 8: numerieke stabiliteit van ufuncs en aggregaties.
Hoofdstuk 10: prestaties en geheugenbeheer met NumPy.
Appendix/verdere bronnen: LAPACK/BLAS, Golub & Van Loan (Matrix Computations).

Keyboard shortcuts

Python Data