PCA: Derive, Fit, and Reuse a Projection

Two sensors may move together so closely that recording both adds little variation. PCA replaces their original coordinates with directions ordered by the variance they capture. Keeping fewer directions minimizes squared reconstruction error among linear orthogonal projections of the centered training data. It does not guarantee preservation of labels, meaning, or useful predictive information.

Before starting, be comfortable with dot products, matrix shapes, means, and variance. This module connects those tools to the pca-implementation exercise using Python lists. By the end, you should be able to derive a component, explain deflation, and transform a new row without fitting again.

Fix the exercise conventions first

Rows are observations. For n rows, d features, and k retained components:

Object	Shape	Meaning
X	(n,d)	Original observations
mean	(d,)	Training feature means
C	(d,d)	Covariance of centered rows
components	(k,d)	One unit direction per row
Z	(n,k)	Coordinates in the retained directions

The exercise returns (Z, components, mean). It uses C = X_centered.T @ X_centered / n. Statistics texts often divide by n-1; for n>1 that rescales all eigenvalues by n/(n-1) without changing directions or explained-variance ratios. Follow the exercise's /n convention when implementing or checking its numbers.

The formulas below use @ as mathematical notation. Python lists require the dot-product and matrix-vector helpers in the starter, not NumPy.

Trace a complete fit by hand

Use training rows [[1,2], [3,4], [5,6]].

mean       = [3,4]
centered   = [[-2,-2], [0,0], [2,2]]
C          = [[8/3,8/3], [8/3,8/3]]

Each covariance entry is a sum of paired centered feature values divided by three. Both sensors vary together, so every entry is 8/3.

For a unit direction v, a centered row's coordinate is its dot product with v. Its mean squared coordinate is v.T @ C @ v. Maximizing that expression chooses the first component. A covariance eigenvector satisfies C @ v = lambda*v; for a unit eigenvector, lambda is the variance along that direction.

Here, choose:

v1         = [1/sqrt(2), 1/sqrt(2)]
lambda1    = 16/3
Z          = [[-2*sqrt(2)], [0], [2*sqrt(2)]]

Check the variance: ((−2*sqrt(2))² + 0² + (2*sqrt(2))²)/3 = 16/3. The second perpendicular direction [1/sqrt(2), -1/sqrt(2)] has zero variance. Every row can be reconstructed exactly from its first coordinate: z*v1 + mean.

Pause. What changes if the implementation returns -v1? Answer: all first coordinates change sign too. Reconstruction and captured variance stay the same. A test should compare directions up to sign or check their projection, not demand identical signs.

Find a direction with power iteration

The exercise asks you to approximate an eigenvector without calling eig or SVD. Start with a unit vector, repeatedly multiply by C, and normalize:

v <- v / ||v||
repeat:
    candidate <- C @ v
    v <- candidate / ||candidate||
lambda <- dot(v, C @ v)

Why does this work? Resolve the initial vector into eigenvector directions. Every multiplication scales each contribution by its eigenvalue. When the largest eigenvalue is strictly larger than the next and the start has a nonzero contribution in its direction, that contribution eventually dominates.

For C = [[4,0],[0,1]], start with [1,1]/sqrt(2). The next directions are [4,1]/sqrt(17), then [16,1]/sqrt(257). The second-to-first coordinate ratio shrinks from 1 to 1/4 to 1/16. Normalization keeps the vector's length bounded; it does not change this ratio.

The approximation is slower when the two largest eigenvalues are close. Equal largest eigenvalues permit any unit direction in their shared eigenspace. A start exactly perpendicular to the dominant eigenspace cannot discover it, so a seeded random start is useful. This teaching solver uses a finite iteration budget; it is not a general accuracy guarantee for ill-conditioned matrices.

A zero product needs a separate branch: do not divide by zero. For an all-zero covariance matrix, every unit direction is an eigenvector with eigenvalue zero.

Remove captured variance, then keep directions distinct

After extracting unit vector v and eigenvalue lambda, deflation removes its contribution:

C_remaining[i][j] = C[i][j] - lambda*v[i]*v[j]

For C=diag(4,1), v=[1,0], and lambda=4, the remaining matrix is diag(0,1). The next component is the second axis.

Exact deflation of exact eigenvectors gives orthogonal subsequent directions. Floating-point deflation can leave tiny residuals along the first direction. On the sensor example, blindly iterating on those residuals can return the first direction again and count the same variation twice.

During each later iteration, remove the candidate's projections onto previously selected components:

for previous in components:
    candidate <- candidate - dot(candidate, previous)*previous

Then normalize. Repeating that projection pass can reduce accumulated roundoff. Apply it to the starting vector too. If a start disappears, try a coordinate axis with a nonzero remaining projection. When no resolvable variance remains, complete the basis with unit orthogonal directions; their training coordinates are zero to numerical tolerance. Do not return duplicate directions or zero vectors. The exercise accepts any valid orientation of this zero-variance subspace.

Pause. Does requesting two components from rank-one data reveal a second signal? Answer: no. It adds a coordinate direction with zero training variance. A full orthonormal basis is still useful for a well-defined transform of new rows.

Separate fitting from transforming

pca learns the mean and components from training data. transform reuses them:

Z[i][r] = sum((X[i][j] - mean[j]) * components[r][j] for j in features)

For the fitted sensor example, new row [7,8] has centered row [4,4] and coordinate 4*sqrt(2). Recomputing the mean from this one new row would center it to zero and erase its displacement from the training data.

In matrix notation, projection is Z = (X-mean) @ components.T. Reconstruction is Z @ components + mean. Keeping all d orthonormal directions reconstructs each row; keeping fewer generally loses some variation.

Fit PCA only on training data. During cross-validation, fit it separately within each training fold. Otherwise the learned directions reveal held-out data structure to the model-selection process.

Choose what variance should mean

Centering removes location. Scaling changes the relative importance of feature variation. If meters and millimeters represent otherwise comparable features, standardization can prevent units from deciding the result. If all features share meaningful units and their variance differences matter, automatic standardization would change the intended objective. The exercise centers but does not standardize. The scikit-learn PCA documentation also distinguishes centering from scaling.

With eigenvalues [6,3,1], two components retain (6+3)/10 = 90% of training variance. This percentage is descriptive, not a universal threshold for selecting k. A low-variance feature may contain the entire class signal while high-variance features contain unrelated noise. Validate the downstream task as well as reconstruction. When total variance is zero, a variance ratio is undefined; do not divide zero by zero.

Build and check in stages

1. Implement dot, matvec, and normalize; give the zero vector an explicit normalization convention.
2. Verify centered column sums and the numeric covariance above.
3. Run power iteration on a diagonal matrix before adding deflation.
4. Check component lengths and pairwise dot products, including constant and rank-one data.
5. Verify projection with the saved training mean and reconstruct a known row.

Storage includes a d*d covariance matrix. Computing it costs O(n*d*d) operations. Each dense power iteration costs O(d*d), with extra projection work for previously selected components. For large or difficult datasets, production PCA implementations commonly use SVD or more sophisticated solvers; this exercise exposes the mechanics on small, well-separated examples.