A machine-learning dataset is often a table: rows are observations and columns are features. A prediction uses several measurements from each row, and a training calculation uses many rows together. Linear algebra gives names to these objects and rules for combining them. It also explains why two pieces of code that both run can nevertheless represent different calculations.
This chapter develops the operations needed to understand regression. We begin with weighted sums, connect them to dot products and matrix multiplication, then examine systems of equations, rank, and projection. The small examples are constructed numerical demonstrations; they do not evaluate a learned model’s performance on unseen observations.
Prerequisites:Chapter 5 and introductory NumPy arrays. No calculus is required for the geometric derivation of least squares here.
Learning objectives:
State the dimensions and units of vectors and data matrices.
Compute dot products, norms, matrix products, and projections.
Explain linear independence, rank, and non-unique parameter values.
Distinguish exact linear solving from approximate least squares.
Recognize why a numerical solver is preferable to an explicitly computed inverse.
6.2 Scalars, vectors, and matrices
A scalar is one number. A vector is an ordered collection of numbers, and a matrix is a rectangular collection with rows and columns. Order is essential: the vector \((\text{speed},\text{cadence})\) must not silently become \((\text{cadence},\text{speed})\) when supplied to a fitted model.
For \(n\) observations with \(p\) features, define
The notation \(\mathbb{R}^{n\times p}\) means real-valued matrices with \(n\) rows and \(p\) columns. Entry \(x_{ij}\) is feature \(j\) for observation \(i\). The target vector is \(\mathbf y\in\mathbb R^n\), and a coefficient vector is \(\mathbf w\in\mathbb R^p\). Bold symbols indicate vectors; uppercase \(X\) names a matrix. These conventions describe the mathematics, while the Python examples use descriptive names.
Mathematical vectors are treated as columns unless stated otherwise. A NumPy array with shape (p,) is one-dimensional: it has neither a row axis nor a column axis. Transposing it does not create a (p, 1) column. A two-dimensional array with shape (p, 1) explicitly represents a column. Both conventions can be used correctly, but mixing them can trigger unintended broadcasting.
The transpose\(X^\top\) exchanges rows and columns, so it has shape \(p\times n\). Entry \((j,i)\) of \(X^\top\) equals entry \((i,j)\) of \(X\). Transposition reorganizes indices; it does not invert the numerical operation performed by a matrix.
6.3 Dot products as weighted sums
For vectors \(\mathbf a,\mathbf b\in\mathbb R^p\), their dot product is
First multiply corresponding entries, then sum the products. The output is a scalar. For \(\mathbf a=(2,3)\) and \(\mathbf b=(4,-1)\), the result is \(2(4)+3(-1)=5\). A linear model’s prediction is exactly this kind of weighted sum: \(\widehat y_i=\mathbf x_i^\top\mathbf w+b\), where \(\mathbf x_i\) contains the features of observation \(i\).
If one feature is speed and another is stride duration, their coefficients must convert those separate units into target units. A raw geometric distance between mixed physical units has no automatic interpretation. Scaling can make comparisons meaningful for a chosen task, but scaling does not remove the need to explain what the representation means.
For \((3,4)\) the length is 5. The \(\ell_1\) norm is \(\|\mathbf a\|_1=\sum_j|a_j|\), which is 7 for the same vector. The subscript names the norm definition, not an observation index. The Euclidean distance between vectors is \(\|\mathbf a-\mathbf b\|_2\); the difference must represent comparable features in the same order.
Two vectors are orthogonal if their dot product is zero. For example, \((1,1)\) and \((1,-1)\) are orthogonal. For nonzero vectors, the identity \(\mathbf a^\top\mathbf b=\|\mathbf a\|_2\|\mathbf b\|_2\cos\theta\) relates the dot product to angle \(\theta\). The angle interpretation assumes ordinary Euclidean coordinates; it is not a substitute for a sensible feature representation.
6.4 Matrix multiplication combines many weighted sums
If \(A\) has shape \(m\times p\) and \(B\) has shape \(p\times q\), then \(AB\) has shape \(m\times q\) with entries
\[
(AB)_{ik}=\sum_{j=1}^p A_{ij}B_{jk}.
\]
The inner dimensions must match because each row of \(A\) is paired with a column of \(B\). In particular, \(X\mathbf w\) has \(n\) entries: one prediction contribution per observation. Matrix multiplication generally depends on order; \(AB\) and \(BA\) may have different shapes or one may be undefined. It differs from elementwise multiplication, which pairs entries without summing.
An intercept can be absorbed into the matrix representation. Add a leading column of ones to \(X\), call the result \(D\in\mathbb R^{n\times(p+1)}\), and define \(\boldsymbol\beta=(b,w_1,\ldots,w_p)^\top\). Then \(\widehat{\mathbf y}=D\boldsymbol\beta\). The column of ones ensures that every row receives the same intercept. We use \(D\) for this augmented design matrix so the distinction from the original feature matrix remains explicit.
6.4.1 From one dot product to all predictions
Use four hypothetical observations with two dimensionless features. Choose known coefficients so that the arithmetic, array dimensions, and intercept can be checked independently.
Store one row per observation and one column per feature.
Compute one weighted sum with an explicit loop.
Use matrix multiplication to calculate every prediction.
Add a ones column and verify the equivalent augmented representation.
# Step 1: Import the numerical tools used throughout this chapter.import numpy as np# Store four artificial observations with two dimensionless features each.feature_matrix = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [2.0, 1.0]])# The first coefficient multiplies feature 1; the second multiplies feature 2.feature_weights = np.array([2.0, -1.0])# Add the same offset, in target units, to every prediction.intercept =3.0# Step 2: Expand the first row's dot product one term at a time.# Accumulate only the weighted feature contribution before adding the intercept.first_weighted_sum =0.0for feature_index inrange(feature_matrix.shape[1]):# Match this feature with the coefficient in the same position. feature_contribution = ( feature_matrix[0, feature_index] * feature_weights[feature_index] ) first_weighted_sum = first_weighted_sum + feature_contribution# Step 3: The @ operator computes all row-by-coefficient dot products.# The resulting one-dimensional array has one target prediction per row.predicted_targets = feature_matrix @ feature_weights + interceptassert np.isclose(predicted_targets[0], first_weighted_sum + intercept)assert np.allclose(predicted_targets, [5.0, 2.0, 4.0, 6.0])# Step 4: Represent the intercept by adding a column of ones.# Each one supplies the same intercept multiplier for its observation.intercept_column = np.ones(feature_matrix.shape[0])# np.column_stack places these one-dimensional arrays beside existing columns.design_matrix = np.column_stack((intercept_column, feature_matrix))# Coefficient order must match the design matrix's intercept-first column order.model_coefficients = np.array([intercept, 2.0, -1.0])assert design_matrix.shape == (4, 3)assert np.allclose(design_matrix @ model_coefficients, predicted_targets)print("Design shape:", design_matrix.shape)print("Predictions:", predicted_targets)
Design shape: (4, 3)
Predictions: [5. 2. 4. 6.]
The addition of a scalar intercept to a one-dimensional prediction array is an intentional NumPy broadcast. By contrast, subtracting an (n, 1) target array from an (n,) prediction array would produce an (n, n) result. The program may run while computing a different quantity. Check shapes before arithmetic, especially when calculating residuals and losses.
6.5 Linear transformations, span, and rank
A matrix transforms a vector by forming weighted combinations of its columns. If \(A=[\mathbf a_1\ \cdots\ \mathbf a_p]\), then \(A\mathbf w=w_1\mathbf a_1+\cdots+w_p\mathbf a_p\). The column space is the set of all such combinations. The span of a set of vectors is the set of all linear combinations they can produce. This viewpoint changes the question “which coefficient values fit?” into “which output vectors are reachable?”
Columns are linearly independent when the only combination producing the zero vector uses all-zero coefficients. The rank is the number of independent column directions, equivalently the dimension of the column space. For an \(n\times p\) matrix it cannot exceed \(\min(n,p)\).
Suppose one feature column is exactly twice another. Its coefficient can be traded against the first without changing predictions. The matrix is rank-deficient relative to its column count, and coefficients need not be unique. More data with the same exact column relationship cannot solve this identifiability problem. Removing redundancy, imposing a justified constraint, or adding regularization changes the situation, but these choices have different interpretations.
Numerical rank also involves a tolerance. Nearly dependent columns are not exactly dependent mathematically, yet small perturbations may produce large coefficient changes. The condition number describes sensitivity for a specified norm; a large value signals potential amplification of error. It does not prove that every prediction is unstable, but it encourages us to inspect scaling, representation, and coefficient interpretation carefully.
Numerical solvers often report singular values: nonnegative factors describing how strongly a matrix stretches independent orthogonal directions. A zero singular value indicates a lost direction, and a very small one indicates a weakly determined direction. We do not need the full singular-value factorization to use these diagnostic outputs responsibly.
6.6 Exact systems and least squares
An exact linear system has the form \(A\mathbf z=\mathbf c\). When \(A\) is square and nonsingular, there is exactly one solution. For example,
The second equation gives \(z_2=z_1-1\). Substitute into the first: \(2z_1+z_1-1=5\), so \(z_1=2\) and \(z_2=1\). Substitution into both original equations checks the answer. For a general data table with more observations than coefficients, an exact fit may not exist; noisy observations need not lie in the design matrix’s column space.
The notation \(\operatorname*{arg\,min}\) means the coefficient vector or vectors attaining the smallest objective value, rather than the minimum number itself. Squared residual norm is the sum of squared errors. Dividing by fixed \(n\) would not change the minimizing coefficients.
NumPy’s solve documentation specifies a square, full-rank system. Its lstsq documentation covers least squares and the minimum-norm solution when minimizers are not unique. The examples below provide rcond explicitly rather than depending on a numerical cutoff default.
6.6.1 Why projection solves least squares
Start with a single nonzero direction \(\mathbf u\) and a target vector \(\mathbf y\). We seek a multiple \(a\mathbf u\) close to \(\mathbf y\). At the closest point, the residual \(\mathbf r=\mathbf y-a\mathbf u\) is orthogonal to \(\mathbf u\). Otherwise, a small movement along \(\mathbf u\) could reduce the remaining component in that direction. Setting the dot product to zero gives
The projection is \(\operatorname{proj}_{\mathbf u}\mathbf y=a\mathbf u\). For \(\mathbf u=(1,1)\) and \(\mathbf y=(2,0)\), \(a=2/2=1\), so the projection is \((1,1)\) and the residual is \((1,-1)\). Their dot product is zero.
The same idea extends to all columns of \(D\): the best residual is orthogonal to every reachable direction. Thus
\[
D^\top(\mathbf y-D\widehat{\boldsymbol\beta})=\mathbf0,
\qquad
D^\top D\widehat{\boldsymbol\beta}=D^\top\mathbf y.
\]
These are the normal equations. If \(D\) has full column rank, \(D^\top D\) is invertible and they specify one coefficient vector. If rank is deficient, the fitted projection is still unique, while multiple coefficient vectors can represent it. The normal equations explain the geometry; forming an inverse is not the recommended computational procedure.
6.6.2 Solving and checking the geometry
# Step 5: Solve the hand-checkable square system without forming an inverse.# Rows store the two equations, and columns store coefficients of z_1 and z_2.system_matrix = np.array([[2.0, 1.0], [1.0, -1.0]])# The right side gives the required value for each equation.system_targets = np.array([5.0, 1.0])# solve returns the vector whose matrix product matches the right side.system_solution = np.linalg.solve(system_matrix, system_targets)assert np.allclose(system_solution, [2.0, 1.0])# Step 6: Fit an approximate target vector using the earlier design matrix.# Deliberately add small inconsistencies so no exact linear fit is available.observed_targets = np.array([5.2, 1.9, 4.1, 5.8])# lstsq returns coefficients, residual summaries, numerical rank, and singular values.# The explicit relative cutoff controls which very small singular values count as zero.estimated_coefficients, residual_sums, numerical_rank, singular_values = np.linalg.lstsq( design_matrix, observed_targets, rcond=1e-12)# Reconstruct the fitted target vector in the design matrix's column space.fitted_targets = design_matrix @ estimated_coefficients# Residuals retain the same target units as the observed data.fitting_residuals = observed_targets - fitted_targets# This product checks orthogonality against every design column at once.orthogonality_check = design_matrix.T @ fitting_residualsprint("Estimated coefficients:", estimated_coefficients)print("Numerical rank:", numerical_rank)print("Column-residual dot products:", orthogonality_check)assert numerical_rank == design_matrix.shape[1]assert np.allclose(orthogonality_check, 0.0, atol=1e-10)
The residual check tests the defining property of the least-squares projection; it is more informative than printing a coefficient vector alone. A tiny floating-point residual dot product is expected. It does not imply that the observations have zero error or that the coefficients have a causal meaning.
6.7 An experiment with redundant features
If one column is duplicated, the model gains a parameter but no new output direction. The next experiment constructs this situation explicitly. We compare two coefficient vectors with the same predictions, then inspect the numerical rank. No random seed is needed because the example is deterministic.
# Step 7: Create an exactly redundant two-column representation.# The base feature has three observations and no intercept for this illustration.base_feature = np.array([1.0, 2.0, 3.0])# The second column is exactly twice the first, so the columns are dependent.redundant_matrix = np.column_stack((base_feature, 2.0* base_feature))# Define two different parameter choices whose combined multipliers both equal 3.# The first places all weight on the first feature.first_coefficients = np.array([3.0, 0.0])# The second places all weight on the doubled feature.second_coefficients = np.array([0.0, 1.5])assert np.allclose( redundant_matrix @ first_coefficients, redundant_matrix @ second_coefficients)# Step 8: Inspect numerical rank and the zero direction directly.# This vector changes coefficients without changing the resulting predictions.null_direction = np.array([-2.0, 1.0])assert np.allclose(redundant_matrix @ null_direction, 0.0)print("Redundant matrix rank:", np.linalg.matrix_rank(redundant_matrix))print("Identical predictions:", redundant_matrix @ first_coefficients)
The rank is one even though there are two columns. Adding any multiple of \((-2,1)\) to a coefficient vector leaves predictions unchanged. This set of directions mapping to zero is the null space. In regression, strong redundancy can make individual coefficients difficult to interpret even when predictions within the observed feature relationships appear reasonable.
6.8 Computation and common mistakes
Multiplying a dense \(n\times p\) matrix by a vector uses order \(np\) multiply- and-add work and stores order \(np\) entries for the matrix. These statements describe how work scales, not a measured runtime. Dense least-squares solving has additional factorization costs; for \(n\geq p\), common dense methods require order \(np^2\) work, with details depending on the solver and matrix structure. Sparse data and repeated solves can change practical choices.
The formula \((D^\top D)^{-1}D^\top\mathbf y\) is useful for analysis when the inverse exists, but explicitly constructing an inverse adds unnecessary work. Moreover, forming \(D^\top D\) squares the spectral condition number for a full-column-rank \(D\), potentially worsening numerical sensitivity. Use lstsq for the general unregularized least-squares calculation and inspect rank or singular values when identification matters.
Frequent mistakes include using * where @ is intended, transposing a one-dimensional array and expecting a new axis, reversing feature order, forgetting the intercept column, and assuming full rank because there are many rows. Another mistake is treating an exactly zero mathematical quantity as if it must equal the floating-point number zero after numerical factorization. Use a scale-appropriate tolerance and understand what it is checking.
6.9 Exercises and worked answers
Check dimensions. If \(X\) is \(20\times3\) and \(\mathbf w\) has three entries, what are the dimensions of \(X\mathbf w\) and \(X^\top X\)? Answer: 20 entries and \(3\times3\), respectively. The first produces one weighted sum per row; the second pairs every feature column with every feature column.
Calculate a norm. Find the Euclidean norm of \((-3,4)\). Answer:\(\sqrt{(-3)^2+4^2}=5\). The signs vanish only because of squaring, not because negative entries should be discarded.
Project a vector. Project \((3,1)\) onto the direction \((1,0)\). Answer: the multiplier is \(3/1=3\), so the projection is \((3,0)\). The residual \((0,1)\) has zero dot product with \((1,0)\).
Detect non-uniqueness. Two identical feature columns have coefficients 2 and 5. Give another pair with the same predictions. Answer: 3 and 4, because only their sum matters. Without additional information the data cannot identify their separate contributions.
Choose a solver. Ten noisy observations constrain two coefficients. Should solve or lstsq be the first choice? Answer:lstsq, because the matrix is rectangular and the equations will generally be inconsistent. A good approximate solution does not turn the system into an exact one.
Interpret orthogonality. Does \(D^\top\mathbf r=0\) imply every residual is zero? Answer: no. The residual can occupy a direction orthogonal to all columns. It vanishes only if the observed target lies in the column space.
Find a broadcasting error. Why is (n, 1) - (n,) dangerous for residual calculation? Answer: NumPy can produce an (n, n) array by broadcasting. Make both target and prediction shapes agree before subtraction.
6.10 One-page revision summary
Object or operation
Dimensions and meaning
\(X\in\mathbb R^{n\times p}\)
Observations in rows, features in columns.
\(\mathbf a^\top\mathbf b\)
Scalar sum of corresponding products.
\(X\mathbf w\)
\(n\) weighted sums; inner dimensions must match.
\(\|\mathbf r\|_2^2\)
Sum of squared residuals.
Column space
All output vectors the matrix can generate.
Rank
Number of independent column directions.
Projection
Closest reachable vector; residual is orthogonal to the space.
Normal equations; uniqueness needs full column rank.
solve / lstsq
Exact square nonsingular solving / general least squares.
For a new calculation, write dimensions beside each object before multiplying. For a fitted model, check predictions and residual geometry as well as coefficients. Chapter 7 next explains how uncertainty enters a model whose linear algebra may otherwise be entirely deterministic.