Linear algebra basics with NumPy

Linear algebra sits under almost every data science technique: regression, PCA, recommendation systems and neural networks are all matrix operations. NumPy's @ operator and the np.linalg module cover the essentials.

Vectors and dot products

import numpy as np

qty = np.array([2, 1, 3])            # kg of three sweets
price = np.array([600, 850, 480])    # ₹ per kg
bill = qty @ price                    # dot product: 2*600 + 1*850 + 3*480
print(bill)                           # 3490

The dot product multiplies matching elements and sums them. It is a weighted total, which is why it appears in every linear model.

Matrix multiplication

For A @ B, the number of columns in A must equal the number of rows in B. The result has shape (rows of A, columns of B).

# 3 orders x 3 products (kg)
orders = np.array([[2, 1, 3],
                   [0, 2, 1],
                   [1, 1, 1]])
# 3 products x 2 columns: price and cost per kg
rates = np.array([[600, 420],
                  [850, 610],
                  [480, 300]])

result = orders @ rates      # (3, 3) @ (3, 2) -> (3, 2)
print(result)                # each row: [revenue, cost] for one order

Note that orders * rates would be element-wise multiplication and fails here because the shapes differ. Use @ (or np.matmul) for matrix products.

Solving a system of equations

Three shops bought mixed boxes. How much did each item cost?

2x + 1y + 1z = 1550
1x + 3y + 2z = 2350
1x + 1y + 4z = 2050

A = np.array([[2, 1, 1],
              [1, 3, 2],
              [1, 1, 4]])
b = np.array([1550, 2350, 2050])

x = np.linalg.solve(A, b)
print(x)                    # [400. 450. 300.]
print(np.allclose(A @ x, b))   # True

Prefer np.linalg.solve(A, b) over np.linalg.inv(A) @ b. It is faster and more numerically accurate.

Other useful functions

Function	Returns
np.linalg.det(A)	Determinant (near 0 means singular, no unique solution)
np.linalg.inv(A)	Inverse matrix
np.linalg.norm(v)	Length of a vector (Euclidean by default)
np.linalg.eig(A)	Eigenvalues and eigenvectors
np.linalg.matrix_rank(A)	Number of independent rows/columns
np.linalg.lstsq(X, y)	Least-squares solution for over-determined systems

Least squares: a regression from scratch

Fit sales = a * ad_spend + b from noisy data:

rng = np.random.default_rng(1)
ad = rng.uniform(10, 100, size=50)                # ₹ thousand
sales = 3.2 * ad + 40 + rng.normal(0, 15, size=50)

X = np.column_stack([ad, np.ones_like(ad)])       # add intercept column
coef, *_ = np.linalg.lstsq(X, sales, rcond=None)
print(coef.round(2))                               # close to [3.2, 40]

This is exactly what linear regression libraries do internally.

Common mistakes

• Using * when you meant matrix multiplication @.
• Inverting matrices just to solve a system; use solve.
• Ignoring a near-zero determinant or a LinAlgError: Singular matrix, which means the equations are not independent.
• Forgetting the column of ones for the intercept in least squares.