NumPy Tutorial for Beginners (2025): Learn NumPy Step by Step

NumPy (Numerical Python) is the foundation of scientific computing in Python. It powers pandas, scikit-learn, TensorFlow, PyTorch, and virtually every data science and machine learning library. This tutorial takes you from zero to confident NumPy user — no prior experience required.

Why Learn NumPy?

NumPy is the single most important library for numerical computing in Python. Here's why it matters:

Installation

# With pip
pip install numpy

# With conda (recommended for data science)
conda install numpy

# Verify
python -c "import numpy as np; print(np.__version__)"
# 1.26.4 (or newer)

NumPy is almost always pre-installed in Jupyter, Google Colab, and Anaconda distributions.

1. Your First NumPy Array

import numpy as np

# Create from a Python list
a = np.array([1, 2, 3, 4, 5])
print(a)        # [1 2 3 4 5]
print(type(a))  # <class 'numpy.ndarray'>

# 2D array (matrix)
matrix = np.array([[1, 2, 3],
                   [4, 5, 6],
                   [7, 8, 9]])
print(matrix.shape)  # (3, 3)
print(matrix.ndim)   # 2
print(matrix.size)   # 9
print(matrix.dtype)  # int64

ndarray vs Python list

2. Creating Arrays

NumPy provides many functions to create arrays without manual typing:

import numpy as np

# Zeros and ones
np.zeros((3, 4))          # 3×4 matrix of 0.0
np.ones((2, 3))           # 2×3 matrix of 1.0
np.full((2, 2), 7)        # 2×2 matrix filled with 7

# Identity matrix
np.eye(3)                 # 3×3 identity matrix

# Ranges
np.arange(0, 10, 2)       # [0 2 4 6 8] — like range()
np.linspace(0, 1, 5)      # [0.   0.25 0.5  0.75 1. ] — evenly spaced

# Random
np.random.rand(3, 3)      # Uniform random [0, 1)
np.random.randn(3, 3)     # Standard normal (mean=0, std=1)
np.random.randint(0, 10, size=(3, 3))   # Random integers

# From existing data
np.array([1.5, 2.5, 3.5], dtype=np.float32)   # Explicit dtype
np.zeros_like(matrix)     # Same shape/dtype as matrix, filled with 0
np.ones_like(matrix)      # Same shape/dtype as matrix, filled with 1

Common dtypes

# Convert dtype
a = np.array([1, 2, 3])
a_float = a.astype(np.float64)
print(a_float.dtype)  # float64

3. Array Indexing and Slicing

1D Arrays

a = np.array([10, 20, 30, 40, 50])

a[0]      # 10   — first element
a[-1]     # 50   — last element
a[1:3]    # [20 30]  — slice (end exclusive)
a[:3]     # [10 20 30]
a[2:]     # [30 40 50]
a[::2]    # [10 30 50]  — every 2nd
a[::-1]   # [50 40 30 20 10]  — reverse

2D Arrays

m = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

m[0, 0]   # 1    — row 0, col 0
m[1, 2]   # 6    — row 1, col 2
m[0]      # [1 2 3]  — entire row 0
m[:, 1]   # [2 5 8]  — entire column 1
m[0:2, 1:3]  # [[2 3] [5 6]]  — submatrix
m[::2, ::2]  # [[1 3] [7 9]]  — every 2nd row and col

Boolean Indexing (Very Important)

a = np.array([10, 25, 3, 47, 8])

mask = a > 10
print(mask)     # [False  True False  True False]
print(a[mask])  # [25 47]

# One-liner
print(a[a > 10])   # [25 47]
print(a[a % 2 == 0])  # [10 8]

# 2D boolean indexing
m = np.array([[1, 2, 3], [4, 5, 6]])
print(m[m > 3])  # [4 5 6]

Fancy Indexing

a = np.array([10, 20, 30, 40, 50])

# Index with an array of indices
indices = [0, 2, 4]
print(a[indices])   # [10 30 50]

# 2D fancy indexing
m = np.array([[1, 2], [3, 4], [5, 6]])
rows = [0, 2]
cols = [1, 0]
print(m[rows, cols])  # [2 5]  — (row 0, col 1) and (row 2, col 0)

4. Array Shape Manipulation

Understanding shapes is critical for ML/data work:

a = np.arange(12)
print(a.shape)  # (12,)

# Reshape
b = a.reshape(3, 4)
print(b.shape)   # (3, 4)

c = a.reshape(2, 2, 3)
print(c.shape)   # (2, 2, 3)

# -1 means "figure it out"
d = a.reshape(3, -1)   # (3, 4)
e = a.reshape(-1, 6)   # (2, 6)

# Flatten (always returns a copy)
flat = b.flatten()     # shape (12,)

# Ravel (returns a view when possible)
flat2 = b.ravel()      # shape (12,)

# Transpose
m = np.array([[1, 2, 3], [4, 5, 6]])
print(m.shape)       # (2, 3)
print(m.T.shape)     # (3, 2)

# Add dimension
a = np.array([1, 2, 3])        # shape (3,)
print(a[np.newaxis, :].shape)  # (1, 3)  — row vector
print(a[:, np.newaxis].shape)  # (3, 1)  — column vector

# Squeeze (remove dimensions of size 1)
x = np.zeros((1, 3, 1))
print(np.squeeze(x).shape)  # (3,)

Stack and Split

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.vstack([a, b])          # [[1 2 3] [4 5 6]] — stack rows
np.hstack([a, b])          # [1 2 3 4 5 6]     — concatenate along axis 1
np.stack([a, b], axis=0)   # [[1 2 3] [4 5 6]] — new axis
np.stack([a, b], axis=1)   # [[1 4] [2 5] [3 6]]
np.concatenate([a, b])     # [1 2 3 4 5 6]

# Split
x = np.arange(9)
np.split(x, 3)             # [array([0,1,2]), array([3,4,5]), array([6,7,8])]
np.array_split(x, 4)       # Unequal split OK

5. Math Operations

Element-wise Operations

NumPy performs operations element-by-element by default:

a = np.array([1, 2, 3, 4])
b = np.array([10, 20, 30, 40])

a + b      # [11 22 33 44]
a - b      # [-9 -18 -27 -36]
a * b      # [10 40 90 160]
b / a      # [10. 10. 10. 10.]
b // a     # [10 10 10 10]  — floor division
a ** 2     # [1 4 9 16]
b % 3      # [1 2 0 1]

# Scalar operations (broadcast to entire array)
a * 2      # [2 4 6 8]
a + 100    # [101 102 103 104]
a > 2      # [False False True True]

Universal Functions (ufuncs)

Ufuncs are fast element-wise functions implemented in C:

a = np.array([1.0, 4.0, 9.0, 16.0])

np.sqrt(a)          # [1. 2. 3. 4.]
np.exp(a)           # [e^1 e^4 e^9 e^16]
np.log(a)           # natural log
np.log2(a)          # base-2 log
np.log10(a)         # base-10 log
np.sin(a)           # sine
np.cos(a)           # cosine
np.abs(np.array([-1, -2, 3]))  # [1 2 3]
np.floor(np.array([1.7, 2.3]))  # [1. 2.]
np.ceil(np.array([1.7, 2.3]))   # [2. 3.]
np.round(np.array([1.456]), 2)  # [1.46]

Aggregation Functions

a = np.array([[1, 2, 3],
              [4, 5, 6]])

np.sum(a)           # 21  — all elements
np.sum(a, axis=0)   # [5 7 9]  — sum each column
np.sum(a, axis=1)   # [6 15]   — sum each row

np.mean(a)          # 3.5
np.mean(a, axis=0)  # [2.5 3.5 4.5]

np.std(a)           # 1.707...
np.var(a)           # 2.916...

np.min(a)           # 1
np.max(a)           # 6
np.argmin(a)        # 0  — index of minimum (flattened)
np.argmax(a)        # 5  — index of maximum (flattened)

np.cumsum(a)        # [ 1  3  6 10 15 21] — cumulative sum
np.cumprod(a)       # [  1   2   6  24 120 720]

np.median(a)        # 3.5
np.percentile(a, 75)  # 5.25

6. Broadcasting

Broadcasting lets NumPy operate on arrays of different shapes — this is one of NumPy's most powerful features:

# Scalar broadcasts to entire array
a = np.array([1, 2, 3])
a + 10          # [11 12 13]

# 1D array broadcasts across 2D matrix rows
matrix = np.array([[1, 2, 3],
                   [4, 5, 6]])
row = np.array([10, 20, 30])

matrix + row    # [[11 22 33]
                #  [14 25 36]]

# Column vector broadcasts across columns
col = np.array([[100],
                [200]])

matrix + col    # [[101 102 103]
                #  [204 205 206]]

Broadcasting Rules

Two dimensions are compatible if:

1. They are equal, OR
2. One of them is 1

Shape A:  (3, 4)   →  (3, 4)
Shape B:     (4)   →  (1, 4)  [auto-padded]
Result:   (3, 4)

Shape A:  (3, 1)   →  (3, 1)
Shape B:  (1, 4)   →  (1, 4)
Result:   (3, 4)

# Practical: subtract mean from each column
data = np.array([[1, 2, 3],
                 [4, 5, 6],
                 [7, 8, 9]])
col_means = data.mean(axis=0)  # [4. 5. 6.]
centered = data - col_means     # [[-3,-3,-3], [0,0,0], [3,3,3]]

# Practical: add bias vector to batch of vectors
batch = np.random.randn(32, 128)  # 32 examples, 128 features
bias = np.random.randn(128)        # 128 biases
result = batch + bias               # (32, 128) — bias broadcast over batch

7. Linear Algebra

NumPy's linalg module covers all essential matrix math:

import numpy as np

A = np.array([[1, 2],
              [3, 4]])
B = np.array([[5, 6],
              [7, 8]])

# Matrix multiplication
np.dot(A, B)         # [[19 22] [43 50]]
A @ B                # Same — preferred in Python 3.5+

# Element-wise vs matrix multiply
A * B                # [[5 12] [21 32]]  — NOT matrix multiply

# Transpose
A.T                  # [[1 3] [2 4]]

# Determinant
np.linalg.det(A)     # -2.0

# Inverse
np.linalg.inv(A)     # [[-2. 1.] [1.5 -0.5]]

# Solve linear system Ax = b
b = np.array([1, 2])
x = np.linalg.solve(A, b)   # [-0.  0.5]

# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)

# Singular Value Decomposition
U, S, Vt = np.linalg.svd(A)

# Matrix norm
np.linalg.norm(A)           # Frobenius norm
np.linalg.norm(A, ord=2)    # Spectral norm (largest singular value)

Practical Linear Algebra: Least Squares

# Fit a line y = mx + c to noisy data
x = np.array([1, 2, 3, 4, 5], dtype=float)
y = 2 * x + 1 + np.random.randn(5) * 0.5

# Build design matrix
X = np.column_stack([x, np.ones_like(x)])  # [[1,1],[2,1],...]

# Least squares solution
coeffs, residuals, rank, sv = np.linalg.lstsq(X, y, rcond=None)
m, c = coeffs
print(f"y ≈ {m:.2f}x + {c:.2f}")

8. Sorting and Searching

a = np.array([3, 1, 4, 1, 5, 9, 2, 6])

np.sort(a)              # [1 1 2 3 4 5 6 9]  — returns sorted copy
a.sort()                # sorts in-place

np.argsort(a)           # indices that would sort the array
np.argsort(a)[::-1]     # indices for descending order

# Top-k elements
top3_indices = np.argsort(a)[-3:][::-1]  # indices of 3 largest

# 2D sorting
m = np.array([[3, 1], [0, 2]])
np.sort(m, axis=1)      # sort each row: [[1 3] [0 2]]
np.sort(m, axis=0)      # sort each column: [[0 1] [3 2]]

# Search
np.where(a > 3)         # indices where condition is True
np.where(a > 3, a, 0)   # a where True, 0 where False (ternary)

np.searchsorted([1, 3, 5, 7], 4)   # 2 — where to insert 4 to keep sorted

np.nonzero(a)           # indices of non-zero elements

9. Set Operations

a = np.array([1, 2, 3, 4])
b = np.array([3, 4, 5, 6])

np.union1d(a, b)        # [1 2 3 4 5 6]
np.intersect1d(a, b)    # [3 4]
np.setdiff1d(a, b)      # [1 2]  — in a but not b
np.setxor1d(a, b)       # [1 2 5 6]  — in one but not both
np.in1d(a, b)           # [False False True True]  — a elements in b?
np.isin(a, b)           # Same as in1d, newer API
np.unique(np.array([1, 2, 2, 3, 1]))   # [1 2 3]

10. Random Module

rng = np.random.default_rng(seed=42)  # Modern API — reproducible

rng.random((3, 3))          # Uniform [0, 1)
rng.standard_normal((3, 3)) # Standard normal N(0,1)
rng.integers(0, 10, (3, 3)) # Random integers [0, 10)
rng.choice([1, 2, 3, 4], size=5)            # With replacement
rng.choice([1, 2, 3, 4], size=3, replace=False)  # Without replacement

# Shuffle (in-place)
a = np.arange(10)
rng.shuffle(a)

# Distributions
rng.normal(loc=0, scale=1, size=1000)     # Gaussian
rng.uniform(low=0, high=10, size=1000)    # Uniform
rng.binomial(n=10, p=0.5, size=1000)     # Binomial
rng.poisson(lam=5, size=1000)             # Poisson
rng.exponential(scale=1, size=1000)       # Exponential

11. File I/O

import numpy as np

# Save and load binary (fast, NumPy format)
a = np.array([1, 2, 3])
np.save('array.npy', a)           # Single array
b = np.load('array.npy')          # Load it back

# Multiple arrays
np.savez('arrays.npz', x=a, y=a*2)
data = np.load('arrays.npz')
print(data['x'], data['y'])

# Compressed (larger arrays)
np.savez_compressed('arrays_compressed.npz', x=a)

# Text files (CSV)
matrix = np.array([[1, 2, 3], [4, 5, 6]])
np.savetxt('matrix.csv', matrix, delimiter=',', fmt='%d')
loaded = np.loadtxt('matrix.csv', delimiter=',')

# Load CSV with header (skip first row)
loaded = np.loadtxt('matrix.csv', delimiter=',', skiprows=1)

# genfromtxt — handles missing values
data = np.genfromtxt('data.csv', delimiter=',', filling_values=0)

3 Practical Projects

Project 1: Image as a NumPy Array

Images are just 3D NumPy arrays — this is the foundation of computer vision:

import numpy as np

# Simulate an RGB image (height=100, width=100, channels=3)
image = np.random.randint(0, 256, (100, 100, 3), dtype=np.uint8)

print(image.shape)    # (100, 100, 3)
print(image.dtype)    # uint8

# Crop (rows 10-50, cols 20-60)
crop = image[10:50, 20:60]
print(crop.shape)  # (40, 40, 3)

# Convert to grayscale (average channels)
grayscale = image.mean(axis=2)  # (100, 100)

# Extract red channel only
red = image[:, :, 0]

# Flip horizontally
flipped = image[:, ::-1, :]

# Normalize to [0, 1]
normalized = image / 255.0

# Adjust brightness (clamp to [0, 255])
brighter = np.clip(image.astype(int) + 50, 0, 255).astype(np.uint8)

# Batch of images (4 images, 100×100, RGB)
batch = np.random.randint(0, 256, (4, 100, 100, 3), dtype=np.uint8)
print(batch.shape)  # (4, 100, 100, 3)

# Normalize entire batch
batch_norm = batch / 255.0

Project 2: Statistics on Student Scores

import numpy as np

# Student scores: 50 students × 5 subjects
rng = np.random.default_rng(42)
scores = rng.integers(40, 101, size=(50, 5))
subjects = ["Math", "Science", "English", "History", "CS"]

# --- Summary Statistics ---
print("=== Class Summary ===")
for i, subj in enumerate(subjects):
    col = scores[:, i]
    print(f"{subj:10s}: mean={col.mean():.1f}  std={col.std():.1f}  "
          f"min={col.min()}  max={col.max()}")

# --- Student Averages ---
student_avgs = scores.mean(axis=1)    # Average per student
top5_idx = np.argsort(student_avgs)[-5:][::-1]
print(f"\nTop 5 students (indices): {top5_idx}")
print(f"Their averages: {student_avgs[top5_idx].round(1)}")

# --- Grade Distribution ---
letter_grades = np.select(
    [scores >= 90, scores >= 80, scores >= 70, scores >= 60],
    ['A', 'B', 'C', 'D'],
    default='F'
)
for grade in ['A', 'B', 'C', 'D', 'F']:
    count = np.sum(letter_grades == grade)
    print(f"Grade {grade}: {count} ({count/scores.size*100:.1f}%)")

# --- Pass/Fail ---
passing_score = 60
passed = (scores >= passing_score).all(axis=1)  # Pass ALL subjects
print(f"\nStudents passing all subjects: {passed.sum()}/{len(passed)}")

# --- Correlation between subjects ---
correlation = np.corrcoef(scores.T)  # 5×5 correlation matrix
print(f"\nMath–CS correlation: {correlation[0, 4]:.3f}")

Project 3: Neural Network Forward Pass from Scratch

import numpy as np

# Simple 2-layer neural network (manual — no PyTorch/TF)
rng = np.random.default_rng(42)

# Dimensions
input_dim, hidden_dim, output_dim = 784, 128, 10
batch_size = 32

# Random weights and biases
W1 = rng.standard_normal((input_dim, hidden_dim)) * 0.01
b1 = np.zeros((1, hidden_dim))
W2 = rng.standard_normal((hidden_dim, output_dim)) * 0.01
b2 = np.zeros((1, output_dim))

# Activation functions
def relu(z):
    return np.maximum(0, z)

def softmax(z):
    e = np.exp(z - z.max(axis=1, keepdims=True))  # subtract max for stability
    return e / e.sum(axis=1, keepdims=True)

# Forward pass
X = rng.standard_normal((batch_size, input_dim))   # input batch
y = rng.integers(0, output_dim, batch_size)          # true labels

Z1 = X @ W1 + b1          # (32, 128)
A1 = relu(Z1)              # (32, 128)
Z2 = A1 @ W2 + b2         # (32, 10)
probs = softmax(Z2)        # (32, 10) — probability distribution

# Cross-entropy loss
log_probs = np.log(probs[np.arange(batch_size), y] + 1e-9)
loss = -log_probs.mean()
print(f"Loss: {loss:.4f}")

# Predictions
predictions = probs.argmax(axis=1)
accuracy = (predictions == y).mean()
print(f"Accuracy (random init): {accuracy:.2%}")

NumPy vs Python Lists: Performance

import numpy as np
import time

n = 1_000_000

# Python list
lst = list(range(n))
start = time.time()
result = [x * 2 for x in lst]
python_time = time.time() - start

# NumPy array
arr = np.arange(n)
start = time.time()
result = arr * 2
numpy_time = time.time() - start

print(f"Python list: {python_time*1000:.1f}ms")
print(f"NumPy array: {numpy_time*1000:.1f}ms")
print(f"Speedup: {python_time/numpy_time:.0f}x")
# Typical output:
# Python list: 45.2ms
# NumPy array: 0.9ms
# Speedup: 50x

Common Mistakes

NumPy vs Related Libraries

Learning Path

Quick Reference

import numpy as np

# Creation
np.array([1,2,3])              # from list
np.zeros((m, n))               # m×n zeros
np.ones((m, n))                # m×n ones
np.eye(n)                      # n×n identity
np.arange(start, stop, step)   # range
np.linspace(start, stop, num)  # evenly spaced
np.random.default_rng(seed).random((m,n))  # random

# Info
a.shape    a.ndim    a.size    a.dtype

# Reshape
a.reshape(m, n)   a.flatten()   a.T   a[np.newaxis,:]

# Index
a[i, j]   a[i:j]   a[a > 0]   a[[0,2,4]]

# Math
a + b   a * b   a @ b   np.dot(a,b)
np.sum(a, axis=0)   np.mean(a)   np.std(a)
np.sort(a)   np.argsort(a)   np.where(cond, x, y)

# Linalg
np.linalg.det(A)   np.linalg.inv(A)
np.linalg.solve(A, b)   np.linalg.eig(A)

FAQ

Q: Is NumPy the same as MATLAB?
Similar concepts (matrices, vectorized ops, broadcasting), but NumPy is open source, Python-native, and widely used in industry. MATLAB is commercial and common in academia/engineering.

Q: When should I use NumPy vs pandas?
NumPy for homogeneous numerical data and matrix math. Pandas for tabular data with labeled columns, mixed types, and CSV/database workflows. Pandas is built on NumPy.

Q: Does NumPy use the GPU?
By default, no — NumPy uses CPU only. For GPU, use CuPy (drop-in NumPy replacement on CUDA), PyTorch tensors, or JAX arrays.

Q: What is a NumPy view vs copy?
A view shares memory with the original array — modifying the view modifies the original. Slices create views. Use .copy() to create an independent copy. a.base is None returns True if a owns its data.

Q: Why does NumPy use C order vs Fortran order?
C order (row-major, default) stores rows contiguously. Fortran order (column-major) stores columns. Operations along rows are faster in C order. Specify with order='C' or order='F'. Most NumPy operations default to C order.

Q: How do I check if two arrays are equal?
Use np.array_equal(a, b) for exact equality (returns a single bool). Use np.allclose(a, b, rtol=1e-5) for floating-point comparison with tolerance. The == operator returns an array of booleans.