Linear Algebra: A Practical Guide to Vector Spaces, Factorizations, and Numerical Computing

A vector space is a collection whose elements can be added and multiplied by scalars while satisfying rules such as associativity and distributivity. Real coordinate vectors are familiar examples, but polynomials and suitable collections of functions also qualify. A subspace must contain zero and remain closed under addition and scalar multiplication. For example, the plane x + y + z = 0 is a subspace of R³; the plane x + y + z = 1 is not, because it excludes zero.

A basis is a linearly independent set that spans the space, so every vector has a unique coordinate representation in that basis. Consider u = (1, 1) and v = (1, −1). They form a basis of R² because neither is a scalar multiple of the other. Solving a + b = 4 and a − b = 2 gives (4, 2) = 3u + v. Thus the same vector has standard coordinates (4, 2) and coordinates (3, 1) in the new basis. This distinction matters whenever practitioners rotate axes, change physical reference frames, or transform features.

A transformation T is linear when T(ax + by) = aT(x) + bT(y). After choosing bases, an m-by-n matrix A represents a map from Rⁿ to Rᵐ. Its columns are the images of the input basis vectors. Matrix multiplication represents composition: ABx applies B first, then A. Order therefore matters, and AB generally differs from BA. An affine transformation Ax + b is not linear when b is nonzero, although it is common in regression and neural network layers.

The rank-nullity theorem states that rank(A) + dim ker(A) = n, where n is the number of columns. For A = [[1, 2, 3], [2, 4, 6]], the second row duplicates twice the first, so the rank is one. Solving Ax = 0 gives x₁ = −2x₂ − 3x₃, making (−2, 1, 0) and (−3, 0, 1) a basis for the kernel. Its dimension is two, and 1 + 2 = 3 verifies the theorem. The system Ax = (5, 10) has infinitely many solutions, whereas Ax = (5, 11) has none. Rank therefore reveals both lost information and constraints on attainable outputs.

For a square, nonsingular system Ax = b, use a linear solver rather than explicitly computing A⁻¹. A common approach is LU factorization with partial pivoting, written PA = LU: P reorders rows, L is lower triangular, and U is upper triangular. Forward and backward substitution then solve the triangular systems. If A is symmetric positive definite, Cholesky factorization A = LLᵀ exploits that structure with less work and storage. When many right-hand sides share A, factor it once and reuse the factors.

As a small example, solve 3x + y = 7 and x + 2y = 5. Subtracting one-third of the first equation from the second gives (5/3)y = 8/3, so y = 1.6 and x = 1.8. In Python, np.linalg.solve(A, b) performs the corresponding solve without forming an inverse. However, a small residual r = b − Ax does not guarantee a small solution error: an ill-conditioned matrix can amplify tiny perturbations. The 2-norm condition number of a nonsingular matrix is σmax/σmin. Report a scaled residual and consider conditioning, data uncertainty, and precision before accepting a solution.

When measurements produce more equations than unknowns, an exact solution may not exist. Least squares chooses x to minimize ||Ax − b||₂². Geometrically, Ax is the orthogonal projection of b onto the column space of A, and the residual is perpendicular to every column. This yields the normal equations AᵀAx = Aᵀb. They are useful for derivation, but explicitly forming AᵀA squares the 2-norm condition number when A has full column rank. QR factorization or singular value decomposition is usually preferable when numerical accuracy matters.

Fit a line y = c + mt to observations (0, 1), (1, 2), and (2, 2). Use A = [[1, 0], [1, 1], [1, 2]] and b = (1, 2, 2). The normal equations become 3c + 3m = 5 and 3c + 5m = 6, giving m = 1/2 and c = 7/6. The residual b − Ax is (−1/6, 1/3, −1/6), with squared norm 1/6. Its dot product with each column of A is zero, confirming orthogonality. In practice, np.linalg.lstsq(A, b, rcond=None) also returns rank and singular values, which help diagnose dependent predictors.

Eigenvalues and eigenvectors describe directions that a square linear transformation scales without changing their span: Av = λv for nonzero v. For A = [[2, 1], [1, 2]], the vectors (1, 1) and (1, −1) have eigenvalues 3 and 1. After normalization, they form the columns of an orthogonal matrix Q, and A = QΛQᵀ with Λ = diag(3, 1). In those coordinates, applying A simply multiplies one coordinate by three and leaves the other unchanged. Repeated applications scale them by 3ᵏ and 1ᵏ.

The spectral theorem guarantees an orthonormal eigenbasis for real symmetric matrices; complex Hermitian matrices have the analogous result using conjugate transpose. It does not say that every square matrix is diagonalizable. For example, [[1, 1], [0, 1]] has only one independent eigenvector. In engineering, symmetric eigenproblems appear in energy models and vibration analysis, sometimes as generalized problems Kv = λMv. Use structure-aware routines such as np.linalg.eigh for symmetric or Hermitian inputs. For a covariance matrix, eigenvectors identify orthogonal directions of variation, while eigenvalues measure variance along those directions.

Singular value decomposition applies to rectangular as well as square matrices. For real A, write A = UΣVᵀ, with orthogonal factors in the full decomposition and nonnegative singular values in descending order. The right singular vectors identify input directions; the left singular vectors identify corresponding output directions. Keeping the first k components gives Aₖ = UₖΣₖVₖᵀ, a best approximation of rank at most k in both spectral and Frobenius norms. For A = diag(4, 1), the rank-one approximation diag(4, 0) has error one in either norm. The discarded singular values quantify the information lost.

Principal component analysis starts with a data matrix X whose rows are observations and columns are features. Subtract each column’s training-set mean, then compute X = UΣVᵀ. The reduced coordinates are XVₖ = UₖΣₖ. With n observations and sample covariance XᵀX/(n − 1), component variances are σᵢ²/(n − 1). Choose k using explained variance, validation performance, and domain constraints rather than a universal threshold. Standardizing features changes the optimization problem and can help when units differ. Fit centering and scaling only on training data, then reuse those parameters on validation and test sets to avoid leakage.

Numerical linear algebra requires matching algorithms to matrix structure and hardware. NumPy and SciPy support dense and sparse workflows; PyTorch and JAX integrate matrix operations with accelerator execution and automatic differentiation. Prefer vectorized products and batched factorizations over Python loops when the workload supports them. Yet accelerators are not automatically faster: transfer costs, matrix sizes, and precision requirements matter. For large sparse systems, use sparse storage and suitable iterative numerical methods, such as conjugate gradients for symmetric positive definite matrices, often with preconditioning.

A reliable workflow checks shapes, units, dtype, finiteness, symmetry, and scale before solving. Afterward, test residuals and relevant structural properties: QᵀQ should approximate I, and reconstructed factors should approximate the original matrix. Numerical rank requires a tolerance because floating-point computations rarely produce exact zeros; choose it with matrix scale and measurement noise in mind. Benchmark representative workloads and avoid unnecessary dense covariance matrices or full decompositions when reduced forms suffice. This combination of theorem-based reasoning, diagnostics, and implementation discipline is the bridge between classroom calculations and production-grade scientific or machine learning systems.