Keep a shadow and measure the missing part

Remove a coordinate and inspect what survives

At k=1, points a and b have the same hollow reconstruction square. Their first scores agree, even though the points differ. The residual segments show which information was discarded. Move either observation and inspect its scores, reconstruction and squared error in the table.

Select Replay k = 0, 1, 2. The first frame keeps no components and reconstructs every observation at the mean. Next keeps the leading direction. The last retains the whole two-feature space. This trace changes the component count; PCA is not taking three optimization iterations. Editing a point, normalization or k during replay starts from the displayed cloud and ends replay.

Reset restores the original cloud with k=1. Try line cloud: a single component suffices in exact arithmetic. Try isotropic cloud: any one line loses the same total energy, although the individual shadows depend on the conventional basis. Try constant cloud: reconstruction error is zero even though retained variance is undefined. A reconstruction can be perfect without a meaningful variance ratio.

Projection changes the representation, not the observation count

Let Vk contain k orthonormal direction columns, so Vk is d×k. Each centered observation has a k×1 score vector $z_i=V_k^T(x_i-\mu)$. The row-batch equation is

The observation count n remains the same. Scores are coordinates in a learned basis, not selected original feature columns. A single principal component can combine every original feature. Reconstructing adds the training mean back into source coordinates; forgetting it would place the cloud around the origin rather than its actual location.

At k=0, Z has n×0 shape. The empty product is zero, so reconstruction is the mean-only baseline. At k=d, VkVkᵀ=I and reconstruction is exactly X in real arithmetic. Both the browser and canonical Python apply that full-space identity directly for k=2, copying the original coordinate values. They still calculate and expose the scores, but avoid an unnecessary floating rotate/unrotate round trip. Thus full-space residuals and SSE are exactly zero in this implementation, rather than accumulated rounding artifacts.

Derive the residual and its optimality

The projector Pk=VkVkᵀ is symmetric and idempotent. Reconstruction of the centered column xci is Pkxci; the residual is

The residual is perpendicular to every retained direction. Pythagoras separates retained and discarded energy. Over the whole batch,

Why is PCA the best rank-k orthogonal reconstruction for this centered batch? A k-dimensional projector retains variance tr(VkᵀCVk), maximized by the leading eigenvectors. Total variance tr(C) is fixed, so maximizing retained variance minimizes the remainder. Equivalently, the SVD keeps the largest singular values and discards the rest. This is a training reconstruction optimum in the declared Euclidean metric, not a guarantee about downstream prediction or future data.

We report SSE, per-observation MSE=SSE/n, and per-entry MSE=SSE/(nd). They answer different questions. The first averages each observation’s squared vector norm; the second averages individual scalar-coordinate errors. For this lab d=2, they differ by a factor of two. None uses n−1: that divisor belongs to sample covariance, not the reconstruction metric.

Work every shadow in the original cloud

For $v_1=(1,1)^T/\sqrt2$ and μ=0, point a=(−2,−1) has

Point b has the same first score and reconstruction, with residual (1/2,−1/2). Points c and d reconstruct at (3/2,3/2), each with squared residual 1/2. Hence SSE=2, per-observation MSE=1/2 and per-entry MSE=1/4. With sample eigenvalues 6 and 2/3, the spectral identity gives qλ2=3(2/3)=2. Population covariance instead gives 4(1/2)=2: the geometric error agrees.

Changing the first eigenvector to its negative flips −3/√2 to +3/√2 and flips the reconstruction direction at the same time. The product is identical. If the spectrum repeats at the truncation boundary, changing the retained line can change the individual reconstructed points while preserving the optimal total SSE.

Use the learned transform on future rows

Fit μ and Vk using training data. A future row xnew uses znew=Vkᵀ(xnew−μ) and reconstructs as μ+Vkznew. Refitting the basis on the future row would be a different operation. The lab refits its displayed descriptive cloud after edits; it does not implement a frozen train/test pipeline or estimate unseen reconstruction performance.

A retained fraction can guide a compression budget, but low-variance directions may contain important predictive information. Consider a class label encoded in a tiny perpendicular offset from a long line. One-component PCA can retain almost all variation while discarding exactly that label signal. Use evaluation suited to the downstream task, rather than assuming a 90% retained fraction is a 90% accuracy promise.

Units matter here too. Raw squared reconstruction error adds squared errors from each feature. Changing one feature’s numerical scale changes the objective and can change the chosen basis. Standardization and inverse scaling must be explicit parts of a larger pipeline; this Part applies neither automatically. It also does not whiten, invent missing values or claim arbitrary-dimensional browser support.

Numerical and shape caveats

The 2D core accepts 1–128 dense finite observations with unique IDs and integral numeric k in [0,2]. Sparse array holes, ragged rows, NaN, infinity, invalid configuration and computed overflow are rejected before degenerate returns. k=0 still validates every coordinate: “not retaining components” is not permission to accept malformed data.

For a singleton, population covariance has zero spread and reconstruction is the observed mean. Sample covariance is undefined, and the lab withholds PCA outputs instead of fabricating them. For a constant batch, covariance is zero, reconstruction error is zero and ratios stay undefined at every k. There is no unique principal direction for zero spread.

Execute the actual Python and NumPy

The reference pca(request) returns each row’s scores, fullScores, reconstructed, residual and squaredError, plus batch SSE and both MSE conventions. scores has k entries; fullScores contains the two inspectable coordinates before truncation. axes is a list of direction columns. The kernel deliberately handles k=2 as the identity and k=0 as the empty-basis mean baseline.

.venv/bin/python -B code/python/p06/pca.py

That terminal invocation uses the four original rows. Complete JSON on stdin selects another cloud or k. Independent committed literals include the k=0/1/2 cloud, a rank-one line, repeated spectrum, zero spread, one-sample conventions and a tiny positive minor eigenvalue. The independent Fraction/SVD/eigh oracle executes the canonical module and checks actual reconstruction energy; it does not trust copied browser output.

"""Inspectable two-feature PCA. X stores observations in rows; directions are columns.

No scaling or whitening. Sample normalization divides by n-1; population by n.
One sample has undefined sample covariance. Zero spread has undefined variance ratios.
"""
import json
import math
import sys
import numpy as np

CLOUD = [{"id": i, "x": x} for i, x in zip("abcd", [[-2, -1], [-1, -2], [1, 2], [2, 1]])]
DEFAULT = dict(points=CLOUD, normalization="sample", k=1)


def finite(value, label):
    if isinstance(value, (bool, np.bool_)) or not isinstance(value, (int, float, np.integer, np.floating)):
        raise ValueError(f"{label} must be a finite real number.")
    try:
        number = float(value)
    except (OverflowError, ValueError) as error:
        raise ValueError(f"{label} must be finite; binary64 conversion overflow is rejected.") from error
    if not math.isfinite(number):
        raise ValueError(f"{label} must be finite; arithmetic overflow is rejected.")
    return number + 0.0


def energy(values):
    total = 0.0
    for value in values:
        square = finite(float(value) * float(value), "Squared component")
        if value != 0 and square == 0:
            raise ValueError("Squared component underflow: rescale the observations.")
        total = finite(total + square, "Sum of squares")
    return total


def normalized_energy(values, divisor):
    total = energy(values)
    value = finite(total / divisor, "Normalized energy")
    if total > 0 and value == 0:
        raise ValueError("Variance underflow: rescale the observations.")
    return value


def pca(request):
    if not isinstance(request, dict):
        raise ValueError("PCA request must be an object.")
    normalization, k, points = request.get("normalization"), request.get("k"), request.get("points")
    if normalization not in ("sample", "population"):
        raise ValueError("Normalization must be sample or population.")
    k = finite(k, "Component count")
    if not k.is_integer() or not 0 <= k <= 2:
        raise ValueError("Component count must be an integer from 0 to 2.")
    k = int(k)
    if not isinstance(points, list) or not 1 <= len(points) <= 128:
        raise ValueError("Use 1 to 128 dense two-feature observations.")
    ids, coordinates = [], []
    for point in points:
        if not isinstance(point, dict) or not isinstance(point.get("id"), str) or not point["id"].strip() or point["id"] in ids:
            raise ValueError("Point IDs must be nonempty unique strings.")
        x = point.get("x")
        if not isinstance(x, list) or len(x) != 2:
            raise ValueError("Each observation must have exactly two dense components.")
        ids.append(point["id"])
        coordinates.append([finite(value, "Observation component") for value in x])
    X = np.asarray(coordinates, dtype=np.float64)  # n×2; boundary checked before degenerate returns
    n = len(points)
    mean = []
    for j in range(2):
        first = float(X[0, j])
        if all(float(x) == first for x in X[:, j]):
            mean.append(first)
        else:
            delta = 0.0
            for x in X[:, j]:
                delta = finite(delta + finite(float(x) - first, "Mean difference") / n, "Mean accumulation")
            mean.append(finite(first + delta, "Mean"))
    centered = np.asarray([[finite(float(x[j]) - mean[j], "Centering") for j in range(2)] for x in X])
    divisor = n - 1 if normalization == "sample" else n
    products = [[finite(float(x[0]) * float(x[0]), "Squared x"), finite(float(x[0]) * float(x[1]), "Cross product"), finite(float(x[1]) * float(x[1]), "Squared y")] for x in centered]
    basic = dict(n=n, divisor=divisor, mean=mean, centered=centered.tolist(), products=products)
    if divisor == 0:
        return dict(basic, status="undefined", reason="Sample covariance needs n ≥ 2; n−1 is zero.")
    a, c = (normalized_energy(centered[:, j], divisor) for j in range(2))
    b = 0.0
    for x in centered:
        b = finite(b + finite(float(x[0]) * float(x[1]), "Cross product"), "Covariance sum")
    b = finite(b / divisor, "Covariance")
    repeated, scale = a == c and b == 0, max(a, c, abs(b))
    v1, large = [1.0, 0.0], 0.0
    if scale > 0:
        A, B, C = a / scale, b / scale, c / scale
        L = (A + C) / 2 + math.hypot((A - C) / 2, B)
        large = finite(scale * L, "Largest eigenvalue")
        if b == 0:
            v1 = [1.0, 0.0] if a >= c else [0.0, 1.0]
        else:
            raw = [L - C, B] if A >= C else [B, L - A]
            norm = math.hypot(*raw)
            if norm == 0:
                raise ValueError("Eigenvector scaling underflow: direction is unresolved in binary64.")
            sign = -1 if raw[0] < 0 or (raw[0] == 0 and raw[1] < 0) else 1
            v1 = [sign * x / norm for x in raw]
    v2 = [-v1[1], v1[0]]
    axes = [v1, v2]  # serialized list of 2×1 columns, not a row-batch basis matrix
    def dot(x, v):
        return finite(finite(float(x[0]) * v[0], "Dot product") + finite(float(x[1]) * v[1], "Dot product"), "Dot product")
    scores = [[dot(x, v) for v in axes] for x in centered]
    small = normalized_energy([z[1] for z in scores], divisor)
    if small > large:  # last-bit ordering in a nearly repeated spectrum
        large, small = small, large
        axes = [v2, v1]
        scores = [[z[1], z[0]] for z in scores]
        v1 = axes[0]
    eigenvalues = [large, small]
    total = finite(large + small, "Total variance")
    ratios = None if total == 0 else [v / total for v in eigenvalues]
    rows = []
    for i, x in enumerate(X):
        z = scores[i]
        # Full-space identity is applied directly, avoiding needless rotate/unrotate roundoff.
        reconstructed = x.tolist() if k == 2 else [finite(mean[j] + (finite(z[0] * v1[j], "Reconstruction product") if k else 0), "Reconstruction") for j in range(2)]
        residual = [finite(float(x[j]) - reconstructed[j], "Residual") for j in range(2)]
        rows.append(dict(id=ids[i], x=x.tolist(), centered=centered[i].tolist(), scores=z[:k], fullScores=z,
                         reconstructed=reconstructed, residual=residual, squaredError=energy(residual)))
    sse = 0.0
    for row in rows:
        sse = finite(sse + row["squaredError"], "Reconstruction SSE")
    return dict(basic, status="defined", covariance=[[a, b], [b, c]], axes=axes, eigenvalues=eigenvalues,
                totalVariance=total, ratios=ratios, retained=None if ratios is None else sum(ratios[:k]),
                repeated=repeated, nearRepeated=not repeated and abs(large - small) <= 64 * sys.float_info.epsilon * total,
                zeroVariance=total == 0, rows=rows, sse=sse, mse=sse / n, entryMse=sse / (2 * n))


if __name__ == "__main__":
    try:
        request = DEFAULT if sys.stdin.isatty() else json.load(sys.stdin)
        print(json.dumps(pca(request), allow_nan=False))
    except (ValueError, TypeError, OverflowError) as error:
        print(json.dumps({"error": str(error)}), file=sys.stderr)
        sys.exit(2)

Five graded investigations

Carry the mechanism forward

You can compute scores and source-space reconstructions, prove residual orthogonality, relate discarded eigenvalues to SSE, distinguish both MSE denominators and handle zero or repeated spectra honestly. III-29 established centered covariance; III-30 selected its directions. This completes the geometry of dimensionality reduction and prepares you to ask how a representation affects later supervised evaluation.