Linear Algebra in Action with NumPy

Why Linear Algebra Matters for ML

Nearly every machine learning model is built on linear algebra. Data is stored as vectors and matrices, model parameters are vectors, and training involves matrix operations. NumPy gives us fast, vectorized tools to do this math efficiently in Python.

Vectors and Matrices

A vector is a 1-D array of numbers; a matrix is a 2-D array. In NumPy both are ndarray objects.

import numpy as np
v = np.array([1, 2, 3])          # shape (3,)
M = np.array([[1, 2], [3, 4]])    # shape (2, 2)
print(v.shape, M.shape)

The .shape attribute is your best friend—most bugs are shape mismatches.

Dot Products and Norms

The dot product sums element-wise products and measures alignment between vectors:

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
print(np.dot(a, b))   # 32
print(a @ b)          # 32 (preferred operator)

The norm (length/magnitude) of a vector:

np.linalg.norm(a)     # sqrt(1+4+9) = 3.742

Cosine similarity measures the angle between vectors (used in recommendations and NLP):

cos = (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))

Matrix Multiplication

Matrix multiplication combines rows and columns. For A @ B, the inner dimensions must match: (m, n) @ (n, p) -> (m, p).

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print(A @ B)          # matrix product
print(A * B)          # element-wise (NOT the same!)
print(A.T)            # transpose

This is exactly how a neural network layer computes outputs = inputs @ weights + bias.

Broadcasting

Broadcasting lets NumPy operate on arrays of different shapes without explicit loops, by stretching smaller arrays.

X = np.array([[1, 2, 3], [4, 5, 6]])
mean = X.mean(axis=0)        # shape (3,)
X_centered = X - mean        # broadcasts row-wise

This is the foundation of feature normalization.

Solving Linear Systems

Linear regression’s closed-form solution uses matrix inversion. Prefer np.linalg.solve over computing the inverse directly for stability.

# Solve Ax = b
A = np.array([[3, 1], [1, 2]])
b = np.array([9, 8])
x = np.linalg.solve(A, b)    # [2., 3.]

Matrix Decompositions

Decompositions break a matrix into useful factors.

Eigendecomposition finds vectors that only get scaled (not rotated) by a matrix. It powers PCA on covariance matrices.

values, vectors = np.linalg.eig(A)

Singular Value Decomposition (SVD) factors any matrix A = U S Vᵀ. It underpins PCA, recommender systems, and data compression.

Xc = X - X.mean(axis=0)          # PCA requires centered data
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
# columns of Vt.T are principal directions; S**2/(n-1) are the explained variances

The largest singular values capture the most variance—keeping only the top ones gives dimensionality reduction. Note: the variance interpretation only holds after mean-centering (sklearn’s PCA centers the data internally for exactly this reason).

Putting It Together: Linear Regression

# Least squares with an intercept (no need to form XᵀX explicitly)
Xb = np.hstack([np.ones((X.shape[0], 1)), X])   # add intercept
w, *_ = np.linalg.lstsq(Xb, y, rcond=None)       # QR/SVD-based; more stable than solving Xt@X

Forming XᵀX squares the condition number; lstsq (QR/SVD) is numerically preferable. This trains a linear regression model—pure linear algebra.

Key Takeaways

• Always check shapes; use @ for matrix multiply, * for element-wise.
• Use np.linalg.solve instead of explicit inverses.
• SVD and eigendecomposition are the math behind PCA and many ML methods.