Eric Campbell 4 min read Python

NumPy - Vectorization and Performance

• Vectorized NumPy operations execute 10-100x faster than Python loops by leveraging pre-compiled C code and SIMD instructions that process multiple data elements simultaneously

Key Insights

• Vectorized NumPy operations execute 10-100x faster than Python loops by leveraging pre-compiled C code and SIMD instructions that process multiple data elements simultaneously • Memory layout matters: contiguous arrays in C-order enable cache-friendly access patterns, while unnecessary copies and mixed data types destroy performance gains • Understanding when to use vectorization versus native Python is critical—small datasets and complex conditional logic often perform better without NumPy overhead

Why Vectorization Outperforms Loops

Vectorization replaces explicit Python loops with array operations that execute in compiled C code. Python’s interpreter overhead—type checking, reference counting, and dynamic dispatch—disappears when NumPy handles iteration internally.

import numpy as np
import time

# Scalar approach with Python loop
def sum_python(arr):
    total = 0
    for x in arr:
        total += x
    return total

# Vectorized approach
def sum_numpy(arr):
    return np.sum(arr)

data = list(range(1_000_000))
arr = np.array(data)

# Python loop
start = time.perf_counter()
result1 = sum_python(data)
python_time = time.perf_counter() - start

# NumPy vectorization
start = time.perf_counter()
result2 = sum_numpy(arr)
numpy_time = time.perf_counter() - start

print(f"Python: {python_time:.4f}s")
print(f"NumPy: {numpy_time:.4f}s")
print(f"Speedup: {python_time/numpy_time:.1f}x")

On typical hardware, this produces 50-100x speedups. The gap widens with mathematical operations beyond simple addition.

Broadcasting: Implicit Vectorization

Broadcasting extends vectorization to arrays with different shapes, eliminating the need for explicit loops or array replication.

# Without broadcasting (inefficient)
matrix = np.random.rand(1000, 1000)
row_means = np.mean(matrix, axis=1)

# Manual broadcasting - creates unnecessary copy
row_means_expanded = np.tile(row_means.reshape(-1, 1), (1, 1000))
centered_manual = matrix - row_means_expanded

# With broadcasting (efficient)
centered = matrix - row_means[:, np.newaxis]

# Verify memory efficiency
print(f"Manual copy size: {row_means_expanded.nbytes / 1e6:.2f} MB")
print(f"Broadcasting overhead: minimal (no copy)")

Broadcasting rules apply when:

  1. Arrays have the same number of dimensions, or
  2. One array can be prepended with dimensions of size 1
  3. Dimensions are compatible (equal or one is 1)
# Practical example: normalize columns
data = np.random.randn(10000, 50)

# Column-wise statistics
col_mean = data.mean(axis=0)  # Shape: (50,)
col_std = data.std(axis=0)    # Shape: (50,)

# Broadcasting automatically aligns shapes
normalized = (data - col_mean) / col_std  # (10000, 50) - (50,) / (50,)

# Verify normalization
print(f"New means: {normalized.mean(axis=0)[:5]}")  # ~0
print(f"New stds: {normalized.std(axis=0)[:5]}")    # ~1

Memory Layout and Contiguity

NumPy arrays store data in contiguous memory blocks. The layout—row-major (C-order) or column-major (Fortran-order)—dramatically affects cache performance.

# Create arrays with different memory layouts
c_order = np.random.rand(5000, 5000)  # Default C-contiguous
f_order = np.asfortranarray(c_order)  # Fortran-contiguous

print(f"C-order flags: {c_order.flags['C_CONTIGUOUS']}")
print(f"F-order flags: {f_order.flags['F_CONTIGUOUS']}")

# Row-wise operations favor C-order
start = time.perf_counter()
row_sum_c = c_order.sum(axis=1)
c_time = time.perf_counter() - start

start = time.perf_counter()
row_sum_f = f_order.sum(axis=1)
f_time = time.perf_counter() - start

print(f"C-order row sum: {c_time:.4f}s")
print(f"F-order row sum: {f_time:.4f}s")
print(f"Ratio: {f_time/c_time:.2f}x slower")

Slicing can create non-contiguous views that degrade performance:

arr = np.random.rand(1000, 1000)

# Contiguous slice
contiguous = arr[:, :500]
print(f"Contiguous: {contiguous.flags['C_CONTIGUOUS']}")

# Non-contiguous slice (every other column)
strided = arr[:, ::2]
print(f"Strided: {strided.flags['C_CONTIGUOUS']}")

# Force contiguity when needed
contiguous_copy = np.ascontiguousarray(strided)

# Performance comparison
start = time.perf_counter()
_ = strided.sum()
strided_time = time.perf_counter() - start

start = time.perf_counter()
_ = contiguous_copy.sum()
contiguous_time = time.perf_counter() - start

print(f"Speedup from contiguity: {strided_time/contiguous_time:.2f}x")

Universal Functions (ufuncs)

Ufuncs are vectorized wrappers around scalar operations, providing element-wise operations with optimal performance.

# Custom ufunc example
def sigmoid_python(x):
    return 1 / (1 + np.exp(-x))

# NumPy's vectorize (not truly compiled, but cleaner)
sigmoid_vec = np.vectorize(sigmoid_python)

# Proper vectorized implementation
def sigmoid_numpy(x):
    return 1 / (1 + np.exp(-x))

x = np.linspace(-10, 10, 1_000_000)

# Compare performance
start = time.perf_counter()
result1 = sigmoid_numpy(x)
numpy_time = time.perf_counter() - start

start = time.perf_counter()
result2 = sigmoid_vec(x)
vectorize_time = time.perf_counter() - start

print(f"Vectorized: {numpy_time:.4f}s")
print(f"np.vectorize: {vectorize_time:.4f}s")
print(f"Note: np.vectorize is syntactic sugar, not performance optimization")

Chaining ufuncs efficiently:

# Inefficient: creates intermediate arrays
x = np.random.rand(1_000_000)
result = np.sqrt(np.square(x) + 1)

# Better: minimize intermediate allocations
result = np.sqrt(x * x + 1)

# Best: use out parameter to avoid allocations
out = np.empty_like(x)
np.square(x, out=out)
out += 1
np.sqrt(out, out=out)

When Vectorization Fails

Not all operations benefit from vectorization. Conditional logic and small datasets often perform worse.

# Complex conditional logic
def process_conditional(arr):
    """Poor vectorization candidate"""
    result = np.empty_like(arr)
    for i, val in enumerate(arr):
        if val < 0:
            result[i] = val ** 2
        elif val < 10:
            result[i] = val * 2
        else:
            result[i] = np.log(val)
    return result

# Vectorized attempt (not necessarily faster)
def process_vectorized(arr):
    result = np.empty_like(arr)
    mask1 = arr < 0
    mask2 = (arr >= 0) & (arr < 10)
    mask3 = arr >= 10
    
    result[mask1] = arr[mask1] ** 2
    result[mask2] = arr[mask2] * 2
    result[mask3] = np.log(arr[mask3])
    return result

# Small dataset benchmark
small = np.random.randn(100)
large = np.random.randn(1_000_000)

for data, label in [(small, "Small"), (large, "Large")]:
    start = time.perf_counter()
    _ = process_conditional(data)
    loop_time = time.perf_counter() - start
    
    start = time.perf_counter()
    _ = process_vectorized(data)
    vec_time = time.perf_counter() - start
    
    print(f"{label} - Loop: {loop_time:.6f}s, Vec: {vec_time:.6f}s")

Practical Optimization Patterns

# Pattern 1: Avoid repeated array creation
def inefficient_accumulation(n):
    result = np.array([])
    for i in range(n):
        result = np.append(result, i ** 2)
    return result

def efficient_accumulation(n):
    result = np.empty(n)
    for i in range(n):
        result[i] = i ** 2
    return result

def best_accumulation(n):
    indices = np.arange(n)
    return indices ** 2

# Pattern 2: Use views instead of copies
large_array = np.random.rand(10000, 10000)

# Copy (expensive)
subset_copy = large_array[1000:2000, 1000:2000].copy()

# View (cheap, but be careful with modifications)
subset_view = large_array[1000:2000, 1000:2000]

# Pattern 3: Leverage einsum for complex operations
a = np.random.rand(100, 50)
b = np.random.rand(50, 80)
c = np.random.rand(80, 60)

# Traditional approach
result1 = a @ b @ c

# Einsum (often faster, more explicit)
result2 = np.einsum('ij,jk,kl->il', a, b, c)

# Einsum for batch operations
batch = np.random.rand(32, 100, 50)
weights = np.random.rand(50, 10)

# Batch matrix multiplication
output = np.einsum('bij,jk->bik', batch, weights)
print(f"Output shape: {output.shape}")  # (32, 100, 10)

Profiling and Measurement

Always measure before optimizing. Use appropriate tools to identify bottlenecks.

import numpy as np
from time import perf_counter

def benchmark(func, *args, iterations=100):
    """Simple benchmark wrapper"""
    times = []
    for _ in range(iterations):
        start = perf_counter()
        func(*args)
        times.append(perf_counter() - start)
    return np.median(times)

# Example usage
data = np.random.rand(10000, 1000)

time1 = benchmark(np.sum, data, axis=0)
time2 = benchmark(np.mean, data, axis=0)

print(f"Sum: {time1*1000:.3f}ms")
print(f"Mean: {time2*1000:.3f}ms")

Vectorization transforms Python from a slow scripting language into a competitive numerical computing platform. The key is understanding when array operations execute in optimized C code versus when Python overhead dominates. Master memory layout, broadcasting, and ufuncs, but always profile real workloads before assuming vectorization improves performance.

Liked this? There's more.

Every week: one practical technique, explained simply, with code you can use immediately.