NumPy - Convolution (np.convolve)

Convolution mathematically combines two sequences by sliding one over the other, multiplying overlapping elements, and summing the results. For discrete sequences, the convolution of arrays `a` and...

Key Insights

  • Convolution is a fundamental mathematical operation that combines two sequences to produce a third sequence, widely used in signal processing, time series analysis, and image filtering
  • NumPy’s np.convolve() implements 1D discrete convolution with three modes (full, same, valid) that control output size, making it essential to choose the right mode for your use case
  • Understanding convolution’s computational complexity O(n*m) and when to switch to FFT-based methods for large arrays can dramatically improve performance in production systems

What is Convolution

Convolution mathematically combines two sequences by sliding one over the other, multiplying overlapping elements, and summing the results. For discrete sequences, the convolution of arrays a and v at position n is:

(a * v)[n] = Σ a[m] * v[n - m]

In NumPy, np.convolve(a, v, mode) computes this operation efficiently. The function slides the second array v (the kernel) across the first array a (the signal).

import numpy as np

# Simple convolution example
signal = np.array([1, 2, 3, 4, 5])
kernel = np.array([0.5, 1, 0.5])

result = np.convolve(signal, kernel, mode='valid')
print(result)  # [2.5 4.  5.5 7. ]

The kernel slides across the signal, computing weighted sums at each position.

Understanding Convolution Modes

The mode parameter controls how edge cases are handled and determines the output size. This is critical for practical applications.

Full Mode

Returns the full convolution at each point of overlap. Output length is len(a) + len(v) - 1.

signal = np.array([1, 2, 3, 4, 5])
kernel = np.array([1, 1, 1])

full = np.convolve(signal, kernel, mode='full')
print(f"Input length: {len(signal)}")
print(f"Kernel length: {len(kernel)}")
print(f"Output length: {len(full)}")
print(f"Output: {full}")

# Output:
# Input length: 5
# Kernel length: 3
# Output length: 7
# Output: [1 3 6 9 12 9 5]

Full mode includes partial overlaps at the edges, useful when you need all convolution products.

Same Mode

Returns output of the same length as the first input array. Centers the kernel over each element.

same = np.convolve(signal, kernel, mode='same')
print(f"Output length: {len(same)}")
print(f"Output: {same}")

# Output:
# Output length: 5
# Output: [3 6 9 12 9]

This mode is ideal for filtering operations where you want to preserve the original signal length.

Valid Mode

Returns only the convolution computed without zero-padding. Output length is max(len(a), len(v)) - min(len(a), len(v)) + 1.

valid = np.convolve(signal, kernel, mode='valid')
print(f"Output length: {len(valid)}")
print(f"Output: {valid}")

# Output:
# Output length: 3
# Output: [6 9 12]

Valid mode ensures all output values use complete kernel coverage, eliminating edge effects.

Moving Average Filter

A practical application of convolution is computing moving averages for smoothing time series data.

import numpy as np
import matplotlib.pyplot as plt

# Generate noisy data
np.random.seed(42)
time = np.linspace(0, 10, 200)
signal = np.sin(time) + 0.3 * np.random.randn(200)

# Create moving average kernel
window_size = 10
kernel = np.ones(window_size) / window_size

# Apply convolution
smoothed = np.convolve(signal, kernel, mode='same')

# Edge correction for 'same' mode
# Edges use fewer points, so we need to adjust
edge_correction = np.ones(len(signal))
edge_correction[:window_size//2] = np.arange(window_size//2, window_size) / window_size
edge_correction[-window_size//2:] = np.arange(window_size, window_size//2, -1) / window_size

smoothed_corrected = np.convolve(signal, kernel, mode='same') / edge_correction

print(f"Original variance: {np.var(signal):.4f}")
print(f"Smoothed variance: {np.var(smoothed_corrected):.4f}")

This moving average implementation handles edge cases properly by accounting for the reduced number of samples at boundaries.

Derivative Approximation

Convolution can approximate derivatives using finite difference kernels.

# First derivative kernel (central difference)
derivative_kernel = np.array([-1, 0, 1]) / 2

# Generate a smooth function
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)

# Compute derivative via convolution
dy_dx = np.convolve(y, derivative_kernel, mode='same')

# Analytical derivative for comparison
analytical_derivative = np.cos(x)

# Compare at center points (edges have boundary effects)
center_slice = slice(10, -10)
error = np.mean(np.abs(dy_dx[center_slice] - analytical_derivative[center_slice]))
print(f"Mean absolute error: {error:.6f}")

The central difference kernel [-1, 0, 1]/2 provides a second-order accurate approximation of the first derivative.

Pattern Detection

Convolution excels at detecting patterns in sequences by using the pattern as a kernel.

# Signal with embedded pattern
signal = np.array([0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0])
pattern = np.array([1, 1, 0])

# Convolve to find pattern matches
correlation = np.convolve(signal, pattern[::-1], mode='valid')

# High values indicate pattern matches
threshold = 2
matches = np.where(correlation >= threshold)[0]

print(f"Pattern found at indices: {matches}")
print(f"Correlation values: {correlation}")

# Output:
# Pattern found at indices: [3 8]
# Correlation values: [1 1 1 2 1 1 1 1 2 1 1 1 1]

Note the kernel reversal (pattern[::-1]) converts convolution to correlation, which is what we want for pattern matching.

Performance Considerations

Convolution has O(n*m) complexity. For large arrays, this becomes expensive.

import time

# Compare performance for different sizes
sizes = [100, 1000, 10000]
kernel = np.ones(50)

for size in sizes:
    signal = np.random.randn(size)
    
    start = time.time()
    result = np.convolve(signal, kernel, mode='same')
    elapsed = time.time() - start
    
    print(f"Size {size:5d}: {elapsed*1000:.2f} ms")

For large convolutions, consider FFT-based methods:

from scipy import signal as sp_signal

# Large arrays
large_signal = np.random.randn(100000)
large_kernel = np.random.randn(1000)

# Direct convolution
start = time.time()
result_direct = np.convolve(large_signal, large_kernel, mode='same')
time_direct = time.time() - start

# FFT-based convolution
start = time.time()
result_fft = sp_signal.fftconvolve(large_signal, large_kernel, mode='same')
time_fft = time.time() - start

print(f"Direct convolution: {time_direct:.3f} s")
print(f"FFT convolution: {time_fft:.3f} s")
print(f"Speedup: {time_direct/time_fft:.1f}x")
print(f"Max difference: {np.max(np.abs(result_direct - result_fft)):.2e}")

FFT-based convolution becomes faster when min(n, m) > ~500, with speedups of 10-100x for large arrays.

2D Convolution Alternative

While np.convolve() handles only 1D arrays, you can extend the concept to 2D using nested convolutions or use scipy.signal.convolve2d for image processing:

from scipy.signal import convolve2d

# Simple edge detection on a synthetic image
image = np.zeros((10, 10))
image[3:7, 3:7] = 1  # Square in center

# Sobel operator for edge detection
sobel_x = np.array([[-1, 0, 1],
                     [-2, 0, 2],
                     [-1, 0, 1]])

edges = convolve2d(image, sobel_x, mode='same', boundary='fill')
print("Vertical edges detected:")
print(edges.astype(int))

For production image processing, libraries like scipy.ndimage or opencv provide optimized implementations with additional features like boundary handling and multi-channel support.

Common Pitfalls

Kernel Ordering: Convolution mathematically reverses the kernel. For symmetric kernels this doesn’t matter, but for asymmetric ones:

signal = np.array([1, 2, 3, 4])
kernel = np.array([1, 10])

conv_result = np.convolve(signal, kernel, mode='valid')
print(f"Convolution: {conv_result}")  # [12 23 34]

# Manual calculation shows kernel reversal
# Position 0: 1*10 + 2*1 = 12
# Position 1: 2*10 + 3*1 = 23

Memory Usage: Full mode creates arrays of size n+m-1, which can exhaust memory for large inputs. Always consider output size before calling np.convolve().

Type Promotion: NumPy promotes integer types to prevent overflow, but this can surprise you:

a = np.array([1, 2, 3], dtype=np.int8)
b = np.array([1, 1], dtype=np.int8)
result = np.convolve(a, b)
print(result.dtype)  # int64 on most systems

Convolution remains a cornerstone operation in numerical computing. Master these modes and performance characteristics to build robust signal processing and data analysis pipelines.

Liked this? There's more.

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