Orthogonality, projections, and least squares

An equation Xw=y can fail to have a solution while still admitting one uniquely closest predicted vector. The distinction matters whenever measured targets do not lie in the design’s column space. Least squares chooses that closest prediction according to a specified distance. It does not by itself make the coefficients unique, the measurements trustworthy, or the predictions reliable on future records.

This module derives the geometry before introducing algorithms. You will prove projection formulas, construct an orthonormal basis, solve with QR, and establish the normal equations without differentiation. The last section describes the whole family of rank-deficient minimisers and its minimum-norm representative. The labs compare mathematical equivalence with floating-point computation and data sensitivity.

Retrieval check: give the column space and null space of a matrix and explain why a null displacement leaves predictions unchanged. Review Module 11 if needed. Dot products and lengths come from Module 09. No calculus is assumed; a later lesson will derive the same normal equations by differentiation. Use the pinned environment in NumPy preparation.

A real inner product assigns a real number ⟨u,v⟩ to two vectors. It is linear in each argument, symmetric, and positive definite: ⟨v,v⟩≥0 with equality only for v=0. These conditions generalise the coordinate dot product while preserving useful geometric reasoning. The associated norm is √⟨v,v⟩. Orthogonality means ⟨u,v⟩=0. It refers to the selected inner product, so changing coordinate weights can change which vectors are orthogonal even though addition and linear dependence remain the same.

The Euclidean product is uᵀv=Σuᵢvᵢ. A weighted product is uᵀWv for a symmetric positive-definite W. For a diagonal W with entries αᵢ>0, this is Σαᵢuᵢvᵢ. Symmetry and bilinearity follow term by term, and a nonzero vector has at least one nonzero coordinate contributing a strictly positive square. Thus this diagonal construction satisfies every condition. It corresponds to scaling coordinates by √αᵢ before using the ordinary product. A zero weight instead allows a nonzero vector to have zero length, so the expression is only semidefinite and not an inner product on the full coordinate space. Negative weights fail even nonnegativity.

If u and v are orthogonal, expansion gives ‖u+v‖²=‖u‖²+2⟨u,v⟩+‖v‖²=‖u‖²+‖v‖². This Pythagorean identity is the engine of projection and least squares. Without orthogonality the cross term remains, and adding squared lengths is generally wrong. The vectors need not be coordinate axes; any pair perpendicular under the specified product works. A decomposition into mutually orthogonal components extends the identity to any finite sum by the same expansion.

An orthonormal list q₁,…,q_r has unit vectors that are pairwise orthogonal. Such a list is automatically independent: dot a zero combination with q_j and its jth coefficient must vanish. For a matrix Q with these columns, QᵀQ=I_r. This identity has dimensions r×r. It does not generally imply QQᵀ=I_m when Q is m×r with r<m; the latter operation will be the projection onto the smaller column space. Calling rectangular Q an orthogonal square matrix would confuse two different contracts.

In data fitting, a weighted distance encodes a modelling choice about which discrepancies count more. Large weights penalise some coordinate errors more strongly. This can reflect known measurement precision or an application-specific objective, but the weights need justification. Simply declaring a weighted inner product does not supply a statistical noise model. For the main derivations below, use Euclidean geometry unless explicitly labelled otherwise. This keeps the objective, residual-orthogonality condition and algorithm consistent.

The zero vector is perpendicular to every vector but has no unit normalisation or defined angle. A procedure asking for a direction must guard against dividing by its norm. Similarly, a pair of nonzero vectors is orthogonal exactly when their angle has zero cosine in the Euclidean geometry, but using cosine as a probability would add an interpretation unsupported by these algebraic facts. Here orthogonality is an exact relation that makes a decomposition provable, rather than a vague claim of feature “unrelatedness.”

Given a nonzero vector u and a target y, seek its closest point p=cu on the line span{u}. Requiring the residual r=y−cu to be perpendicular to u gives uᵀy−c uᵀu=0. Since uᵀu>0, the unique coefficient is c=(uᵀy)/(uᵀu). This coefficient is a scalar; p is a vector in the original target space. They coincide numerically only in particular coordinate conventions, so an implementation should name them separately.

Why does the perpendicular choice minimise distance? For any other scalar a, y−au=r+(c−a)u. The two terms are orthogonal, so ‖y−au‖²=‖r‖²+(c−a)²‖u‖². The second term is nonnegative and vanishes exactly at a=c. This is a complete optimality and uniqueness proof without calculus. If the generator is zero, its span is the zero space and the projection is simply zero; the fraction is undefined and unnecessary. Rescaling a nonzero generator changes c inversely and leaves p unchanged, as a subspace projection should.

For a finite-dimensional subspace S with orthonormal basis columns Q, define p=QQᵀy. Its coefficients in that basis are Qᵀy. The residual r=y−p satisfies Qᵀr=Qᵀy−QᵀQ Qᵀy=0, so it is perpendicular to every basis vector and hence every vector of S. For any other s∈S, y−s=r+(p−s), an orthogonal sum. Pythagoras gives ‖y−s‖²=‖r‖²+‖p−s‖², proving existence of a closest point and uniqueness of the predicted vector p. Section 3 supplies the orthonormal-basis construction needed for any finite-dimensional S.

The projection matrix P=QQᵀ is symmetric and idempotent: Pᵀ=P and P²=Q(QᵀQ)Qᵀ=P. Idempotence means projecting an already projected vector changes nothing. Symmetry reflects Euclidean perpendicularity. The image of P is S, since Py∈S and Ps=s for s∈S. Its kernel is S⊥ because Qᵀy=0 exactly when y is perpendicular to the basis. The complementary projection I−P sends y to its residual. These operations express one unique orthogonal decomposition y=p+r, not an invertible change of coordinates.

Pythagoras also proves that projection cannot increase Euclidean length: ‖Py‖²≤‖y‖². Applying this to a difference gives ‖Py−Pz‖≤‖y−z‖, since P is linear. Thus fitted observation vectors vary no more than their targets for a fixed subspace, even though recovering coefficient coordinates from them can be very sensitive. This gives a precise reason to discuss stability of predictions separately from stability of parameters.

If U has independent but nonorthonormal columns spanning S, perpendicularity gives Uᵀ(y−Uc)=0. The matrix UᵀU is invertible: cᵀUᵀUc=‖Uc‖² is positive for every nonzero c by independence. Thus c solves (UᵀU)c=Uᵀy, and formally P=U(UᵀU)⁻¹Uᵀ. The formula explains the map; it does not prescribe explicitly calculating an inverse as the best numerical algorithm. QR will compute the same projection without forming that Gram matrix.

For diagonal positive weights, the line coefficient becomes (uᵀWy)/(uᵀWu), and the residual satisfies uᵀWr=0. With y=(3,2),u=(1,1),W=diag(4,1), it is 14/5, giving p=(2.8,2.8). The residual (0.2,−0.8) is weighted-orthogonal but not Euclidean-orthogonal to u. The closest point changed because the objective changed. To project onto an affine set p₀+S, project y−p₀ onto S then add p₀; an offset feasible set should not be silently treated as a subspace through zero.

Start with independent vectors a₁,…,a_d in Rᵐ. Gram–Schmidt constructs orthonormal vectors with the same successive spans. First take q₁=a₁/‖a₁‖. At step j, remove the components along earlier q’s: v_j=a_j−Σ_{i<j}(qᵢᵀa_j)qᵢ, then set q_j=v_j/‖v_j‖. Dotting v_j with any earlier q_k gives q_kᵀa_j−q_kᵀa_j=0, because the orthonormal cross products remove every other term. Thus the remainder is perpendicular to the existing span.

The remainder cannot be zero in exact arithmetic for an independent input list. If it were, a_j would already be a combination of earlier q’s and therefore of earlier a’s, contradicting independence. Normalisation consequently creates a well-defined unit vector. Each a_j is a combination of q₁,…,q_j, and each q_j is a combination of a₁,…,a_j, so the successive spans are equal. Induction establishes both the construction and its promised span preservation. This proves every finite-dimensional Euclidean subspace admits an orthonormal basis, starting from any ordinary basis.

Writing all columns together gives A=QR. For m≥d and independent input columns, reduced Q is m×d with QᵀQ=I_d, while R is d×d upper triangular. Its entries above the diagonal are rᵢⱼ=qᵢᵀa_j and its positive diagonal entries are the remainder lengths. Nonzero diagonal makes R invertible. This triangular structure records the order in which spans were built. Changing the sign of one q column and the matching row of R leaves QR unchanged. A QR returned by a library can therefore differ in signs from the hand construction without being incorrect.

For dependent columns, the exact zero remainder is a diagnostic. Do not divide by it or claim that d orthonormal directions were obtained. One can skip redundant vectors to construct an r-column basis for the rank-r column space, or use a factorisation designed to handle dependence. In floating-point arithmetic, a very small remainder needs a stated resolution policy; it may reflect exact dependence, near dependence or accumulated error. A fixed absolute threshold changes meaning when column units change. The simple lab guard illustrates the decision but does not provide a general rank-revealing solver.

Classical Gram–Schmidt computes all removal coefficients against the original a_j and subtracts their combined contribution. Modified Gram–Schmidt removes one component at a time and computes the next coefficient against the current remainder. They give the same exact formulas, because earlier q’s are exactly orthogonal. Floating-point rounding can make the procedures behave differently; sequential removal often preserves orthogonality better, but it is not immune to severe near dependence. Reorthogonalisation or Householder QR can address algorithmic error. Neither creates information missing from an ill-conditioned design.

An orthonormal basis makes subspace coordinates particularly simple: qᵢᵀy reads the component in direction i, and the sum of squared components equals the squared norm of the projection. The residual norm accounts for the missing perpendicular part. For a square orthogonal Q, these coordinates preserve the entire vector norm because there is no missing part. For a rectangular Q, Qᵀy retains only its coordinates in the smaller column space. This is the geometric difference between a full coordinate change and dimension reduction.

For a full-column-rank design X of shape m×d with m≥d, use reduced X=QR. Decompose y as Q Qᵀy+r with Qᵀr=0. Then y−Xw=r+Q(Qᵀy−Rw). Pythagoras and Q’s length preservation on d-dimensional coordinates give ‖y−Xw‖²=‖r‖²+‖Qᵀy−Rw‖². The first term is fixed. Because R is invertible, the second is zero exactly when Rw=Qᵀy. Back substitution solves this d-dimensional triangular system and gives a unique w.

The shapes are explicit: y is length m, Qᵀy is length d, R is d×d, and w is length d. The predicted vector Q Qᵀy is length m. A library implementation may store y as a column or allow several target columns, but those storage conventions must be declared separately. For multiple targets Y of shape m×k, the same R can solve for d×k weights with RHS QᵀY; each column is its own target problem. Factorisation reuse avoids recomputing geometry unnecessarily.

Householder QR offers another orthonormal construction. For a unit vector v, H=I−2vvᵀ is symmetric and HᵀH=I, since expansion cancels −4vvᵀ with 4v(vᵀv)vᵀ. It is a norm-preserving reflection. For a nonzero column segment x, take v proportional to x+s‖x‖e₁ where s matches the sign of the first coordinate, choosing s=1 at zero. This avoids subtracting nearly equal leading quantities. The reflection sends x to −s‖x‖e₁. Applying successive reflections to trailing column segments zeros the entries below successive diagonals; their orthogonal product supplies Q. This explains why a practical QR need not use the simple educational Gram–Schmidt routine.

Avoid forming a Gram matrix when its loss of resolution matters. For X with columns (1,1,1) and (1,1+δ,1−δ), exact XᵀX=[[3,3],[3,3+2δ²]]. The independent difference δ becomes a much smaller squared term in the Gram matrix. With δ=10⁻⁸ and local double-precision arithmetic, 3+2δ² rounds to three in the lab, producing a singular computed matrix. The original floating-point X still retains its distinct column entries, so QR and a least-squares library can use that information. The exact formulas are equivalent; their evaluated numerical intermediates are not.

This example does not establish that every normal-equation computation fails. A well-conditioned small problem may work accurately, as the hand line-fit calculation does. It demonstrates an avoidable source of damage for near-collinear features. The familiar statement that a Gram matrix squares the two-norm condition number will be derived from singular values in Module 14. Here the δ² term already shows the mechanism explicitly without assuming that theorem.

Algorithmic stability and problem sensitivity remain separate. A sound QR calculation can return accurate coefficients for the supplied numbers yet produce very different coefficients after a tiny target perturbation along a weak direction. Lab 3 adds (0,δ,−δ) to y, changing weights by approximately (−1,1) while retaining a tiny residual. A small fitted discrepancy establishes agreement with the current targets, not stable recovery of physically true coefficients. Appropriate feature scaling, measurement assumptions and later regularisation can alter the modelling or resolution decisions; they do not change this distinction.

For arbitrary X, minimise the Euclidean objective F(w)=‖y−Xw‖². This is the least-squares criterion; the word squares refers to the sum of coordinate residual squares. Dividing by the fixed number of observations or multiplying by a positive constant changes the objective value but not its minimisers. Choosing another metric, adding a penalty or restricting w can change the minimising set, so each altered problem needs its own derivation.

The column space S=C(X) is finite dimensional. Its orthogonal projection p of y exists and is unique by Section 2. Because p belongs to S, at least one w satisfies Xw=p. These and only these weights minimise F: for any candidate, ‖y−Xw‖²=‖y−p‖²+‖p−Xw‖². Thus existence of minimising weights holds even for a rank-deficient design. The fitted vector is uniquely determined by the subspace and the chosen distance, while the map from weights to that vector may be many-to-one.

Residual perpendicularity gives the normal equations Xᵀ(y−Xw)=0, equivalently XᵀXw=Xᵀy. Necessity can also be proved by a direct scalar variation. If residual r has a nonzero dot product with a column a_j, change only weight j by t. The objective change is t²‖a_j‖²−2t a_jᵀr. The column is nonzero in this situation. Choosing t=(a_jᵀr)/‖a_j‖² gives a strictly negative change, contradicting minimality. Consequently every column dot product is zero. This argument has no derivative and covers all unconstrained real weights.

For sufficiency, suppose Xᵀr=0. A weight change h gives residual r−Xh, with rᵀXh=0. Expansion yields F(w+h)=F(w)+‖Xh‖²≥F(w). Hence every normal-equation solution is a global minimiser. Equality holds exactly for h∈N(X). We have proved both directions and the entire minimiser family, rather than merely finding a stationary candidate. With full column rank, Xh=0 forces h=0 and the minimiser is unique. Otherwise all weights in w+N(X) minimise with the same fitted vector.

With an intercept column of ones, residual perpendicularity implies their sum is zero. Without that column, there is no such automatic guarantee. An exact-fit target has residual zero; zero is perpendicular to every column, so the normal equations include interpolation as a special case. The residual can be nonzero while Xᵀr is zero, as in the line-fit example. Confusing the two arrays would make a correct approximate fit appear exact.

Weighted least squares with symmetric positive-definite W minimises rᵀWr. The same perpendicular-decomposition or variation argument gives XᵀWr=0 and XᵀWXw=XᵀWy. For diagonal weights, one can scale rows of X and y by the square roots of the weights and solve an ordinary Euclidean problem. Zero or negative weights require separate treatment, because the inner-product proof no longer applies unchanged. The numerical shape and interpretation of weights are part of the specification.

Normal equations describe optimality; they do not require using their Gram matrix as the algorithm. A QR or SVD solution can satisfy them to numerical tolerance while avoiding that intermediate. When testing a computed fit, inspect the residual and Xᵀr with a scale-aware tolerance and the declared rank. A good optimality check establishes that the specified objective was solved well. It still supplies no independent evidence about noise, omitted variables, causal meaning or extrapolation.

For a rank-deficient design, a least-squares minimiser w₀ determines the family w₀+N(X). To choose one representative, a common additional rule minimises the Euclidean parameter norm. Module 11 gave the orthogonal decomposition of parameter space into R(X) and N(X). Write w₀=u+k with u in the row space and k in the null space. Since Xk=0, Xu=Xw₀, so u is also a least-squares minimiser. Every equivalent minimiser has form u+n with n∈N(X). Orthogonality gives ‖u+n‖²=‖u‖²+‖n‖², uniquely minimised at n=0.

Thus the minimum-norm least-squares solution is the unique minimiser that belongs to the row space. Its uniqueness comes from an added choice of parameter metric, not from information newly supplied by the observations. Reparameterising or rescaling features can change which coefficient vector has smallest ordinary Euclidean norm, even when the possible predictions are unchanged. The coordinate convention therefore matters for this selection rule.

The map taking y to this unique representative is the pseudoinverse X⁺. Its geometry has two stages: project y onto C(X), then apply the inverse of X restricted from R(X) to C(X). That restricted map is bijective by Module 11: its kernel in the row space is zero, and it reaches the whole column space. Projection and this inverse are linear, so their composition is a linear map from Rᵐ to Rᵈ, represented by a d×m matrix. This construction supplies the pseudoinverse without requiring a square inverse of X.

It satisfies XX⁺=P_C and X⁺X=P_R, the orthogonal projections onto the column and row spaces. The first identity follows because projecting and lifting then remapping returns the projected target. The second follows because any parameter splits into row-space and null-space parts, and the lift selects exactly its row-space part. Therefore XX⁺X=X and X⁺XX⁺=X⁺, and both XX⁺ and X⁺X are symmetric. These four identities are the Moore–Penrose conditions. The geometric construction proves them and explains why the rectangular pseudoinverse is different from an ordinary two-sided inverse. Module 14 will derive its singular-value formula and numerical cutoff issues.

A numerical pseudoinverse may discard sufficiently weak directions according to a tolerance. That changes the effective subspaces used in the calculation. A returned minimum-norm vector concerns that numerical policy and the supplied array; keep the exact mathematical construction separate. In the pinned least-squares interface, rank-deficient designs can also yield an empty returned residual-summary array even when the actual residual is nonzero. Compute y−Xw explicitly to assess the discrepancy rather than treating an empty summary as a proof of exact fit.

Least squares finally separates three questions. What vector is closest to the observed targets under the objective? Which coefficient representative is selected among equivalent fits? What prediction or interpretation should transfer to new records? Orthogonal projection answers the first, and full rank or an added convention answers the second. The third needs a feature-domain assumption and later statistical evidence. A zero training residual can arise from interpolation and still accompany poor future behaviour; a small residual can coexist with unstable coefficients. Keep these limits with the useful algebra, rather than attaching unsupported conclusions to a successful solver call.

Use the Python 3.11 and NumPy baseline in array preparation. Outputs below are captured by running the downloadable scripts. Predict the shapes and expected identities first; use tolerances for floating-point equalities and distinguish them from the exact proofs.

Lab A Projections and orthonormalisation

Compare the Euclidean and weighted diagonal-line projections. Predict which residual dot product should be zero in each geometry. Trace modified Gram–Schmidt, verify QᵀQ and QR, and explain the duplicate-column guard rather than removing it to continue execution.

Download lab1_projections_and_bases.py

"""Projection depends on the inner product; orthonormalisation preserves span."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)

y = np.array([3., 2.])
u = np.array([1., 1.])
coefficient = (u@y)/(u@u)
projection = coefficient*u
residual = y-projection
assert np.allclose(projection, [2.5, 2.5]) and np.isclose(u@residual, 0)
print("Euclidean coefficient / projection / residual:", coefficient, projection, residual)
weights = np.array([4., 1.])
weighted_coefficient = (u@(weights*y))/(u@(weights*u))
weighted_residual = y-weighted_coefficient*u
print("Weighted coefficient / residual:", weighted_coefficient, weighted_residual)
print("Weighted / Euclidean residual dot:", u@(weights*weighted_residual), u@weighted_residual)
assert np.isclose(weighted_coefficient, 2.8) and np.isclose(u@(weights*weighted_residual), 0)

def modified_gram_schmidt(A, tolerance=1e-12):
    Q = np.zeros_like(A, dtype=float)
    R = np.zeros((A.shape[1], A.shape[1]))
    for j in range(A.shape[1]):
        remaining = A[:, j].astype(float).copy()
        for i in range(j):
            R[i, j] = Q[:, i]@remaining
            remaining -= R[i, j]*Q[:, i]
        R[j, j] = np.linalg.norm(remaining)
        if R[j, j] <= tolerance:
            raise ValueError("Dependent or unresolved column: do not normalise zero")
        Q[:, j] = remaining/R[j, j]
    return Q, R

A = np.array([[1., 1.], [1., 0.], [0., 1.]])
Q, R = modified_gram_schmidt(A)
print("Q:", Q.tolist())
print("R:", R.tolist())
print("Orthogonality / reconstruction checks:", np.allclose(Q.T@Q, np.eye(2)), np.allclose(Q@R, A))
assert np.allclose(Q.T@Q, np.eye(2)) and np.allclose(Q@R, A)
try:
    modified_gram_schmidt(np.ones((3, 2)))
except ValueError as error:
    print("Duplicate-column guard:", error)
else:
    raise AssertionError("A zero remainder must not be normalised")

Euclidean coefficient / projection / residual: 2.5 [2.5 2.5] [ 0.5 -0.5]
Weighted coefficient / residual: 2.8 [ 0.2 -0.8]
Weighted / Euclidean residual dot: 8.881784197001252e-16 -0.5999999999999996
Q: [[0.7071067811865475, 0.40824829046386313], [0.7071067811865475, -0.40824829046386296], [0.0, 0.8164965809277261]]
R: [[1.4142135623730951, 0.7071067811865475], [0.0, 1.224744871391589]]
Orthogonality / reconstruction checks: True True
Duplicate-column guard: Dependent or unresolved column: do not normalise zero

Deliverable: the two objectives, projected vectors and residual conditions, followed by a two-column factorisation with shapes and reconstruction. Explain why tiny roundoff residual dots are compatible with an exact orthogonality theorem.

Lab B QR fitting and the minimum-norm family

Fit the three records with reduced QR and the library least-squares solver. Check residual perpendicularity explicitly. Duplicate the x feature, predict the weight family, and verify that the minimum-norm representative is perpendicular to the null direction. Inspect the returned residual summary separately from the computed residual.

Download lab2_qr_and_minimum_norm.py

"""A unique fitted vector may have an affine family of minimising weights."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)
x = np.array([0., 1., 2.])
X = np.column_stack((np.ones(3), x))
y = np.array([1., 2., 2.])
Q, R = np.linalg.qr(X, mode="reduced")
w_qr = np.linalg.solve(R, Q.T@y)
w_library, returned_residuals, rank, values = np.linalg.lstsq(X, y, rcond=None)
residual = y-X@w_qr
assert np.allclose(w_qr, [7/6, 1/2]) and np.allclose(w_qr, w_library)
assert np.allclose(X.T@residual, 0)
print("X / Q / R shapes:", X.shape, Q.shape, R.shape)
print("QR / library weights:", w_qr, w_library)
print("Fitted values / residual:", X@w_qr, residual)
print("Residual sum of squares:", residual@residual)
print("Normal-equation check:", np.allclose(X.T@residual, 0))

duplicate = np.column_stack((np.ones(3), x, x))
w_min, library_residuals, duplicate_rank, _ = np.linalg.lstsq(duplicate, y, rcond=None)
null = np.array([0., -1., 1.])
alternative = w_min + 2*null
assert np.allclose(w_min, [7/6, 1/4, 1/4])
assert np.allclose(duplicate@alternative, duplicate@w_min)
assert np.isclose(w_min@null, 0)
assert np.allclose(np.linalg.pinv(duplicate)@y, w_min)
print("Duplicate design rank / minimum-norm weights:", duplicate_rank, w_min)
print("Alternative weights / same fit:", alternative, np.allclose(duplicate@alternative, duplicate@w_min))
print("Squared norms:", w_min@w_min, alternative@alternative)
print("Rank-deficient returned residuals:", library_residuals.tolist())
print("Explicit rank-deficient SSE:", np.linalg.norm(y-duplicate@w_min)**2)
assert library_residuals.size == 0 and np.linalg.norm(y-duplicate@w_min) > 0
print("An empty library residual array does not assert a zero actual residual.")

X / Q / R shapes: (3, 2) (3, 2) (2, 2)
QR / library weights: [1.166667 0.5     ] [1.166667 0.5     ]
Fitted values / residual: [1.166667 1.666667 2.166667] [-0.166667  0.333333 -0.166667]
Residual sum of squares: 0.16666666666666674
Normal-equation check: True
Duplicate design rank / minimum-norm weights: 2 [1.166667 0.25     0.25    ]
Alternative weights / same fit: [ 1.166667 -1.75      2.25    ] True
Squared norms: 1.4861111111111098 9.48611111111111
Rank-deficient returned residuals: []
Explicit rank-deficient SSE: 0.1666666666666668
An empty library residual array does not assert a zero actual residual.

Deliverable: the unique predicted vector, the full equivalent coefficient family, its minimum-norm member and SSE. Explain why the library’s empty rank-deficient summary does not mean that all records were fitted exactly.

Lab C Gram rounding and target sensitivity

Write the exact Gram matrix before running. Compare it with the rounded matrix, then compare QR and library weights. Perturb the target along the weak direction and measure the coefficient change and residual. The captured normal-equation result is for this baseline; another arithmetic backend may produce inaccurate solvable coefficients instead of the same exception.

Download lab3_gram_rounding.py

"""Forming a Gram matrix can erase a weak direction; QR cannot cure noisy data."""
import numpy as np
np.set_printoptions(precision=9, suppress=True)
delta = 1e-8
X = np.array([[1., 1.], [1., 1.+delta], [1., 1.-delta]])
truth = np.array([1., 2.])
y = X@truth
gram = X.T@X
Q, R = np.linalg.qr(X, mode="reduced")
w_qr = np.linalg.solve(R, Q.T@y)
w_library, _, rank, _ = np.linalg.lstsq(X, y, rcond=None)
print("Rounded Gram matrix:", gram.tolist())
print("QR / library weights:", w_qr, w_library)
print("Library rank / condition number of X:", rank, format(np.linalg.cond(X), ".6e"))
assert rank == 2 and np.allclose(w_qr, truth, atol=1e-6)
try:
    w_normal = np.linalg.solve(gram, X.T@y)
except np.linalg.LinAlgError:
    print("Normal equations: singular after Gram rounding")
else:
    print("Normal-equation weights / coefficient error:", w_normal, np.linalg.norm(w_normal-truth))
perturbed = y + np.array([0., delta, -delta])
changed = np.linalg.solve(R, Q.T@perturbed)
print("Target perturbation norm:", format(np.linalg.norm(perturbed-y), ".6e"))
print("QR weights after target perturbation:", changed)
print("Weight change norm:", format(np.linalg.norm(changed-w_qr), ".6e"))
print("Perturbed residual norm:", format(np.linalg.norm(perturbed-X@changed), ".6e"))
assert np.linalg.norm(changed-w_qr) > 1
assert np.linalg.norm(perturbed-X@changed) < 1e-12
print("A stable algorithm avoids avoidable Gram damage; a tiny residual does not guarantee insensitive coefficients.")

Rounded Gram matrix: [[3.0, 3.0], [3.0, 3.0]]
QR / library weights: [0.99999999 2.00000001] [1.000000011 1.999999989]
Library rank / condition number of X: 2 2.449490e+08
Normal equations: singular after Gram rounding
Target perturbation norm: 1.414214e-08
QR weights after target perturbation: [0.000000001 2.999999999]
Weight change norm: 1.414214e+00
Perturbed residual norm: 4.440892e-16
A stable algorithm avoids avoidable Gram damage; a tiny residual does not guarantee insensitive coefficients.

Deliverable: separate the avoidable Gram-matrix rounding loss from the design’s intrinsic sensitivity to target changes. A stable solver avoids the first issue and cannot remove the second. Explain what a small residual actually certifies.

Exercises 1–12 are required; 13–14 add 35 optional minutes. Unless labelled weighted, every inner product is Euclidean and residual means y−Xw.

Use orthogonality, projection and linear-regression selections in the authors’ Mathematics for Machine Learning companion. For numerical contracts, read official NumPy 1.26 documentation for QR, least squares, and pseudoinverse. These supplement the original proofs and examples above.

Without notes, derive the line projection, prove subspace uniqueness, and construct a two-column QR. Derive both directions of the normal equations without calculus. Distinguish fit uniqueness, weight uniqueness and minimum-norm selection, then describe one numerical and one modelling limitation.

Exit task: X=[[1],[1],[1]],y=(1,2,3) has weight two, fit (2,2,2), residual (−1,0,1) and SSE two. The residual dot with X’s column vanishes. Duplicating that column yields weights (t,2−t), with minimum norm at (1,1). Explain why the fitted prediction did not change and why neither result supplies a future-data guarantee.

Ready to move on: you can identify the relevant spaces, objective and algorithm before interpreting a fit. Module 13 develops eigenvalues and quadratic forms; Module 14 will derive singular values and the computational pseudoinverse. Check the course overview for availability.