III-10 · Multiple Linear Regression

Read a prediction as one matrix row

Edit design cells and coefficients, inspect XW row by row, and solve full-rank or dependent designs with a minimum norm SVD fit.

Fit and optimization laboratory

Design rows → shared coefficients → predictions

Hand-authored toy observations, editable in this sheet. No seed and no measured population.

● Training observation · ◇ Prediction · dashed residual. Axis domain adjusts to the data; x/y units share one scale.

Change the active design and coefficient column

The intercept column is fixed at 1. All other design cells, targets and coefficients are editable.

Four sample rows of X and target column y
ID / bias	x₁	x₂	y
s1 / 1	X row 1 feature 1-1000000 to 1000000 · Enter to apply	X row 1 feature 2-1000000 to 1000000 · Enter to apply	Target row 1-1000000 to 1000000 · Enter to apply
s2 / 1	X row 2 feature 1-1000000 to 1000000 · Enter to apply	X row 2 feature 2-1000000 to 1000000 · Enter to apply	Target row 2-1000000 to 1000000 · Enter to apply
s3 / 1	X row 3 feature 1-1000000 to 1000000 · Enter to apply	X row 3 feature 2-1000000 to 1000000 · Enter to apply	Target row 3-1000000 to 1000000 · Enter to apply
s4 / 1	X row 4 feature 1-1000000 to 1000000 · Enter to apply	X row 4 feature 2-1000000 to 1000000 · Enter to apply	Target row 4-1000000 to 1000000 · Enter to apply

Intercept w0-1000000 to 1000000 · Enter to applyCoefficient w1-1000000 to 1000000 · Enter to applyCoefficient w2-1000000 to 1000000 · Enter to apply

Rank 3 / 3. All coefficient directions are identified.

Singular values: ≈ 2.9032, ≈ 1.7093, ≈ 0.8061. Cutoff: ≈ 2.9032e−10.

X · augmented designSamples in rows; first column is 1.: 4 × 3
W · coefficientsIntercept first, then feature weights.: 3 × 1
y, predictions, residuals: 4 × 1
NormalizationNo half factor; evaluation excludes penalties.: MSE / n

W = [≈ 0, ≈ 0, ≈ 0]ᵀ

Live state

Data MSE: ≈ 6.5
Numerical rank: ≈ 3
Singular-value cutoff: ≈ 2.9032e−10

\hat y=XW,\quad r=\hat y-y

L_{\mathrm{data}}=\frac1{4}\sum_i r_i^2\approx 6.5

W_*=\sum_{\sigma_j>\tau}v_j\frac{u_j^\mathsf{T}y}{\sigma_j}

Exact inputs; approximate computed results from the active state
ID	Index	Design row	Target	Prediction	Residual
s1	1	[1, -1, 0]	-1	≈ 0	≈ 1
s2	2	[1, 0, 1]	0	≈ 0	≈ 0
s3	3	[1, 1, 0]	3	≈ 0	≈ −3
s4	4	[1, 2, 1]	4	≈ 0	≈ −4

Change the controls first, then inspect the active equation, numeric rows and shapes. Exact inputs and approximate computed values are distinguished; Reset restores the hand-authored experiment.

Make the design do visible work

The opening sheet begins with three zero coefficients, so every prediction is zero. Edit an $x_1$ or $x_2$ design cell or a target, and inspect the changed prediction/residual rows. Editing a design cell changes the data used by every later solve; editing a coefficient changes the current model.

Press Solve minimum norm. For the initial design, the fit gives $W=[1,2,-1]^\mathsf T$ and zero residuals up to floating-point rounding. Then press Duplicate feature column and solve again. Rank falls from three to two. The sheet reports the ambiguity and selects one explicit minimum norm coefficient column.

Define the augmented model and its shapes

Each sample is a row, and each parameter vector is a column:

X=\begin{bmatrix}1&-1&0\\1&0&1\\1&1&0\\1&2&1\end{bmatrix}\in\mathbb R^{4\times3}, \quad W=\begin{bmatrix}b\\w_1\\w_2\end{bmatrix}\in\mathbb R^{3\times1}, \quad y=\begin{bmatrix}-1\\0\\3\\4\end{bmatrix}\in\mathbb R^{4\times1}.

The fixed first column of ones encodes an intercept; it is not another measured feature. Prediction $i$ is the dot product of row $i$ with $W$ . Here the model is $\hat y_i=b+w_1x_{i1}+w_2x_{i2}$ , while $\hat y=XW$ and residual $r=XW-y$ both have shape $4\times1$ .

These are hand-authored unitless toy rows with no random seed. The observation plot uses sample index horizontally and target/prediction vertically. It does not pretend that two features form a scalar horizontal input. The actual features remain in the design table.

Derive the least-squares condition

Full MSE is $L(W)=\|XW-y\|_2^2/n$ . Differentiating each coefficient gives

\nabla_WL=\frac2nX^\mathsf T(XW-y)\in\mathbb R^{d\times1}.

At an optimum, $X^\mathsf T(XW-y)=0$ : the residual is orthogonal to every design column. This stationary condition characterizes the fit. It does not require forming $(X^\mathsf TX)^{-1}$ . Squaring the design’s condition number and then trying to invert a singular or nearly singular matrix is a poor numerical implementation.

Instead decompose $X=U\Sigma V^\mathsf T$ and retain numerically resolvable singular directions:

W_* = \sum_{\sigma_j>\tau}v_j\frac{u_j^\mathsf Ty}{\sigma_j},\qquad \tau=10^{-10}\sigma_{\max}.

Zero or discarded singular directions receive zero coefficient in singular coordinates. This chooses the minimum Euclidean norm among fits to the retained design. Full column rank gives a unique parameter solution. Rank deficiency identifies predictions but leaves a family of parameter columns. Numerical rank is a declared tolerance policy, not an exact algebraic theorem about noisy measurements.

Work every row

For $W=[1,2,-1]^\mathsf T$ , the four dot products are

1+2(-1)-1(0)=-1,\quad 1+2(0)-1(1)=0, \quad 1+2(1)-1(0)=3,\quad 1+2(2)-1(1)=4.

Thus the canonical data MSE changes from $(1+0+9+16)/4=6.5$ at zero coefficients to approximately zero. This exact fit reflects how this toy dataset was authored, not a general property of multivariate regression.

If two columns are identical, only the sum of their coefficients matters. For $X=[[1,1],[2,2]]$ and targets $[2,4]^\mathsf T$ , every pair with $w_1+w_2=2$ fits exactly. Minimum norm gives $[1,1]^\mathsf T$ . A one-row, three-column design $[1,1,1]$ targeting $3$ similarly selects $[1,1,1]^\mathsf T$ from infinitely many fits.

Inspect rank before interpreting coefficients

The browser and canonical Python use a small scaled one-sided Jacobi SVD: orthogonalize design columns through rotations, recover singular values and right singular directions, then apply the retained-direction formula. They never invert normal equations. A distinct oracle uses NumPy’s robust lstsq with the same explicit cutoff, including wide, dependent, nearly dependent and very small/large-scale designs.

The first intercept column remains fixed in this sheet. Numeric fields label their row and feature, preserve exact input text, and accept Enter/arrow/touch operation. After solving, change a coefficient to see that current predictions can leave the optimal solution while numerical rank stays a property of the design.

Mistakes and numerical limits

A feature coefficient is conditional on the other columns and their units. Duplicate or correlated features make individual coefficient interpretation fragile; a perfect residual does not establish causality. Omitting the intercept changes the model family. Transposing $X$ would confuse samples with features and break $XW$ shapes.

The public pure solver supports 1–64 rows and 1–7 columns. Internal Ridge augmentation is trusted separately so 64 observations can still acquire up to seven penalty rows. Inputs must be finite and rectangular. The solver scales the design, uses relative column-orthogonality tolerance $10^{-14}$ and a 100-sweep budget, then rejects nonconvergence or nonfinite results. It reports rank and the cutoff rather than silently regularizing a near-singular direction.

At a singular value extremely near the cutoff, a small perturbation can change numerical rank; inspect singular values and data precision. Parity tolerance is $2\cdot10^{-9}$ absolute plus relative for supported small fixtures. This does not certify arbitrary ill-conditioned data. The interactive table fixes four rows and two editable features; the kernel and oracle cover wider shapes.

Connections and checkpoints

This is Book II’s SVD/projection machinery applied to fitting. Polynomial features in III-11 change the design columns while retaining this linear solve. You should be able to account for the intercept, calculate a row prediction, identify rank ambiguity and explain why a minimum norm solve remains defined when inverse normal equations do not.

Execute the actual Python and NumPy

The imported source below is the runnable file used in canonical comparisons. Read the named mechanism from this chapter; the other functions coordinate the same Part experiment. From the book root:

.venv/bin/python -B code/python/p02/fit_reference.py

code/python/p02/fit_reference.py

"""III-08–12: centered fits, frozen regression updates, small SVD and explicit penalties.

Run from the book root: .venv/bin/python -B code/python/p02/fit_reference.py
Rows are samples; W and y are columns. Intercept is first in augmented X/W.
"""
import json
import math
import numpy as np
from linear_regression import CANONICAL_REGRESSION, derive_regression, regression_trace

RCOND = 1e-10
POLYNOMIAL_DATA = {
    "trainX": [-1, -.6, -.2, .2, .6, 1], "trainY": [1.1, .46, 0, .16, .26, 1.2],
    "heldX": [-.8, -.4, 0, .4, .8], "heldY": [.64, .16, 0, .16, .64],
}


def finite(values, name):
    if not np.isfinite(values).all():
        raise ValueError(f"{name} must remain finite")
    return values


def design(X, y):
    # Validate before coercion: booleans, strings and ragged rows are not numeric data.
    if not isinstance(X, (list, tuple)) or not 1 <= len(X) <= 64:
        raise ValueError("X requires 1–64 rows")
    if not isinstance(X[0], (list, tuple)) or not 1 <= len(X[0]) <= 7:
        raise ValueError("X requires 1–7 columns")
    d = len(X[0])
    for row in X:
        if not isinstance(row, (list, tuple)) or len(row) != d:
            raise ValueError("X must be rectangular")
        for v in row:
            if isinstance(v, bool) or not isinstance(v, (int, float)) or not math.isfinite(v):
                raise ValueError("X must contain finite numbers")
    column(y, len(X), "Targets")
    return finite(np.array(X, dtype=float), "X"), np.array(y, dtype=float).reshape(-1, 1)


def column(values, size, name):
    if not isinstance(values, (list, tuple)) or len(values) != size:
        raise ValueError(f"{name} requires {size} entries (a column)")
    if any(isinstance(v, bool) or not isinstance(v, (int, float)) or not math.isfinite(v) for v in values):
        raise ValueError(f"{name} must contain finite numbers")
    return np.array(values, dtype=float).reshape(-1, 1)


def linear(X, y, weights):
    A, target = design(X, y)
    W = column(weights, A.shape[1], "Weights")
    with np.errstate(over="ignore", invalid="ignore", under="ignore"):
        prediction = finite(A @ W, "Predictions")
        residual = finite(prediction - target, "Residuals")
        mse = float(finite(np.sum(residual * residual) / len(A), "MSE"))
        gradient = finite(2 / len(A) * A.T @ residual, "Gradient")
    n, d = A.shape
    return {"predictions": prediction[:, 0].tolist(), "residuals": residual[:, 0].tolist(),
            "mse": mse, "gradient": gradient[:, 0].tolist(),
            "shapes": {"X": [n, d], "W": [d, 1], "y": [n, 1], "prediction": [n, 1], "gradient": [d, 1]}}


def _least_norm_validated(A, target):
    # Trusted Ridge augmentation may add ≤7 rows beyond the 64 public observation rows.
    d = A.shape[1]
    scale = float(np.max(np.abs(A)))
    if scale == 0:
        return {"weights": [0.] * d, "rank": 0, "singularValues": [0.] * d, "cutoff": 0., "condition": None}
    B = A / scale
    V = np.eye(d)
    for sweep in range(100):
        rotated = False
        for p in range(d):
            for q in range(p + 1, d):
                a, b, c = float(B[:, p] @ B[:, p]), float(B[:, q] @ B[:, q]), float(B[:, p] @ B[:, q])
                if a == 0 or b == 0 or abs(c) <= 1e-14 * math.sqrt(a) * math.sqrt(b):
                    continue
                tau = (b - a) / (2 * c)
                t = 1. if tau == 0 else math.copysign(1., tau) / (abs(tau) + math.hypot(1., tau))
                cs, sn = 1 / math.hypot(1., t), t / math.hypot(1., t)
                for matrix in (B, V):
                    left, right = matrix[:, p].copy(), matrix[:, q].copy()
                    matrix[:, p], matrix[:, q] = cs * left - sn * right, sn * left + cs * right
                rotated = True
        if not rotated:
            break
    else:
        raise ValueError("SVD did not converge in 100 sweeps; no fit returned")
    norms = np.array([math.hypot(*B[:, j]) for j in range(d)])
    cutoff = RCOND * float(np.max(norms))
    weights = np.zeros(d)
    for j in range(d):
        if norms[j] > cutoff:
            coefficient = float((B[:, j] / norms[j]) @ target[:, 0]) / norms[j] / scale
            weights += V[:, j] * coefficient
    finite(weights, "Fitted weights")
    residual = finite(A @ weights.reshape(-1, 1) - target, "Fitted residuals")
    finite(np.sum(residual * residual) / len(A), "Fitted MSE")
    values = sorted(finite(norms * scale, "Singular values").tolist(), reverse=True)
    rank = int(np.sum(norms > cutoff))
    condition = values[0] / values[-1] if rank == d else None
    return {"weights": weights.tolist(), "rank": rank, "singularValues": values,
            "cutoff": cutoff * scale, "condition": condition}


def least_norm(X, y):
    A, target = design(X, y)
    return _least_norm_validated(A, target)


def closed_form(samples):
    derive_regression({"samples": samples, "weight": 0., "bias": 0., "learningRate": .1})
    x = np.array([r["x"] for r in samples], dtype=float)
    y = np.array([r["y"] for r in samples], dtype=float)
    with np.errstate(over="ignore", invalid="ignore", under="ignore"):
        mean_x, mean_y = float(finite(np.mean(x), "Mean x")), float(finite(np.mean(y), "Mean y"))
        dx, dy = finite(x - mean_x, "Centered x"), finite(y - mean_y, "Centered y")
        sxx, sxy = float(finite(dx @ dx, "Sxx")), float(finite(dx @ dy, "Sxy"))
    if sxx == 0 and any(v != x[0] for v in x):
        raise ValueError("Input variance underflowed; slope is numerically unresolved")
    weight = 0. if sxx == 0 else sxy / sxx
    bias = mean_y - weight * mean_x
    mse = derive_regression({"samples": samples, "weight": weight, "bias": bias, "learningRate": .1})["mse"]
    return {"weight": weight, "bias": bias, "meanX": mean_x, "meanY": mean_y, "sxx": sxx, "sxy": sxy, "mse": mse, "unique": sxx != 0}


def descent(state):
    result = derive_regression(state)  # Frozen III-07 is the predictor, gradient and update authority.
    x = np.array([r["x"] for r in state["samples"]])
    a, b, c = 2 * float(np.mean(x * x)), 2 * float(np.mean(x)), 2.
    maximum = (a + c + math.hypot(a - c, 2 * b)) / 2
    factor = 1 - state["learningRate"] * maximum
    behavior = ("undamped oscillation" if factor == -1 else "divergence") if factor <= -1 else ("convergent oscillation" if factor < 0 else "convergence without fast-mode sign flips")
    return {**result, "lambdaMax": maximum, "criticalRate": 2 / maximum, "factor": factor, "behavior": behavior}


def polynomial_design(x, degree):
    if type(degree) is not int or not 1 <= degree <= 5:
        raise ValueError("Degree must be an integer from 1 to 5")
    column(x, len(x), "Polynomial input")
    return finite(np.array([[v ** j for j in range(degree + 1)] for v in x]), "Polynomial features").tolist()


def polynomial_fit(samples, degree):
    closed_form(samples)
    X, y = polynomial_design([r["x"] for r in samples], degree), [r["y"] for r in samples]
    solve = least_norm(X, y)
    return {**solve, "train": linear(X, y, solve["weights"]), "held": linear(polynomial_design(POLYNOMIAL_DATA["heldX"], degree), POLYNOMIAL_DATA["heldY"], solve["weights"]), "X": X}


def regularized(X, y, lam, kind, penalize_intercept):
    A, target = design(X, y)
    if isinstance(lam, bool) or not isinstance(lam, (int, float)) or not math.isfinite(lam) or not 0 <= lam <= 10 or kind not in ("l1", "l2") or type(penalize_intercept) is not bool:
        raise ValueError("Penalty requires L1/L2, lambda in [0,10] and an explicit intercept policy")
    if not np.all(A[:, 0] == 1):
        raise ValueError("First design column must be the intercept ones")
    n, d = A.shape
    mask = np.array([j != 0 or penalize_intercept for j in range(d)])
    iterations = 0
    if kind == "l2" or lam == 0:
        extra = np.diag(np.full(d, math.sqrt(n * lam)))[mask] if lam > 0 else np.empty((0, d))
        solve = _least_norm_validated(np.vstack([A, extra]), np.vstack([target, np.zeros((len(extra), 1))]))
        weights = np.array(solve["weights"])
    else:
        weights = np.zeros(d)
        for iterations in range(1, 3001):
            old = weights.copy()
            for j in range(d):
                partial = target[:, 0] - A @ weights + A[:, j] * weights[j]
                rho, curvature = float(A[:, j] @ partial / n), float(A[:, j] @ A[:, j] / n)
                threshold = lam / 2 if mask[j] else 0
                weights[j] = 0 if curvature == 0 else math.copysign(max(0, abs(rho) - threshold), rho) / curvature
            finite(weights, "L1 weights")
            if np.max(np.abs(weights - old)) < 1e-12:
                break
    result = linear(X, y, weights.tolist())
    penalty = float(finite(lam * np.sum((weights ** 2 if kind == "l2" else np.abs(weights))[mask]), "Penalty"))
    gradient = np.array(result["gradient"])
    kkt = []
    for j, w in enumerate(weights):
        g = gradient[j]
        kkt.append(abs(g) if not mask[j] else abs(g + 2 * lam * w) if kind == "l2" else max(0, abs(g) - lam) if w == 0 else abs(g + lam * np.sign(w)))
    residual = float(max(kkt))
    return {**result, "weights": weights.tolist(), "penalty": penalty, "objective": float(finite(result["mse"] + penalty, "Objective")), "kkt": residual, "iterations": iterations, "converged": residual <= 1e-8}


def demonstration():
    samples = CANONICAL_REGRESSION["samples"]
    poly = [dict(id=f"s{i+1}", x=x, y=POLYNOMIAL_DATA["trainY"][i]) for i, x in enumerate(POLYNOMIAL_DATA["trainX"])]
    return {"closedForm": closed_form(samples), "descent": descent(CANONICAL_REGRESSION),
            "descentTrace": regression_trace(CANONICAL_REGRESSION, 4),
            "multiple": least_norm([[1, -1, 0], [1, 0, 1], [1, 1, 0], [1, 2, 1]], [-1, 0, 3, 4]),
            "polynomial": {str(d): polynomial_fit(poly, d) for d in (1, 2, 5)},
            "regularization": {kind: regularized([[1, r["x"]] for r in samples], [r["y"] for r in samples], .5, kind, False) for kind in ("l1", "l2")}}


if __name__ == "__main__":
    print(json.dumps(demonstration(), indent=2, allow_nan=False))