Hoofdstuk 10 — Prestaties en geheugenbeheer met NumPy

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

Doel van dit hoofdstuk:

• Begrijpen waar NumPy-snelheid vandaan komt en hoe je die benut.
• Leren hoe geheugenlayout, dtypes, views/kopieën en broadcasting de performance en het geheugengebruik beïnvloeden.
• Praktische strategieën voor minder allocaties, betere cache-lokaliteit en controle over multithreading.
• Werken met grote data (memory mapping, chunking) en tips om out-of-memory te vermijden.

Aanbevolen voorkennis: Hoofdstuk 6–9 (arrays, indexeren/broadcasting, ufuncs/aggregaties, lineaire algebra en numerieke stabiliteit).

10.1 Waarom is NumPy snel?

• Vectorisatie: operaties zijn geïmplementeerd in C (ufuncs/BLAS), vermijden Python-loop overhead.
• Contigu geheugen en eenvoudige dtypes bevorderen cache-hits en SIMD-vectorisatie.
• Aanroepen naar geoptimaliseerde bibliotheken (BLAS/LAPACK) voor lineaire algebra.

Vergelijking loop vs vectorisatie:

import numpy as np
from time import perf_counter

n = 2_000_000
x = np.random.default_rng(0).random(n)
y = np.random.default_rng(1).random(n)

# Python loop (for illustration; avoid in production)
start = perf_counter()
res_py = [xi + yi for xi, yi in zip(x, y)]
t_loop = perf_counter() - start

# Vectorized
start = perf_counter()
res_np = x + y
t_vec = perf_counter() - start

len(res_py), np.allclose(res_np, np.array(res_py)), t_loop, t_vec

10.2 Dtypes, alignement en geheugenvoetafdruk

• Kies het kleinste dtype dat nauwkeurig genoeg is (bijv. float32 i.p.v. float64) om geheugen en bandbreedte te besparen.
• itemsize en size bepalen de bytes: bytes = size * itemsize.
• Integers kunnen over-/onderlopen; floats hebben precisie/NaN/Inf.

import numpy as np

n = 10_000_000
x32 = np.empty(n, dtype=np.float32)
x64 = np.empty(n, dtype=np.float64)
mb32 = x32.nbytes / (1024**2)
mb64 = x64.nbytes / (1024**2)
mb32, mb64

10.3 Views vs. kopieën: kostprijs en detectie

• Views zijn gratis in geheugen (gedeelde data) en snel te maken; kopieën kosten O(n) tijd en geheugen.
• Fancy indexing (integer/boolean) maakt altijd een kopie; slicing meestal een view.
• Detecteer met np.shares_memory(a, b) en inspecteer a.flags voor contiguïteit.

import numpy as np

A = np.arange(12).reshape(3, 4)
B = A[:, 1:3]          # view
C = A[:, [1, 2]]       # copy (fancy indexing)
np.shares_memory(A, B), np.shares_memory(A, C)

10.4 Contiguïteit, strides en order (C vs Fortran)

• C-order: laatste index varieert het snelst (standaard). Fortran-order: eerste index varieert het snelst.
• Contigue arrays doorlopen in geheugenvolgorde bevordert cache-lokaliteit.
• Converteer layout indien nodig voor bibliotheken: np.ascontiguousarray, np.asfortranarray.

import numpy as np

A = np.arange(12).reshape(3, 4)          # C-order by default
AF = np.asfortranarray(A)
A.flags['C_CONTIGUOUS'], AF.flags['F_CONTIGUOUS'], A.strides, AF.strides

10.5 Broadcasting: krachtig maar let op explosies

• Broadcasting maakt tijdelijk grotere logische shapes; het resultaat wordt niet gematerialiseerd, maar sommige combinaties creëren enorme tijdelijke arrays bij vervolgbewerkingen.
• Tip: reken in de juiste volgorde en reduceer vroeg (bijv. som over de gebroadcastte as binnen dezelfde uitdrukking).

import numpy as np

X = np.random.default_rng(0).random((10_000, 128))   # (N, d)
w = np.random.default_rng(1).random(128)             # (d,)
# Goed: reduceer binnen de uitdrukking
s = np.sum(X * w, axis=1)    # (N,)

Valkuil: onbedoelde broadcast naar (N, N) voor paargewijze afstanden; gebruik efficiëntere vormen (@, einsum) of chunking.

10.6 Minder allocaties: in-place en out=

• Gebruik in-place operators (+=, -=, *=, np.add(..., out=...)) om tijdelijke arrays te vermijden.
• Let op dtype-compatibiliteit; anders forceert NumPy een kopie.

import numpy as np

x = np.arange(5, dtype=np.float64)
np.add(x, 1.0, out=x)    # in-place add
# Voorwaardelijk schrijven zonder extra allocaties
y = np.empty_like(x)
np.sqrt(x, out=y)
np.sqrt(x, out=y, where=(x > 0))
y

10.7 Operationele fusie en ufunc-idiomen

• Combineer ufuncs zodat NumPy intern minder tussenresultaten aanmaakt (niet alle ufuncs fuseren automatisch, maar out= en where= helpen).
• Gebruik ufunc-methoden (reduce, accumulate, outer) en aggregatieparameters (axis, keepdims) om werk te beperken.

import numpy as np

A = np.random.default_rng(0).standard_normal((1000, 1000))
# Minder allocaties: centreer en schaal in twee passes maar zonder extra kopieën
col_mean = A.mean(axis=0)
A -= col_mean
col_std = A.std(axis=0, ddof=0)
A /= col_std

10.8 BLAS/LAPACK en threading

• Matrixvermenigvuldiging en decomposities gebruiken BLAS/LAPACK (OpenBLAS, MKL, Accelerate).
• Multithreading wordt door de backend geregeld; te veel threads kan tot oversubscripie leiden bij nested parallelisme.
• Omgevingsvariabelen (bij processtart instellen): OPENBLAS_NUM_THREADS, MKL_NUM_THREADS, OMP_NUM_THREADS.

import numpy as np
np.show_config()   # laat zien welke BLAS/LAPACK eronder zit

Tip: in multiprocessing of binnen joblib/tqdm: beperk BLAS-threads per proces om thrash te vermijden.

10.9 NumExpr/Numba/Cython (wanneer verder gaan)

• np.vectorize is geen performance-tool (het is syntactic sugar over Python-loops). Vermijd voor snelheid.
• numexpr kan grote expresies sneller uitvoeren door chunking en multi-threading op cache-vriendelijke blokken.
• Numba JIT-compileert Python naar machinecode; nuttig voor custom kernels die niet goed uit te drukken zijn met ufuncs/broadcasting.
• Cython voor nog fijnmaziger controle; meer setup.

Doel van dit boek: eerst optimaal NumPy gebruiken; pas daarna JIT/extern als nodig.

10.10 Geheugenmapping voor grote datasets (np.memmap)

• Memory-mapped bestanden laten je arrays gebruiken die groter dan RAM kunnen zijn door paginagebaseerd IO.
• Werkt goed voor sequentiële/geslicede toegang; random access blijft schijfgebonden.

import numpy as np

# Create a memmap file
filename = 'data.bin'
N = 10_000_000
mm = np.memmap(filename, dtype='float32', mode='w+', shape=(N,))
mm[:100] = np.arange(100, dtype='float32')
mm.flush()
# Reopen read-only
mm_ro = np.memmap(filename, dtype='float32', mode='r', shape=(N,))
mm_ro[10:15]

Tips:

• Gebruik eenvoudige dtypes en contigue layouts.
• Sluit/flush (of laat GC doen) voor consistentie.

10.11 Chunking en streaming

• Verwerk grote arrays in stukken om piekgeheugen te beperken.
• Kies chunkgrootte die in cache/L3 past of afgestemd is op IO-doorvoer.

import numpy as np

rng = np.random.default_rng(0)
X = rng.random((1_000_000, 64), dtype=np.float64)
w = rng.random(64)

def chunked_dot(X, w, chunk=50_000):
    out = np.empty(X.shape[0], dtype=X.dtype)
    i = 0
    while i < X.shape[0]:
        j = min(i + chunk, X.shape[0])
        out[i:j] = X[i:j] @ w
        i = j
    return out

res = chunked_dot(X, w)
res.shape

10.12 Cache-lokaliteit en tiling

• Door iteratie in geheugenvolgorde te houden en blokken (tiles) te gebruiken, verbeter je cache-hits.
• Voor pure NumPy is dit vooral relevant bij eigen implementaties of bij het kiezen van C- vs Fortran-order.

import numpy as np

A = np.asfortranarray(np.arange(12).reshape(3, 4))
# Kolomoperaties kunnen dan cache-vriendelijker zijn
col_means = A.mean(axis=0)
A.flags['F_CONTIGUOUS'], col_means.shape

10.13 Profileren en meten

• Micro-benchmarking: IPython/Jupyter %timeit; of time.perf_counter.
• Inspecteer NumPy/BLAS-config: np.__version__, np.show_config().
• Memory-profielen: tracemalloc, memory_profiler (extern), array.nbytes en gc-inzichten.

import numpy as np
import sys

A = np.empty((2000, 2000), dtype=np.float64)
A.nbytes / (1024**2), sys.getsizeof(A)

10.14 Veelvoorkomende patronen en checklist

• Vermijd Python-elementloops; gebruik ufuncs/aggregaties en @ voor algebra.
• Controleer shapes vroeg; voorkom onbedoelde (N, N) constructies.
• Minimaliseer allocaties met in-place en out=; hergebruik buffers.
• Kies passende dtypes; overweeg float32 bij bandbreedte-gebonden workloads.
• Zorg voor contigu arrays als bibliotheken dat vereisen; converteer expliciet.
• Beperk BLAS-threads om oversubscriptie te voorkomen in parallelle workloads.
• Voor grote data: memmap, chunking en streaming.

10.15 Korte oefeningen

1. Meet het verschil in tijd tussen een Python-loop en een vectorized optelling op 5 miljoen elementen. Gebruik perf_counter en rapporteer de ratio.
2. Maak een grote array (5_000_000,) in float32 en float64. Bereken de geheugenvoetafdruk in MB en bespreek een scenario waarin float32 volstaat.
3. Toon met een klein voorbeeld dat fancy indexing een kopie maakt en slicing een view. Gebruik np.shares_memory en verander een element.
4. Implementeer een chunked dot-product y = X @ w voor X van shape (1_000_000, 128) en vergelijk de piekgeheugenbehoefte met de naïeve broadcast np.sum(X*w, axis=1).
5. Gebruik np.memmap om een bestand van 1e7 float32’s aan te maken, vul de eerste 100 waarden, sluit en heropen, en verifieer de inhoud.
6. Profileringsopdracht: gebruik %timeit of perf_counter om de impact van OPENBLAS_NUM_THREADS=1 vs standaard op A @ B (vorm (2000,2000)) te vergelijken. (Beschrijf hoe je de omgeving zou instellen.)

Werkboek en uitwerkingen:

• Oefeningen bij dit hoofdstuk kunnen later toegevoegd worden in ../workbook/hoofdstuk-10/oefeningen.md.
• Uitwerkingen kunnen later toegevoegd worden in ../solutions/hoofdstuk-10/uitwerkingen.md.

Verder lezen en vervolg:

• Hoofdstuk 6: basis over arrays, dtypes en geheugen.
• Hoofdstuk 7: broadcasting en vormmanipulatie (prestaties/vallenkuilen).
• Hoofdstuk 8: ufuncs, out=, where en aggregaties (allocaties beperken).
• Hoofdstuk 9: BLAS/SVD en numerieke stabiliteit.
• Hoofdstuk 11: interoperabiliteit en memory-views voor zero-copy patronen.