Linear Systems and Least Squares for AI Engineers

A model often has more observations than parameters and noisy data that cannot satisfy every equation exactly. Linear systems tell us when an exact answer exists. Least squares gives us a principled best approximation when it does not.

This is the mathematical core of linear regression, calibration, projection, feature fitting, and many optimization subproblems. It also teaches an engineering habit that matters far beyond linear algebra: inspect whether a problem is identifiable and numerically stable before trusting an answer.

By the end, you will be able to:

• read Ax=b as a collection of equations, a transformation, and a geometric intersection;
• distinguish unique, nonexistent, and non-unique exact solutions using rank;
• derive and interpret the least-squares objective and residual;
• build a regression design matrix and calculate a fitted line by hand;
• use solve and lstsq appropriately without explicitly computing an inverse;
• diagnose rank deficiency and ill-conditioning from residuals, singular values, and condition numbers.

Prerequisites: matrix multiplication, transpose, rank, vector norms, and basic NumPy arrays from the previous lesson.

One notation, three useful views#

Consider the system:

Write it as Ax=b:

The shapes are:

A: [m equations, n unknowns] = [2, 2]
x: [n unknowns]              = [2]
b: [m observations]          = [2]
A @ x                        = [2]

The same expression supports three mental models:

Eliminate x by subtracting half of the first equation from the second, or substitute y=x-2 into the first:

Verification is part of the solution, not an optional flourish:

import numpy as np

A = np.array([[2.0, 1.0], [1.0, -1.0]])
b = np.array([7.0, 2.0])
x = np.linalg.solve(A, b)

assert np.allclose(A @ x, b)
print(x)  # [3. 1.]

np.linalg.solve is for a square, full-rank coefficient matrix. It solves the required system directly. Computing np.linalg.inv(A) @ b does extra work and is generally a worse numerical interface.

Rank determines what the equations can identify#

The number of rows and columns alone does not tell us whether a system has one answer. Rank counts independent directions or independent constraints.

Unique solution

When a square n by n matrix has rank n, its columns span the entire target space and no column is redundant. Every b has exactly one solution.

No exact solution

Suppose two equations have identical left sides but different targets:

No pair (x,y) satisfies both. In column language, b is outside the column space of A. Least squares can still find the closest reachable prediction.

Infinitely many solutions

If the second equation merely repeats the first, such as 2x+2y=4, one constraint cannot identify two unknowns. Every point on x+y=2 is an exact solution.

Compare the rank of A with the augmented matrix formed by appending b:

This is not just classroom classification. Duplicate or derived features can make model parameters non-identifiable. A training pipeline may return coefficients, but different coefficient vectors can make the same predictions.

Overdetermined and underdetermined systems#

Let A have shape [m, n].

• Overdetermined: m>n. There are more equations or observations than unknowns. With noise, an exact solution usually does not exist.
• Square: m=n. A unique solution still requires full rank.
• Underdetermined: m<n. There are fewer independent constraints than unknowns, so multiple exact or equally good solutions may exist.

In machine learning, a design matrix commonly has one row per training example and one column per feature. A dataset with 100,000 examples and 20 features is overdetermined. A small dataset with millions of learned parameters is underdetermined, which is common in modern deep learning. The optimization objective and inductive biases then determine which of many fitting solutions training reaches.

NumPy's lstsq handles over-, well-, and underdetermined systems. If several solutions minimize the error, it returns the one with the smallest Euclidean norm.

The pseudoinverse unifies the non-square cases

For a square full-rank matrix, the inverse gives x=A^{-1}b. Non-square or rank-deficient matrices do not have that ordinary inverse. The Moore–Penrose pseudoinverse, written A^+, provides a useful generalization:

Its meaning depends on the system:

• for an inconsistent overdetermined system, A^+b is a least-squares

solution minimizing ||Ax-b||_2;

• when multiple least-squares solutions exist, it selects the one with minimum

Euclidean norm ||x||_2;

• for a consistent underdetermined system, it selects the minimum-norm exact

solution.

This is a mathematical contract, not a recommendation to form a pseudoinverse explicitly in routine code. numpy.linalg.lstsq computes a solution through a factorization and returns rank and singular-value diagnostics. Use numpy.linalg.pinv when the pseudoinverse matrix itself is genuinely required, not merely to solve one system.

Least squares chooses the closest reachable target#

For a candidate parameter vector x, define the prediction and residual:

Ordinary least squares selects:

Every symbol has a job:

• x_hat is the fitted parameter vector;
• Ax is the target predicted by those parameters;
• r is what the model failed to explain;
• ||r||_2^2 adds the squared residual components.

Squaring prevents positive and negative errors from canceling and creates a differentiable objective. It also makes large errors influential, which is useful under some noise assumptions but sensitive to outliers.

Geometrically, all possible predictions Ax lie in the column space of A. Least squares projects b onto that space. At the optimum, the remaining residual is perpendicular to every column of A:

Rearranging gives the normal equations:

This derivation explains the answer, but it is not usually the best numerical algorithm. Forming A^TA can amplify conditioning problems. Production libraries instead use factorizations such as QR or SVD-backed least-squares drivers.

Text fallback: parameters pass through the design matrix to produce the closest reachable prediction. Subtracting that prediction from the observed target produces a residual orthogonal to the design matrix's columns.

Fit a line completely by hand#

Fit y=c+mx to three observations:

The parameter vector is theta=[c,m]: intercept first, slope second. The constant column of ones lets matrix multiplication represent the intercept:

The system has three observations and two parameters. No line passes through all three points, so solve the normal equations:

Therefore:

Solving gives c=7/6 and m=1/2. The fitted values are:

Using r=b-y_hat, the residual is:

Check the geometry:

The residuals sum to zero because they are orthogonal to the intercept column of ones. Their feature-weighted sum is zero because they are also orthogonal to the x column.

x_values = np.array([0.0, 1.0, 2.0])
targets = np.array([1.0, 2.0, 2.0])
A = np.column_stack((np.ones_like(x_values), x_values))

theta, squared_residuals, rank, singular_values = np.linalg.lstsq(
    A, targets, rcond=None
)
residual = targets - A @ theta

print(theta)              # [1.16666667 0.5]
print(A.T @ residual)     # approximately [0. 0.]

The outputs are diagnostics, not clutter:

• theta contains one coefficient per design column;
• squared_residuals reports the residual sum of squares only in cases documented by NumPy;
• rank reports the effective design rank;
• singular_values reveal how strongly the matrix represents independent directions.

Linear regression is a least-squares design problem#

For N examples with D features:

X:       [N, D]
weights: [D]
target:  [N]
X @ weights → [N]

To fit an intercept, either add a column of ones or let a regression library manage centering/intercept behavior. Do not do both accidentally.

The word linear refers to being linear in the fitted parameters. A polynomial feature such as x^2 can appear as another design column while the model remains linear in its coefficients:

design row = [1, x, x²]
parameters = [β₀, β₁, β₂]

This separates two concerns: feature construction defines what patterns are expressible; least squares selects coefficients for those features.

[age=2, usage=10]
[age=5, usage=4]

Conditioning: a valid answer can still be unreliable#

Rank asks whether a direction is effectively missing. Conditioning asks how sensitive the answer is before a direction disappears entirely.

The 2-norm condition number is the ratio of the largest singular value to the smallest:

A condition number near one means comparable sensitivity across directions. A very large value means small input perturbations or floating-point errors can create large coefficient changes. Exactly singular matrices have an infinite condition number in exact arithmetic.

condition = np.linalg.cond(A)

Common causes include:

• features with radically different numerical scales;
• duplicate or nearly duplicate columns;
• polynomial features over a large raw range;
• insufficient data variation;
• forming A^TA, whose condition number is roughly squared in the 2-norm when A has full column rank.

Conditioning changes how you interpret results. A fit can have a small training residual while individual coefficients are unstable. Predictions near the observed data may look reasonable, yet extrapolation or coefficient explanations can be untrustworthy.

Useful responses include rescaling or centering features, removing redundant variables, collecting more informative data, using an SVD/QR-based solver, or applying regularization when the modeling goal supports it. Regularization changes the objective; it is not merely a numerical switch.

Choose the solver that matches the problem#

Avoid these habits:

• Computing an inverse to solve: use solve or lstsq.
• Assuming square means invertible: check rank/solver behavior.
• Using normal equations blindly: they are valuable for derivation but can worsen conditioning.
• Reading coefficients without feature units: a coefficient changes with feature scale.
• Treating low residual as proof of generalization: least squares optimizes the supplied observations, not unseen data.
• Ignoring outliers: squared loss gives large residuals disproportionate influence.
• Comparing train and production features built differently: design-column order, centering, encoding, and units are part of the model contract.

A repeatable debugging checklist#

When a linear solve or regression looks wrong:

1. Print A.shape, b.shape, and the meaning and unit of every axis.
2. Reject NaN and infinite inputs before invoking a solver.
3. Check effective rank and singular values.
4. Inspect the condition number in the same dtype used by production.
5. Reconstruct predictions with A @ x.
6. Inspect residuals, not only their total squared norm.
7. For ordinary least squares, verify A.T @ residual is near zero within a scale-aware tolerance.
8. Perturb inputs slightly and see whether coefficients or predictions move dramatically.
9. Test on held-out data before making a generalization claim.

For float comparisons, use tolerances such as np.allclose; exact equality is rarely the right invariant after numerical linear algebra. Choose tolerances relative to the scale and precision of the data rather than copying one universal threshold.

Memory chart#

Ax = b
A [m,n] maps parameters x [n] to observations [m]

square + full rank       → solve(A,b)
noisy or rectangular     → lstsq(A,b)
underdetermined          → multiple candidates; lstsq returns minimum norm

prediction = Ax
residual   = b - Ax
least squares minimizes ||residual||²
at optimum: Aᵀ residual ≈ 0

small residual ≠ stable coefficients
check rank + singular values + condition number
avoid inverse; avoid normal equations as the default implementation

The next lesson studies norms, distances, and similarity. The residual norm used here will become one member of a larger family of ways to measure vector size and difference—and you will learn how that choice changes retrieval, clustering, losses, and model behavior.

References#

• Deisenroth, Faisal, and Ong: Mathematics for Machine Learning, Sections 3.8 and 9.2 — orthogonal projections, least squares, and linear-regression parameter estimation.
• Goodfellow, Bengio, and Courville: Deep Learning, Sections 2.3 and 2.9 — linear systems, inverses, and the Moore–Penrose pseudoinverse.
• NumPy: linalg.solve — exact solutions for square, full-rank systems.
• NumPy: linalg.pinv — SVD-based Moore–Penrose pseudoinverse behavior and cutoff policy.
• NumPy: linalg.lstsq — least-squares behavior for under-, well-, and overdetermined systems, including returned diagnostics.
• NumPy: linalg.cond — matrix condition numbers under different norms.
• SciPy: linalg.lstsq — solver drivers, effective-rank cutoff, batched behavior, and residual outputs.
• scikit-learn: LinearRegression — ordinary least-squares regression behavior and estimator shape conventions.