NumPy (Numerical Python) is the foundational library for numerical computing in Python. If you’ve worked with data science, machine learning, image processing, or scientific computing in Python, you’ve almost certainly used NumPy — whether directly or through libraries built on top of it like pandas, scikit-learn, TensorFlow, and OpenCV.
Here’s why NumPy matters:
In this tutorial, we’ll go deep on NumPy arrays — from creation to manipulation, from indexing to linear algebra. By the end, you’ll have a solid, practical understanding of the library that underpins nearly all of Python’s data stack.
NumPy is available via pip. If you don’t have it installed yet:
pip install numpy
If you’re using Anaconda, NumPy comes pre-installed. You can verify your installation:
import numpy as np print(np.__version__)
The convention of importing NumPy as np is universal in the Python ecosystem. Stick with it — every tutorial, Stack Overflow answer, and library documentation assumes this alias.
NumPy arrays (ndarray objects) are the core data structure. There are several ways to create them, each suited to different situations.
The most straightforward way to create a NumPy array is from an existing Python list or tuple:
import numpy as np
# 1D array
a = np.array([1, 2, 3, 4, 5])
print(a)
# Output: [1 2 3 4 5]
# 2D array (matrix)
b = np.array([[1, 2, 3],
[4, 5, 6]])
print(b)
# Output:
# [[1 2 3]
# [4 5 6]]
# 3D array
c = np.array([[[1, 2], [3, 4]],
[[5, 6], [7, 8]]])
print(c.shape)
# Output: (2, 2, 2)
# Specifying data type explicitly
d = np.array([1, 2, 3], dtype=np.float64)
print(d)
# Output: [1. 2. 3.]
When you need arrays pre-filled with zeros or ones (common for initializing weight matrices, accumulators, or masks):
# 1D array of zeros zeros_1d = np.zeros(5) print(zeros_1d) # Output: [0. 0. 0. 0. 0.] # 2D array of zeros (3 rows, 4 columns) zeros_2d = np.zeros((3, 4)) print(zeros_2d) # Output: # [[0. 0. 0. 0.] # [0. 0. 0. 0.] # [0. 0. 0. 0.]] # 1D array of ones ones_1d = np.ones(4) print(ones_1d) # Output: [1. 1. 1. 1.] # 2D array of ones with integer type ones_int = np.ones((2, 3), dtype=np.int32) print(ones_int) # Output: # [[1 1 1] # [1 1 1]] # Full array with a custom fill value filled = np.full((2, 3), 7) print(filled) # Output: # [[7 7 7] # [7 7 7]] # Identity matrix eye = np.eye(3) print(eye) # Output: # [[1. 0. 0.] # [0. 1. 0.] # [0. 0. 1.]]
np.arange() works like Python’s range() but returns an array. np.linspace() creates evenly spaced values between two endpoints — extremely useful for plotting and numerical methods.
# arange: start, stop (exclusive), step a = np.arange(0, 10, 2) print(a) # Output: [0 2 4 6 8] # arange with float step b = np.arange(0, 1, 0.2) print(b) # Output: [0. 0.2 0.4 0.6 0.8] # linspace: start, stop (inclusive), number of points c = np.linspace(0, 1, 5) print(c) # Output: [0. 0.25 0.5 0.75 1. ] # linspace is ideal for generating x-values for plots x = np.linspace(0, 2 * np.pi, 100) # 100 points from 0 to 2π
NumPy’s random module is essential for simulations, testing, and machine learning initialization:
# Uniform random values between 0 and 1 rand_uniform = np.random.rand(3, 3) print(rand_uniform) # Output: 3x3 matrix of random floats in [0, 1) # Standard normal distribution (mean=0, std=1) rand_normal = np.random.randn(3, 3) print(rand_normal) # Output: 3x3 matrix of values from normal distribution # Random integers rand_int = np.random.randint(1, 100, size=(2, 4)) print(rand_int) # Output: 2x4 matrix of random ints between 1 and 99 # Reproducible random numbers with seed np.random.seed(42) reproducible = np.random.rand(3) print(reproducible) # Output: [0.37454012 0.95071431 0.73199394] # Using the newer Generator API (recommended for new code) rng = np.random.default_rng(seed=42) values = rng.random(5) print(values) # Output: [0.77395605 0.43887844 0.85859792 0.69736803 0.09417735] # Random choice from an array choices = rng.choice([10, 20, 30, 40, 50], size=3, replace=False) print(choices) # Output: 3 random elements without replacement
Understanding array properties is essential for debugging and writing correct NumPy code. Every ndarray carries metadata about its structure:
import numpy as np
arr = np.array([[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12]])
# shape: dimensions as a tuple (rows, columns)
print(f"Shape: {arr.shape}")
# Output: Shape: (3, 4)
# ndim: number of dimensions (axes)
print(f"Dimensions: {arr.ndim}")
# Output: Dimensions: 2
# size: total number of elements
print(f"Total elements: {arr.size}")
# Output: Total elements: 12
# dtype: data type of elements
print(f"Data type: {arr.dtype}")
# Output: Data type: int64
# itemsize: size of each element in bytes
print(f"Bytes per element: {arr.itemsize}")
# Output: Bytes per element: 8
# nbytes: total memory consumed
print(f"Total bytes: {arr.nbytes}")
# Output: Total bytes: 96
# Practical example: understanding memory usage
large_arr = np.zeros((1000, 1000), dtype=np.float64)
print(f"Memory: {large_arr.nbytes / 1024 / 1024:.1f} MB")
# Output: Memory: 7.6 MB
# Same array with float32 uses half the memory
small_arr = np.zeros((1000, 1000), dtype=np.float32)
print(f"Memory: {small_arr.nbytes / 1024 / 1024:.1f} MB")
# Output: Memory: 3.8 MB
The dtype attribute is particularly important. NumPy supports many data types: int8, int16, int32, int64, float16, float32, float64, complex64, complex128, bool, and more. Choosing the right dtype can significantly impact both memory usage and computation speed.
NumPy’s indexing is more powerful than Python list indexing. Mastering it will save you from writing unnecessary loops.
arr = np.array([10, 20, 30, 40, 50, 60, 70, 80]) # Basic indexing (0-based) print(arr[0]) # 10 print(arr[-1]) # 80 print(arr[-2]) # 70 # Slicing: start:stop:step print(arr[2:5]) # [30 40 50] print(arr[:3]) # [10 20 30] print(arr[5:]) # [60 70 80] print(arr[::2]) # [10 30 50 70] — every other element print(arr[::-1]) # [80 70 60 50 40 30 20 10] — reversed
matrix = np.array([[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
[13, 14, 15, 16]])
# Single element: [row, col]
print(matrix[0, 0]) # 1
print(matrix[2, 3]) # 12
# Entire row
print(matrix[1]) # [5 6 7 8]
print(matrix[1, :]) # [5 6 7 8] — equivalent
# Entire column
print(matrix[:, 2]) # [ 3 7 11 15]
# Sub-matrix (rows 0-1, columns 1-2)
print(matrix[0:2, 1:3])
# Output:
# [[2 3]
# [6 7]]
# Every other row, every other column
print(matrix[::2, ::2])
# Output:
# [[ 1 3]
# [ 9 11]]
Boolean indexing is one of NumPy’s most powerful features. You create a boolean mask and use it to filter elements:
arr = np.array([15, 22, 8, 41, 3, 67, 29, 55]) # Elements greater than 20 mask = arr > 20 print(mask) # Output: [False True False True False True True True] print(arr[mask]) # Output: [22 41 67 29 55] # Shorthand — most common pattern print(arr[arr > 20]) # Output: [22 41 67 29 55] # Combining conditions (use & for AND, | for OR, ~ for NOT) print(arr[(arr > 10) & (arr < 50)]) # Output: [15 22 41 29] print(arr[(arr < 10) | (arr > 50)]) # Output: [ 8 3 67 55] # Boolean indexing on 2D arrays matrix = np.array([[1, 2], [3, 4], [5, 6]]) print(matrix[matrix % 2 == 0]) # Output: [2 4 6] — returns a flat array of even numbers
Fancy indexing lets you use arrays of indices to access multiple elements at once:
arr = np.array([10, 20, 30, 40, 50])
# Select elements at indices 0, 2, and 4
indices = np.array([0, 2, 4])
print(arr[indices])
# Output: [10 30 50]
# Works with 2D arrays too
matrix = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9],
[10, 11, 12]])
# Select specific rows
print(matrix[[0, 2, 3]])
# Output:
# [[ 1 2 3]
# [ 7 8 9]
# [10 11 12]]
# Select specific elements: (row0,col1), (row1,col2), (row2,col0)
rows = np.array([0, 1, 2])
cols = np.array([1, 2, 0])
print(matrix[rows, cols])
# Output: [2 6 7]
NumPy’s real power shows up in array operations. Everything is vectorized — no loops needed.
a = np.array([1, 2, 3, 4]) b = np.array([10, 20, 30, 40]) # Arithmetic is element-wise print(a + b) # [11 22 33 44] print(a - b) # [ -9 -18 -27 -36] print(a * b) # [ 10 40 90 160] print(b / a) # [10. 10. 10. 10.] print(a ** 2) # [ 1 4 9 16] # Comparison operators return boolean arrays print(a > 2) # [False False True True] print(a == b) # [False False False False] # Scalar operations are broadcast to every element print(a + 100) # [101 102 103 104] print(a * 3) # [ 3 6 9 12]
Broadcasting is the mechanism that lets NumPy perform operations on arrays of different shapes. It’s one of the most important concepts to understand:
# Broadcasting a scalar across an array
arr = np.array([[1, 2, 3],
[4, 5, 6]])
print(arr * 10)
# Output:
# [[10 20 30]
# [40 50 60]]
# Broadcasting a 1D array across rows of a 2D array
row = np.array([100, 200, 300])
print(arr + row)
# Output:
# [[101 202 303]
# [104 205 306]]
# Broadcasting a column vector across columns
col = np.array([[10],
[20]])
print(arr + col)
# Output:
# [[11 12 13]
# [24 25 26]]
# Practical example: centering data (subtracting column means)
data = np.array([[1.0, 200, 3000],
[2.0, 400, 6000],
[3.0, 600, 9000]])
col_means = data.mean(axis=0)
print(f"Column means: {col_means}")
# Output: Column means: [2.000e+00 4.000e+02 6.000e+03]
centered = data - col_means
print(centered)
# Output:
# [[-1.000e+00 -2.000e+02 -3.000e+03]
# [ 0.000e+00 0.000e+00 0.000e+00]
# [ 1.000e+00 2.000e+02 3.000e+03]]
Broadcasting rules:
ValueError.arr = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
# Global aggregations
print(f"Sum: {arr.sum()}") # 45
print(f"Mean: {arr.mean()}") # 5.0
print(f"Min: {arr.min()}") # 1
print(f"Max: {arr.max()}") # 9
print(f"Std Dev: {arr.std():.4f}") # 2.5820
# Aggregation along axes
# axis=0 → collapse rows (compute across rows → one value per column)
# axis=1 → collapse columns (compute across columns → one value per row)
print(f"Column sums: {arr.sum(axis=0)}") # [12 15 18]
print(f"Row sums: {arr.sum(axis=1)}") # [ 6 15 24]
print(f"Column means: {arr.mean(axis=0)}") # [4. 5. 6.]
print(f"Row means: {arr.mean(axis=1)}") # [2. 5. 8.]
# Other useful aggregations
print(f"Cumulative sum: {np.array([1,2,3,4]).cumsum()}")
# Output: [ 1 3 6 10]
print(f"Product: {np.array([1,2,3,4]).prod()}")
# Output: 24
# argmin and argmax — index of min/max value
scores = np.array([82, 91, 76, 95, 88])
print(f"Best score index: {scores.argmax()}") # 3
print(f"Worst score index: {scores.argmin()}") # 2
Reshaping lets you change the dimensions of an array without changing its data. This is critical when preparing data for machine learning models or matrix operations.
arr = np.arange(12) print(arr) # Output: [ 0 1 2 3 4 5 6 7 8 9 10 11] # Reshape to 3 rows × 4 columns reshaped = arr.reshape(3, 4) print(reshaped) # Output: # [[ 0 1 2 3] # [ 4 5 6 7] # [ 8 9 10 11]] # Reshape to 4 rows × 3 columns print(arr.reshape(4, 3)) # Output: # [[ 0 1 2] # [ 3 4 5] # [ 6 7 8] # [ 9 10 11]] # Use -1 to let NumPy infer one dimension print(arr.reshape(2, -1)) # 2 rows, auto-compute columns → (2, 6) print(arr.reshape(-1, 3)) # auto-compute rows, 3 columns → (4, 3) # Reshape to 3D print(arr.reshape(2, 2, 3).shape) # Output: (2, 2, 3) # IMPORTANT: total elements must match # arr.reshape(3, 5) # ValueError: cannot reshape array of size 12 into shape (3,5)
matrix = np.array([[1, 2, 3],
[4, 5, 6]])
# flatten() — always returns a copy
flat = matrix.flatten()
print(flat)
# Output: [1 2 3 4 5 6]
flat[0] = 999
print(matrix[0, 0]) # 1 — original unchanged (it's a copy)
# ravel() — returns a view when possible (more memory efficient)
raveled = matrix.ravel()
print(raveled)
# Output: [1 2 3 4 5 6]
raveled[0] = 999
print(matrix[0, 0]) # 999 — original IS changed (it's a view)
matrix = np.array([[1, 2, 3],
[4, 5, 6]])
print(f"Original shape: {matrix.shape}")
# Output: Original shape: (2, 3)
transposed = matrix.T
print(f"Transposed shape: {transposed.shape}")
# Output: Transposed shape: (3, 2)
print(transposed)
# Output:
# [[1 4]
# [2 5]
# [3 6]]
# np.transpose() and .T are equivalent for 2D arrays
# For higher dimensions, np.transpose() lets you specify axis order
arr_3d = np.arange(24).reshape(2, 3, 4)
print(arr_3d.shape) # (2, 3, 4)
print(np.transpose(arr_3d, (1, 0, 2)).shape) # (3, 2, 4)
Combining and dividing arrays is a common operation when preparing datasets or assembling results.
a = np.array([1, 2, 3]) b = np.array([4, 5, 6]) # Vertical stack — adds rows vs = np.vstack([a, b]) print(vs) # Output: # [[1 2 3] # [4 5 6]] # Horizontal stack — concatenates side by side hs = np.hstack([a, b]) print(hs) # Output: [1 2 3 4 5 6] # 2D stacking m1 = np.array([[1, 2], [3, 4]]) m2 = np.array([[5, 6], [7, 8]]) print(np.vstack([m1, m2])) # Output: # [[1 2] # [3 4] # [5 6] # [7 8]] print(np.hstack([m1, m2])) # Output: # [[1 2 5 6] # [3 4 7 8]] # np.concatenate — general purpose (specify axis) print(np.concatenate([m1, m2], axis=0)) # same as vstack print(np.concatenate([m1, m2], axis=1)) # same as hstack # Column stack — treats 1D arrays as columns c1 = np.array([1, 2, 3]) c2 = np.array([4, 5, 6]) print(np.column_stack([c1, c2])) # Output: # [[1 4] # [2 5] # [3 6]]
arr = np.arange(16).reshape(4, 4)
print(arr)
# Output:
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]
# [12 13 14 15]]
# Split into 2 equal parts along rows (axis=0)
top, bottom = np.vsplit(arr, 2)
print("Top:\n", top)
# Output:
# [[0 1 2 3]
# [4 5 6 7]]
print("Bottom:\n", bottom)
# Output:
# [[ 8 9 10 11]
# [12 13 14 15]]
# Split into 2 equal parts along columns (axis=1)
left, right = np.hsplit(arr, 2)
print("Left:\n", left)
# Output:
# [[ 0 1]
# [ 4 5]
# [ 8 9]
# [12 13]]
# Split at specific indices
first, second, third = np.split(arr, [1, 3], axis=0)
print(f"First (row 0): {first}")
print(f"Second (rows 1-2):\n{second}")
print(f"Third (row 3): {third}")
NumPy provides a comprehensive set of mathematical functions — all vectorized and optimized.
arr = np.array([1, 4, 9, 16, 25]) # Square root print(np.sqrt(arr)) # Output: [1. 2. 3. 4. 5.] # Exponential (e^x) print(np.exp(np.array([0, 1, 2]))) # Output: [1. 2.71828183 7.3890561 ] # Natural logarithm print(np.log(np.array([1, np.e, np.e**2]))) # Output: [0. 1. 2.] # Log base 10 and base 2 print(np.log10(np.array([1, 10, 100, 1000]))) # Output: [0. 1. 2. 3.] print(np.log2(np.array([1, 2, 4, 8]))) # Output: [0. 1. 2. 3.] # Trigonometric functions angles = np.array([0, np.pi/6, np.pi/4, np.pi/3, np.pi/2]) print(np.sin(angles)) # Output: [0. 0.5 0.70710678 0.8660254 1. ] print(np.cos(angles)) # Output: [1.00000000e+00 8.66025404e-01 7.07106781e-01 5.00000000e-01 6.12323400e-17] # Absolute value print(np.abs(np.array([-3, -1, 0, 2, 5]))) # Output: [3 1 0 2 5] # Rounding vals = np.array([1.23, 2.67, 3.5, 4.89]) print(np.round(vals, 1)) # [1.2 2.7 3.5 4.9] print(np.floor(vals)) # [1. 2. 3. 4.] print(np.ceil(vals)) # [2. 3. 4. 5.]
# Dot product of 1D arrays (scalar result)
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
print(np.dot(a, b))
# Output: 32 (1*4 + 2*5 + 3*6)
# Matrix multiplication
A = np.array([[1, 2],
[3, 4]])
B = np.array([[5, 6],
[7, 8]])
# Three equivalent ways to multiply matrices
print(np.dot(A, B))
print(A @ B) # @ operator (Python 3.5+)
print(np.matmul(A, B))
# All output:
# [[19 22]
# [43 50]]
# IMPORTANT: * is element-wise, NOT matrix multiplication
print(A * B)
# Output:
# [[ 5 12]
# [21 32]]
# Cross product
print(np.cross(np.array([1, 0, 0]), np.array([0, 1, 0])))
# Output: [0 0 1]
A = np.array([[1, 2],
[3, 4]])
# Determinant
print(f"Determinant: {np.linalg.det(A):.1f}")
# Output: Determinant: -2.0
# Inverse
A_inv = np.linalg.inv(A)
print(f"Inverse:\n{A_inv}")
# Output:
# [[-2. 1. ]
# [ 1.5 -0.5]]
# Verify: A × A_inv = Identity
print(np.round(A @ A_inv))
# Output:
# [[1. 0.]
# [0. 1.]]
# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)
print(f"Eigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")
# Matrix rank
print(f"Rank: {np.linalg.matrix_rank(A)}")
# Output: Rank: 2
# Norm
print(f"Frobenius norm: {np.linalg.norm(A):.4f}")
# Output: Frobenius norm: 5.4772
Understanding why NumPy is faster than Python lists is important for making good design decisions.
import numpy as np
import time
size = 1_000_000
# Python list approach
py_list = list(range(size))
start = time.time()
py_result = [x ** 2 for x in py_list]
py_time = time.time() - start
print(f"Python list: {py_time:.4f} seconds")
# NumPy approach
np_arr = np.arange(size)
start = time.time()
np_result = np_arr ** 2
np_time = time.time() - start
print(f"NumPy array: {np_time:.4f} seconds")
print(f"NumPy is {py_time / np_time:.0f}x faster")
# Typical output:
# Python list: 0.1654 seconds
# NumPy array: 0.0012 seconds
# NumPy is 138x faster
import sys
# Python list of 1000 integers
py_list = list(range(1000))
py_size = sys.getsizeof(py_list) + sum(sys.getsizeof(x) for x in py_list)
print(f"Python list: {py_size:,} bytes")
# NumPy array of 1000 integers
np_arr = np.arange(1000, dtype=np.int64)
print(f"NumPy array: {np_arr.nbytes:,} bytes")
print(f"Python list uses {py_size / np_arr.nbytes:.1f}x more memory")
# Typical output:
# Python list: 36,056 bytes
# NumPy array: 8,000 bytes
# Python list uses 4.5x more memory
Why is NumPy faster?
Digital images are just NumPy arrays. A grayscale image is a 2D array; a color image is 3D (height × width × channels).
import numpy as np
# Simulate a small 5x5 grayscale image (values 0-255)
image = np.array([
[50, 80, 120, 160, 200],
[55, 85, 125, 165, 205],
[60, 90, 130, 170, 210],
[65, 95, 135, 175, 215],
[70, 100, 140, 180, 220]
], dtype=np.uint8)
print(f"Image shape: {image.shape}")
print(f"Pixel value range: {image.min()} - {image.max()}")
# Invert the image (negative)
inverted = 255 - image
print(f"Inverted:\n{inverted}")
# Increase brightness (clamp to 255)
brightened = np.clip(image.astype(np.int16) + 50, 0, 255).astype(np.uint8)
print(f"Brightened:\n{brightened}")
# Threshold to binary (black/white)
threshold = 128
binary = (image > threshold).astype(np.uint8) * 255
print(f"Binary:\n{binary}")
# Normalize to [0, 1] range (common preprocessing step)
normalized = image.astype(np.float32) / 255.0
print(f"Normalized range: {normalized.min():.2f} - {normalized.max():.2f}")
# Simulate RGB image processing
rgb_image = np.random.randint(0, 256, size=(100, 100, 3), dtype=np.uint8)
print(f"RGB shape: {rgb_image.shape}") # (100, 100, 3)
# Convert to grayscale using weighted average
weights = np.array([0.2989, 0.5870, 0.1140]) # Standard luminance weights
grayscale = np.dot(rgb_image[...,:3], weights).astype(np.uint8)
print(f"Grayscale shape: {grayscale.shape}") # (100, 100)
import numpy as np
# Simulate exam scores for 5 subjects, 100 students
np.random.seed(42)
scores = np.random.normal(loc=72, scale=12, size=(100, 5))
scores = np.clip(scores, 0, 100).round(1)
subjects = ['Math', 'Science', 'English', 'History', 'Art']
print("=== Class Statistics ===\n")
# Per-subject statistics
for i, subject in enumerate(subjects):
col = scores[:, i]
print(f"{subject:>10}: mean={col.mean():.1f}, "
f"std={col.std():.1f}, "
f"min={col.min():.1f}, "
f"max={col.max():.1f}, "
f"median={np.median(col):.1f}")
print(f"\n{'Overall':>10}: mean={scores.mean():.1f}, std={scores.std():.1f}")
# Find top 5 students by average score
student_averages = scores.mean(axis=1)
top_5_indices = np.argsort(student_averages)[-5:][::-1]
print(f"\nTop 5 students (by index): {top_5_indices}")
for idx in top_5_indices:
print(f" Student {idx}: avg = {student_averages[idx]:.1f}")
# Correlation between subjects
correlation = np.corrcoef(scores.T)
print(f"\nCorrelation matrix shape: {correlation.shape}")
print(f"Math-Science correlation: {correlation[0, 1]:.3f}")
# Percentile analysis
print(f"\n90th percentile per subject:")
for i, subject in enumerate(subjects):
p90 = np.percentile(scores[:, i], 90)
print(f" {subject}: {p90:.1f}")
# Students scoring above 90 in all subjects
high_achievers = np.all(scores > 90, axis=1)
print(f"\nStudents scoring >90 in ALL subjects: {high_achievers.sum()}")
Solving systems of linear equations is a fundamental operation in engineering and data science. Consider:
import numpy as np
# Solve the system:
# 2x + 3y - z = 1
# 4x + y + 2z = 2
# -2x + 7y - 3z = -1
# Coefficient matrix
A = np.array([[2, 3, -1],
[4, 1, 2],
[-2, 7, -3]])
# Constants vector
b = np.array([1, 2, -1])
# Solve using np.linalg.solve (faster and more stable than computing inverse)
x = np.linalg.solve(A, b)
print(f"Solution: x={x[0]:.4f}, y={x[1]:.4f}, z={x[2]:.4f}")
# Verify the solution
residual = A @ x - b
print(f"Residual (should be ~0): {residual}")
print(f"Max error: {np.abs(residual).max():.2e}")
# Least squares solution for overdetermined systems
# (more equations than unknowns — common in data fitting)
# Fit y = mx + c to noisy data
np.random.seed(42)
x_data = np.linspace(0, 10, 50)
y_data = 2.5 * x_data + 1.3 + np.random.normal(0, 1, 50)
# Set up matrix A for y = mx + c
A_fit = np.column_stack([x_data, np.ones(len(x_data))])
# Solve via least squares
result, residuals, rank, sv = np.linalg.lstsq(A_fit, y_data, rcond=None)
m, c = result
print(f"\nLeast squares fit: y = {m:.4f}x + {c:.4f}")
print(f"(True values: y = 2.5000x + 1.3000)")
Normalization and standardization are essential preprocessing steps in machine learning. NumPy makes them trivial:
import numpy as np
# Sample dataset: 5 samples with 3 features of different scales
data = np.array([
[25.0, 50000, 3.5],
[30.0, 60000, 4.2],
[22.0, 45000, 3.1],
[35.0, 80000, 4.8],
[28.0, 55000, 3.9]
])
feature_names = ['Age', 'Salary', 'GPA']
print("Original data:")
print(data)
# Min-Max Normalization: scale to [0, 1]
min_vals = data.min(axis=0)
max_vals = data.max(axis=0)
normalized = (data - min_vals) / (max_vals - min_vals)
print(f"\nMin-Max Normalized (range [0, 1]):")
for i, name in enumerate(feature_names):
print(f" {name}: min={normalized[:, i].min():.2f}, max={normalized[:, i].max():.2f}")
print(normalized)
# Z-Score Standardization: mean=0, std=1
mean_vals = data.mean(axis=0)
std_vals = data.std(axis=0)
standardized = (data - mean_vals) / std_vals
print(f"\nZ-Score Standardized (mean≈0, std≈1):")
for i, name in enumerate(feature_names):
print(f" {name}: mean={standardized[:, i].mean():.4f}, std={standardized[:, i].std():.4f}")
print(standardized)
# Robust scaling (using median and IQR — resistant to outliers)
median_vals = np.median(data, axis=0)
q75 = np.percentile(data, 75, axis=0)
q25 = np.percentile(data, 25, axis=0)
iqr = q75 - q25
robust_scaled = (data - median_vals) / iqr
print(f"\nRobust Scaled (using median and IQR):")
print(robust_scaled)
Even experienced developers trip over these. Save yourself the debugging time.
This is the single most common source of bugs in NumPy code:
import numpy as np original = np.array([1, 2, 3, 4, 5]) # Slicing creates a VIEW, not a copy view = original[1:4] view[0] = 999 print(original) # Output: [ 1 999 3 4 5] — original is modified! # To create an independent copy, use .copy() original = np.array([1, 2, 3, 4, 5]) safe_copy = original[1:4].copy() safe_copy[0] = 999 print(original) # Output: [1 2 3 4 5] — original is safe # How to check: use np.shares_memory() a = np.array([1, 2, 3, 4, 5]) b = a[1:4] c = a[1:4].copy() print(np.shares_memory(a, b)) # True — b is a view print(np.shares_memory(a, c)) # False — c is a copy # Boolean and fancy indexing ALWAYS return copies d = a[a > 2] print(np.shares_memory(a, d)) # False
import numpy as np
a = np.array([[1, 2, 3],
[4, 5, 6]]) # shape (2, 3)
# This works — (3,) broadcasts to (2, 3)
row = np.array([10, 20, 30])
print(a + row)
# This FAILS — shapes (2, 3) and (2,) are incompatible
col_wrong = np.array([10, 20])
try:
print(a + col_wrong)
except ValueError as e:
print(f"Error: {e}")
# Error: operands could not be broadcast together with shapes (2,3) (2,)
# Fix: reshape to column vector (2, 1)
col_right = np.array([[10], [20]]) # shape (2, 1)
print(a + col_right)
# Output:
# [[11 12 13]
# [24 25 26]]
# Alternatively, use np.newaxis (or None — they're the same)
col_also_right = np.array([10, 20])[:, np.newaxis]
print(col_also_right.shape) # (2, 1)
print(a + col_also_right) # same result
import numpy as np # int8 can only hold values from -128 to 127 arr = np.array([100, 120, 130], dtype=np.int8) print(arr) # Output: [100 120 -126] — 130 overflowed silently! result = arr + np.int8(50) print(result) # Output: [-106 -86 -76] — completely wrong, no warning! # Fix: use a larger dtype arr_safe = np.array([100, 120, 130], dtype=np.int32) result_safe = arr_safe + 50 print(result_safe) # Output: [150 170 180] — correct # Watch out with uint8 (common for image data, range 0-255) img_pixel = np.array([250], dtype=np.uint8) print(img_pixel + np.uint8(10)) # Output: [4] — wrapped around! (250 + 10 = 260 → 260 % 256 = 4) # Fix: cast before arithmetic print(img_pixel.astype(np.int16) + 10) # Output: [260] — correct
import numpy as np
arr = np.array([[1, 2, 3],
[4, 5, 6]])
# DON'T: Chained indexing may not work for setting values
# arr[arr > 3][0] = 99 # This might NOT modify arr
# DO: Use direct indexing
arr[arr > 3] = 99
print(arr)
# Output:
# [[ 1 2 3]
# [99 99 99]]
# Or use np.where for conditional replacement
arr2 = np.array([[1, 2, 3],
[4, 5, 6]])
result = np.where(arr2 > 3, 99, arr2)
print(result)
# Output:
# [[ 1 2 3]
# [99 99 99]]
Follow these guidelines to write efficient, maintainable NumPy code.
import numpy as np
data = np.random.rand(1_000_000)
# BAD: Python loop
result_slow = np.empty(len(data))
for i in range(len(data)):
result_slow[i] = data[i] ** 2 + 2 * data[i] + 1
# GOOD: Vectorized operation (10-100x faster)
result_fast = data ** 2 + 2 * data + 1
# For custom functions, use np.vectorize (still not as fast as native ufuncs)
def custom_func(x):
if x > 0.5:
return x ** 2
else:
return 0
vectorized_func = np.vectorize(custom_func)
result = vectorized_func(data)
# BEST: Use np.where instead of vectorize
result_best = np.where(data > 0.5, data ** 2, 0)
import numpy as np # Use the smallest dtype that fits your data # Integers small_ints = np.array([1, 2, 3, 4], dtype=np.int8) # -128 to 127 medium_ints = np.array([1, 2, 3, 4], dtype=np.int32) # -2B to 2B big_ints = np.array([1, 2, 3, 4], dtype=np.int64) # default, but 2x memory # Floats — float32 is usually sufficient for ML weights = np.random.randn(1000, 1000).astype(np.float32) # 3.8 MB # vs np.float64 which would be 7.6 MB # Boolean arrays for masks mask = np.zeros(1000, dtype=np.bool_) # 1 byte per element vs 8 for int64
import numpy as np data = np.random.rand(1000, 3) means = data.mean(axis=0) # shape (3,) # BAD: manually tiling to match shapes means_tiled = np.tile(means, (1000, 1)) # creates unnecessary copy centered_slow = data - means_tiled # GOOD: let broadcasting handle it (no extra memory) centered_fast = data - means # (1000, 3) - (3,) → broadcasting
import numpy as np
n = 10000
# BAD: growing an array with append (copies entire array each time)
result = np.array([])
for i in range(n):
result = np.append(result, i ** 2)
# GOOD: preallocate and fill
result = np.empty(n)
for i in range(n):
result[i] = i ** 2
# BEST: vectorize completely
result = np.arange(n) ** 2
import numpy as np arr = np.random.rand(1_000_000) # Creates a new array (uses extra memory) arr = arr * 2 # In-place operation (modifies existing array, saves memory) arr *= 2 # NumPy also provides in-place functions np.multiply(arr, 2, out=arr) np.add(arr, 1, out=arr)
.copy() when you need independence.float32 instead of float64 halves memory usage. Watch out for integer overflow with small dtypes like int8 and uint8.np.linalg.solve() is faster and more numerically stable than computing matrix inverses manually.np.append() in a loop. Preallocate with np.empty() or np.zeros(), or better yet, vectorize the computation entirely.NumPy is one of those libraries where the investment in learning it well pays dividends across your entire Python career. The patterns and concepts here — vectorization, broadcasting, memory-aware programming — are transferable to GPU computing, distributed computing, and any high-performance numerical work.