Linear Algebra with NumPy — Core Concepts

Why linear algebra matters in Python

Almost every field that uses data — machine learning, physics simulations, computer graphics, economics — reduces its core problems to linear algebra. NumPy’s linalg module gives Python developers direct access to the same battle-tested routines (LAPACK, BLAS) that power MATLAB and Fortran scientific code.

Vectors and matrices

A vector is a one-dimensional array of numbers. A matrix is two-dimensional. NumPy represents both as ndarray objects.

Key vector operations:

• Dot product — np.dot(a, b) or a @ b — measures how aligned two vectors are. Returns a scalar.
• Norm — np.linalg.norm(v) — the “length” of a vector. Euclidean norm by default.
• Cross product — np.cross(a, b) — for 3D vectors, returns a vector perpendicular to both inputs.

Solving systems of linear equations

A system like:

2x + 3y = 8
 x -  y = 1

Can be written as Ax = b, where A is the coefficient matrix, x is the unknown vector, and b is the result vector. NumPy solves it directly:

A = np.array([[2, 3], [1, -1]])
b = np.array([8, 1])
x = np.linalg.solve(A, b)  # array([2.2, 1.2])

This is faster and more numerically stable than computing inv(A) @ b.

Key decompositions

Decompositions break a matrix into simpler pieces. Three are essential:

LU decomposition

Factors A into a lower-triangular matrix L and upper-triangular matrix U. Used internally by solve for general systems. Available via SciPy:

from scipy.linalg import lu
P, L, U = lu(A)

QR decomposition

Factors A into an orthogonal matrix Q and upper-triangular R. Used for least-squares problems and stable numerical algorithms:

Q, R = np.linalg.qr(A)

Eigendecomposition

Finds eigenvalues λ and eigenvectors v such that Av = λv. Reveals the “natural axes” of a transformation:

eigenvalues, eigenvectors = np.linalg.eig(A)

For symmetric matrices, use eigh — it is faster and guarantees real results.

Singular Value Decomposition (SVD)

SVD is arguably the most important decomposition. It works on any matrix (not just square ones) and factors A into three matrices: U, Σ, and Vᵀ.

Applications include:

• PCA — principal component analysis for dimensionality reduction
• Image compression — keep only the top k singular values
• Recommender systems — matrix factorization for collaborative filtering
• Pseudoinverse — solving least-squares problems for non-square systems

Norms and condition numbers

The norm of a matrix measures its “size” in a specific sense. The condition number (np.linalg.cond(A)) measures how sensitive solutions are to input perturbations. A high condition number means small errors in your data produce large errors in results — the matrix is “ill-conditioned.”

Common misconception

People often assume that computing the matrix inverse is the right way to solve linear systems. In practice, np.linalg.solve is almost always better — it is faster (O(n³) vs the same for inverse plus an extra multiplication) and introduces fewer floating-point errors. Computing an explicit inverse should be reserved for cases where you genuinely need the inverse matrix itself.

One thing to remember: NumPy’s linalg.solve should be your first instinct for linear systems — it is faster and more stable than computing the inverse.