SVD, low-rank approximation, and PCA

Eigenvectors compare a square map’s input and output in the same space. A rectangular design, image matrix or embedding table instead connects different spaces. Even a square nonsymmetric map can lack an orthonormal eigenbasis. The singular value decomposition supplies orthonormal input and output directions for every real finite matrix, with nonnegative strengths connecting them. It extends the useful geometry from the symmetric spectral theorem without claiming that arbitrary eigenvectors are orthogonal.

This module derives that decomposition, the pseudoinverse and best-rank approximation, then applies them to centred principal component analysis. You will distinguish full, reduced and compact shapes, calculate discarded reconstruction energy, and fit preprocessing on training records before applying it elsewhere. High explained variance measures a particular geometric objective; it does not establish predictive usefulness.

Retrieval check: state the real symmetric spectral theorem and explain why AᵀA is positive semidefinite. Revisit Module 13. Projection, minimum-norm selection and the four spaces come from Modules 12 and 11. The PCA scatter calculations here are algebraic; probability and statistical interpretation return in Module 23. Use NumPy preparation for labs.

For A of shape m×n, AᵀA is an n×n symmetric positive-semidefinite matrix. The spectral theorem supplies an orthonormal basis v₁,…,v_n with nonnegative eigenvalues, sorted from largest to smallest. Define σᵢ as their nonnegative square roots. For each positive σᵢ, set uᵢ=Avᵢ/σᵢ. Then uᵢᵀuⱼ=vᵢᵀAᵀAvⱼ/(σᵢσⱼ), which equals one for i=j and zero otherwise. Thus the nonzero image directions are orthonormal. A zero eigenvalue gives ‖Avᵢ‖²=0, hence Avᵢ=0.

Expanding any input x in the full v basis now gives Ax=Σ_{σᵢ>0}σᵢuᵢ(vᵢᵀx). This proves A=Σ_{i=1}^rσᵢuᵢvᵢᵀ, where r is the number of positive values. The u’s span the column space and the corresponding v’s span the row space; the remaining v’s span the kernel. Extend the u list to an orthonormal basis of Rᵐ. These constructions yield the singular value decomposition, or SVD, A=UΣVᵀ. Its existence follows from the symmetric theorem even if A itself is rectangular or defective.

The full SVD has U of shape m×m, V of shape n×n, and rectangular diagonal Σ of shape m×n. Its first q=min(m,n) diagonal entries are singular values, including zeros when rank is smaller than q. The remaining positions are zero. Full U and V describe complete output and input bases, including directions not used by the map. The formula does not multiply an m×m U by an n×n diagonal array unless m=n; the middle factor’s rectangular shape matters.

The reduced SVD retains q columns on each side: U_q is m×q, V_q is n×q, and the diagonal middle factor is q×q. It still includes zero singular values if r<q. The compact rank-r SVD retains only positive directions: U_r is m×r,V_r is n×r,Σ_r is r×r. Reduced and compact are therefore different when the matrix is rank deficient. Library terminology can vary; specify shapes and whether zero directions are present rather than relying on one word.

The rank-zero case also fits the compact formula: the positive lists are empty and the sum is the zero matrix. Full or reduced factors can still contain arbitrary orthonormal zero-direction bases. Repeated positive singular values permit paired orthonormal rotations of the corresponding u and v lists. Flipping both a u vector and its matching v vector leaves their outer product unchanged. Comparing signed columns without allowing these ambiguities is not an appropriate correctness criterion.

The construction through AᵀA proves existence, but does not prescribe forming that Gram matrix as a numerical SVD algorithm. As Module 12 showed, squaring weak direction strengths can destroy resolution in floating-point intermediates. A library SVD works through specialised numerical routines. Validate the factor shapes, orthogonality and reconstruction, and record dtype and tolerance when deciding which values are effectively zero.

An SVD also separates coordinates from invariant quantities. If orthogonal input and output basis changes replace A by PᵀAQ, the Gram matrix becomes QᵀAᵀAQ and has the same eigenvalues. Singular values, rank and the two matrix norms therefore remain unchanged; singular-vector tuples change with the bases. Multiplying the whole map by a scalar α multiplies every strength by |α|. Scaling just one feature generally changes them in a less uniform way because that operation is not an orthogonal coordinate change. These identities let you distinguish a reorientation that preserves Euclidean geometry from a unit conversion that changes its metric.

The full formula acts in three stages. Vᵀ expresses the input in an orthonormal basis, preserving its Euclidean length. Σ scales the matched coordinates and inserts or removes coordinates as required by the rectangular dimensions. U expresses the resulting output coordinates in an orthonormal output basis, again preserving length. Orthogonal factors may be rotations or reflections; calling them both rotations would omit determinant-minus-one possibilities. The middle diagonal strengths are nonnegative, unlike potentially signed or complex eigenvalues.

For each positive pair, Avᵢ=σᵢuᵢ and Aᵀuᵢ=σᵢvᵢ. The second relation follows from AᵀAvᵢ=σᵢ²vᵢ. Hence uᵢ is also an eigenvector of AAᵀ with value σᵢ². AᵀA and AAᵀ share the positive eigenvalues but generally have different numbers of zeros because their dimensions differ. These are input and output versions of one map’s direction strengths, not a claim that A has the same ordinary eigenvalues.

With rank r<n, some input directions vanish. Their coordinates remain part of a full V basis even though they contribute no output. With r<m, the image occupies a proper subspace of the output space and some output directions are unattainable. The last full U directions span the left null space. This makes all four fundamental spaces visible: the first r right directions span the row space, remaining right directions the null space; the first r left directions span the column space, remaining left directions the left null space.

Singular values directly give the Euclidean operator norm. For a unit input x=Vc, ‖Ax‖²=Σ_{i=1}^rσᵢ²cᵢ²≤σ₁²Σcᵢ²=σ₁². The first right singular vector attains the bound, so ‖A‖₂=σ₁, with norm zero for the zero map. For a square invertible matrix, the smallest value bounds the least stretching and the inverse norm is 1/σ_n. Therefore κ₂(A)=σ₁/σ_n without a symmetry assumption. This repairs the unqualified eigenvalue-ratio shortcut diagnosed in Module 13.

For a full-column-rank rectangular design X, the ratio σ₁/σ_d measures its least-to-greatest input stretching and the norm of its pseudoinverse. Its Gram matrix has eigenvalues σᵢ², so κ₂(XᵀX)=(σ₁/σ_d)². This proves the conditioning-square mechanism introduced earlier. A rank-deficient design has zero input directions and no finite inverse sensitivity over the full input space, although one can discuss a ratio restricted to positive singular directions if that restriction is stated.

Storage conventions deserve the same care as geometry. Real NumPy output has s as a one-dimensional array, U as left vectors in columns, and Vh as right vectors in rows. Multiplying U*s broadcasts the values across U’s columns when U has the matching reduced shape. Scaling its rows would represent a different matrix. For complex inputs the returned Vh is a conjugate transpose; this course’s main examples are real, where ordinary transpose suffices. A reconstruction test with the expected shape detects many orientation mistakes before they become misleading interpretations.

For compact A=U_rΣ_rV_rᵀ, define A⁺=V_rΣ_r⁻¹U_rᵀ, where only positive strengths are inverted. Its shape is n×m. It projects an output target onto the column space by taking u coordinates, divides each retained coordinate by its strength, and constructs an input in the row space. This is precisely Module 12’s geometric minimum-norm least-squares map. An ordinary inverse appears as the special square full-rank case, not as a requirement for the formula.

Multiplying gives AA⁺=U_rU_rᵀ and A⁺A=V_rV_rᵀ, the column- and row-space projections. They are symmetric and idempotent. Consequently AA⁺A=A and A⁺AA⁺=A⁺, recovering all Moore–Penrose conditions. For any y, w_min=A⁺y has prediction U_rU_rᵀy, the unique closest reachable output. Every least-squares minimiser is w_min+n with n in the kernel, and orthogonality makes w_min uniquely smallest in Euclidean norm.

A small positive σ can amplify target noise: a disturbance εuᵢ produces coefficient change (ε/σᵢ)vᵢ when that direction is retained. The input and output norms of those changes are |ε|/σᵢ and |ε| respectively. This is an exact directional sensitivity statement. It explains why a numerical algorithm can be accurate for the supplied targets and yet return unstable coefficients after a measurement perturbation. Singular vectors make the mechanism visible rather than hiding it behind one scalar condition number.

For diag(1,10⁻¹⁰), the target (1,10⁻¹⁰) gives inverse weights (1,1). Adding (0,10⁻⁸) changes them to (1,101). A numerical pseudoinverse that discards the second direction returns (1,0) instead and leaves a second-coordinate residual. It has deliberately altered the effective inverse problem to limit amplification. That choice is not a proof that the exact second singular value vanished. It trades resolution and fit against parameter response, requiring a stated tolerance and purpose.

The pinned pseudoinverse interface uses a relative cutoff based on the largest singular value. An absolute rank threshold, a relative pseudoinverse threshold and a measurement-noise policy are not automatically interchangeable. Global scaling interacts with absolute thresholds differently from relative ones, and scaling individual columns also changes the geometry. Report the positive directions retained, the feature units, dtype and tolerance. If a value is close to the cutoff, small input changes can change the selected rank and the returned representative.

Do not form 1/s blindly when the reduced list includes zeros. Even “multiply by zero afterwards” is unsafe because an infinite intermediate can yield NaN. Build a retained-direction mask or use a library with an explicit cutoff policy. An all-zero matrix has a zero pseudoinverse; there is no positive direction to divide by. A meaningful verification checks the selected projected output and minimum-norm geometry, not just that a numerical call finished without an exception.

Keep the k strongest terms to define A_k=Σ_{i=1}^kσᵢuᵢvᵢᵀ, with 0≤k≤r. Its rank is k for positive retained values, and A₀=0. This is a low-rank approximation in the original matrix shape. Keeping factors rather than constructing all entries can save storage or computation when k is small. The factor scalar count is k(m+n+1) if both direction factors and k strengths are stored. This must be compared with mn and with actual file-format overhead; low rank alone does not guarantee a smaller stored file.

The Frobenius norm is ‖M‖F=√ΣᵢⱼMᵢⱼ², equivalently the square root of the sum of squared column lengths. Left multiplication by an orthogonal matrix preserves each column length; right orthogonal multiplication also preserves the total, by applying the same reasoning to the transpose. Hence Frobenius norm is unchanged by full orthogonal basis changes. In the SVD coordinates, A−A_k has only discarded strengths on its diagonal, so ‖A−A_k‖F²=Σ{i>k}σᵢ². The operator norm error is the largest discarded value σ{k+1}, with zero after every positive direction is kept.

The best-rank-k approximation theorem says no matrix B of rank at most k has smaller Frobenius or operator error than A_k. Here is a complete Frobenius proof. Let P project onto C(B), of dimension t≤k. In each column, the error splits into its component in that space and the perpendicular component, so ‖A−B‖_F²≥‖(I−P)A‖_F²=‖A‖_F²−‖PA‖_F². Using the left singular basis gives ‖PA‖_F²=Σσᵢ²wᵢ, where wᵢ=‖Puᵢ‖² lies between zero and one. Include all full left-basis vectors, assigning zero strength to unused directions. Their weights sum to t, because summing squared projections of an orthonormal basis onto a t-dimensional space counts t unit directions.

The weighted strength sum is at most Σ_{i≤k}σᵢ². To see this when k>0, bound every tail strength by σ_k² and use the remaining weight budget k−Σ_{i≤k}wᵢ. Then the sum is at most Σ_{i≤k}[σᵢ²wᵢ+σ_k²(1−wᵢ)], no greater than Σ_{i≤k}σᵢ². For k=0 the projection is zero. Thus every candidate error is at least the discarded square sum, which A_k attains. The argument handles arbitrary candidate images rather than assuming their directions already align with A.

For operator error with k<r, restrict B to the span of the first k+1 right singular vectors. Its rank is at most k, so rank–nullity gives a unit vector z in this subspace with Bz=0. Then ‖(A−B)z‖=‖Az‖≥σ_{k+1}, since all its included strengths are at least that large. The truncated matrix attains equality. This proves the second norm claim. With repeated values at the truncation boundary, several retained subspaces can be optimal; best error need not identify a unique approximation.

Let a data matrix X contain m records as rows and d features as columns. Compute the feature mean μ using training rows and form Z=X−1μᵀ, so each column sums to zero. Its right singular vectors are principal feature directions. For k retained columns V_k, PCA scores are T=ZV_k, shape m×k. Reconstruction is X_hat=TV_kᵀ+1μᵀ. Subtraction and addition of the mean make this an affine reconstruction of the original records, with a linear projection of their centred displacements.

The sample scatter matrix is ZᵀZ. Its eigendirections are V and its eigenvalues are σᵢ². With m>1, dividing by m−1 gives the customary sample covariance matrix and eigenvalues σᵢ²/(m−1). At this stage this is an algebraic normalisation convention; distributional estimation claims come later. Scores along orthonormal principal directions have diagonal cross-product matrix diag(σ₁²,…,σ_k²), since TᵀT=V_kᵀZᵀZV_k. This explains the decorrelation of these centred coordinates without claiming statistical independence.

Explained variance ratios are σᵢ²/Σσⱼ². Their retained sum measures the fraction of centred training squared spread preserved by the projection. If all training records are identical, Z=0 and the denominator is zero; the ratios are undefined rather than automatically perfect. With one record, the sample-covariance denominator is also undefined. State these boundary cases instead of quietly emitting NaN or changing the mathematical definition.

PCA chooses an orthonormal k-dimensional feature subspace that minimises centred squared reconstruction error. For any candidate orthonormal W, its projection is ZWWᵀ, and retained energy is the sum of the selected direction Rayleigh values of ZᵀZ. The same projection-weight argument from the best-rank theorem bounds this by the top k squared strengths. Taking V_k attains the bound. Its scores reconstruct Z_k=U_kΣ_kV_kᵀ, exactly the truncated SVD of the centred data. Thus variance maximisation and squared-error minimisation are two algebraic views of the same objective.

Because centring forces 1ᵀZ=0, rank is at most min(d,m−1), not merely min(d,m). This matters for small training sets. A wide dataset cannot produce more than m−1 positive principal components after centring. Repeated scatter eigenvalues also make individual principal directions nonunique, while the corresponding whole eigenspace remains well defined. A sign flip changes scores and basis signs together and leaves reconstruction unchanged.

For a new mathematical column x, use t=V_kᵀ(x−μ), with t of length k. For a stored row, compute (x−μ)V_k. Reuse the training μ and V_k. The new row need not project with zero residual, and held-out scores need not have mean zero; neither indicates failure. They are measurements in a coordinate system fixed by training data. Refitting that system on each test batch creates a different transformation and changes what the evaluation means.

For one new centred displacement z, the reconstruction is V_kV_kᵀz. Its discarded part is perpendicular to the retained subspace, so its error square is the sum of its discarded coordinate squares. Training tail strengths give the aggregate training error, but they do not predict this new displacement’s error automatically. A new record can lie largely outside the retained subspace even when training reconstruction was excellent. Evaluate held-out reconstruction separately if that is the intended task.

PCA depends on feature units because its Euclidean reconstruction objective adds coordinate squares. Multiplying one feature by a hundred multiplies its contribution to squared spread by ten thousand and alters the scatter cross terms. A largest-spread direction can therefore change even when the underlying physical measurements are the same. In the lab, changing the second feature’s scale rotates the leading direction almost onto that axis and changes its explained-energy share. Those are meaningful consequences of a changed coordinate metric, not automatically a computational error.

Standardising features by training means and scales can create another justified objective, but it is a modelling choice. A zero-variance feature cannot be divided by its zero standard deviation; remove it or define a documented policy. Scales and any feature-selection decisions belong to the fitted transformation and must also be estimated from training data. A centre, scale and basis have to be kept together for subsequent scores and inverse reconstruction in original units.

Leakage occurs when information reserved for evaluation influences the fitted model or preprocessing. PCA uses no labels, yet fitting its mean or directions on combined training and held-out features still changes the learning procedure. The usual held-out question concerns a transformation trained without those records. In Lab C, adding two distant held-out rows changes the mean from (13,0.5) to about(58.666667,5.333333). Replacing the training centre with that value is already a different transformation, even before its directions are recomputed.

The operational contract is simple: split records first; fit mean, any scales and components on training only; apply those saved parameters to validation and test records; use validation to choose k or other modelling choices; evaluate the final selected procedure on the reserved test data. In cross-validation, each fold needs its own training-fitted transform. Reusing components fitted on the whole dataset across folds gives each fold access to its held-out features. If a deliberately transductive setting permits using unlabelled evaluation features, state that different task explicitly rather than presenting its result as ordinary inductive evaluation.

Uncentred SVD is also a different objective. It approximates raw vectors about the origin, combining mean offset and variation. With records far from zero, a leading direction can primarily represent their common location rather than deviations around their mean. Uncentred factorisation can be useful for a stated application, but should not be silently interpreted as centred PCA. The lab compares line projectors to distinguish actual directional changes from arbitrary vector-sign changes.

Low-rank factors can support image approximation, embedding compression and compact parameter updates. The best-rank theorem guarantees the closest matrix under its specified norms, not the best downstream model, fastest implementation or minimum file size. An update UVᵀ has rank at most k when both factors have k columns, but choosing such an update for a learning task introduces optimisation and statistical questions beyond SVD truncation. Keep the algebraic rank statement separate from a claim of task equivalence.

Finally, numerical rank and principal-component interpretation need a declared resolution. Repeated or near-repeated strengths can make individual directions unstable even when their combined subspace and reconstruction remain stable. A weak inverse direction can magnify noise, while discarding it changes the effective operator. Record selected k, retained strengths, training centre/scales, dtype and tolerances alongside reconstruction error. Those details make the approximation reproducible and its limits assessable without treating one successful decomposition as evidence for every intended use.

Use the Python 3.11 and NumPy baseline in array preparation. Outputs below come from executing the downloadable scripts. Predict shapes, retained directions and error budgets first. Lab values can differ in harmless signed basis conventions; compare reconstructed maps and subspaces.

Lab A Generated-matrix compression

The script constructs an 8×6 matrix from three orthonormal direction pairs and strengths six, three and one. Verify those generating directions, then compare reduced dimensions with exact rank. For each retained rank, predict the discarded squared error and factor scalar count before running.

Download lab1_truncation.py

"""Generated low-rank matrix with known spectrum and measured truncation error."""
import numpy as np
left = np.column_stack((np.ones(8), [1,-1,1,-1,1,-1,1,-1], [1,1,-1,-1,1,1,-1,-1]))/np.sqrt(8)
right = np.column_stack((np.ones(6)/np.sqrt(6), np.array([1,-1,1,-1,1,-1])/np.sqrt(6), np.array([1,1,-2,1,1,-2])/np.sqrt(12)))
assert np.allclose(left.T@left, np.eye(3)) and np.allclose(right.T@right, np.eye(3))
A = (left*np.array([6., 3., 1.]))@right.T
U, s, Vh = np.linalg.svd(A, full_matrices=False)
print("A / reduced U / s / Vh shapes:", A.shape, U.shape, s.shape, Vh.shape)
print("Singular values:", [round(float(v), 9) for v in s])
print("Original scalar entries:", A.size)
assert np.allclose(s[:3], [6,3,1]) and np.allclose(s[3:], 0)
for k in range(4):
    approximation = (U[:, :k]*s[:k])@Vh[:k, :]
    error_squared = np.linalg.norm(A-approximation, "fro")**2
    tail_squared = s[k:]@s[k:]
    scalar_storage = k*(A.shape[0]+A.shape[1]+1)
    print("Rank / squared Frobenius error / discarded sum / factor scalars:", k, round(float(error_squared), 9), round(float(tail_squared), 9), scalar_storage)
    assert np.allclose(error_squared, tail_squared)
assert np.allclose((U*s)@Vh, A)
print("Scalar-count savings are distinct from file-format overhead or predictive usefulness.")

A / reduced U / s / Vh shapes: (8, 6) (8, 6) (6,) (6, 6)
Singular values: [6.0, 3.0, 1.0, 0.0, 0.0, 0.0]
Original scalar entries: 48
Rank / squared Frobenius error / discarded sum / factor scalars: 0 46.0 46.0 0
Rank / squared Frobenius error / discarded sum / factor scalars: 1 10.0 10.0 15
Rank / squared Frobenius error / discarded sum / factor scalars: 2 1.0 1.0 30
Rank / squared Frobenius error / discarded sum / factor scalars: 3 0.0 0.0 45
Scalar-count savings are distinct from file-format overhead or predictive usefulness.

Deliverable: rank, factor shapes, reconstruction error and storage table, with a proof of the tail formula. Explain why a printed zero after rounding does not itself establish exact rank, and why the known construction supplies that exact claim here.

Lab B Training-centred PCA

Compute the training mean and scatter by hand. Compare their spectrum with squared singular values, then reconstruct one-component training data. Apply the saved training parameters to two held-out batches without refitting. Check score/reconstruction shapes and the aggregate tail-error identity.

Download lab2_training_pca.py

"""Fit mean and principal directions on training records; reuse on held-out rows."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)
training = np.array([[10.,0.], [12.,1.], [14.,0.], [16.,1.]])
test = np.array([[20.,2.], [8.,-1.]])
mean = training.mean(axis=0)
centered = training-mean
U, s, Vh = np.linalg.svd(centered, full_matrices=False)
directions = Vh.T
scores = centered@directions[:, :1]
reconstruction = scores@directions[:, :1].T+mean
test_scores = (test-mean)@directions[:, :1]
scatter = centered.T@centered
variances = s**2/(len(training)-1)
ratio = s**2/np.sum(s**2)
assert np.allclose(centered.mean(axis=0), 0)
assert np.allclose(scatter@directions, directions*s**2)
print("Training mean:", mean)
print("Principal directions in columns:", directions.tolist())
print("Sample covariance eigenvalues:", variances)
print("Explained variance ratios:", ratio)
print("Training rank-one reconstruction:", reconstruction.tolist())
print("Held-out scores using training transform:", test_scores.ravel())
error_squared = np.linalg.norm(training-reconstruction, "fro")**2
print("Rank-one reconstruction squared error:", error_squared)
assert np.allclose(error_squared, s[1]**2)
changed_test = np.array([[2000.,200.], [-800.,-100.]])
changed_scores = (changed_test-mean)@directions[:, :1]
print("Changed held-out scores using the same training transform:", changed_scores.ravel())
print("Mean reused for both held-out batches:", mean)
print("All-zero centred data variance ratio: undefined; guard division by zero.")

Training mean: [13.   0.5]
Principal directions in columns: [[0.9946229901291114, 0.10356209493161814], [0.10356209493161814, -0.9946229901291114]]
Sample covariance eigenvalues: [6.736081 0.263919]
Explained variance ratios: [0.962297 0.037703]
Training rank-one reconstruction: [[9.980672702257417, 0.18562172467192445], [12.062227727769086, 0.4023573132283918], [13.937772272230914, 0.5976426867716081], [16.019327297742585, 0.8143782753280755]]
Held-out scores using training transform: [ 7.117704 -5.128458]
Rank-one reconstruction squared error: 0.7917560805262006
Changed held-out scores using the same training transform: [1996.976519 -819.036482]
Mean reused for both held-out batches: [13.   0.5]
All-zero centred data variance ratio: undefined; guard division by zero.

Deliverable: saved mean and component directions, training scores/reconstruction, held-out scores and explained-energy ratios. State the undefined ratio policy for all-zero centred data and the difference between training and new-record error.

Lab C Centring leakage scaling and weak inverse modes

Compare centred and uncentred line projectors, then compute the leaked combined mean. Change feature units and explain the rotated principal line. Finally compare two relative pseudoinverse cutoffs under a small target perturbation.

Download lab3_centering_and_pinv_faults.py

"""Uncentred directions, evaluation leakage, feature units and weak inverse modes."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)
training = np.array([[10.,0.], [12.,1.], [14.,0.], [16.,1.]])
test = np.array([[100.,10.], [200.,20.]])
mean = training.mean(axis=0)
centered = training-mean
_, _, Vh = np.linalg.svd(centered, full_matrices=False)
_, _, uncentered_Vh = np.linalg.svd(training, full_matrices=False)
leaked_mean = np.vstack((training, test)).mean(axis=0)
print("Correct training mean / leaked combined mean:", mean, leaked_mean)
print("Centred / uncentred leading line projectors:", np.outer(Vh[0], Vh[0]).tolist(), np.outer(uncentered_Vh[0], uncentered_Vh[0]).tolist())
assert not np.allclose(mean, leaked_mean)
assert not np.allclose(np.outer(Vh[0], Vh[0]), np.outer(uncentered_Vh[0], uncentered_Vh[0]))
rescaled = centered*np.array([1.,100.])
_, scaled_s, scaled_Vh = np.linalg.svd(rescaled, full_matrices=False)
print("After feature two multiplied by 100, leading projector:", np.outer(scaled_Vh[0], scaled_Vh[0]).tolist())
print("Scaled leading explained ratio:", scaled_s[0]**2/(scaled_s@scaled_s))

A = np.diag([1., 1e-10])
target = np.array([1., 1e-10])
perturbed = target+np.array([0.,1e-8])
for cutoff in [1e-12,1e-8]:
    inverse = np.linalg.pinv(A, rcond=cutoff)
    base = inverse@target
    changed = inverse@perturbed
    print("Relative cutoff / base / changed weights / new residual:", cutoff, base, changed, format(np.linalg.norm(perturbed-A@changed), ".6e"))
assert np.allclose(np.linalg.pinv(A, rcond=1e-12)@perturbed, [1,101])
assert np.allclose(np.linalg.pinv(A, rcond=1e-8)@perturbed, [1,0])
print("Dropping a weak singular direction trades resolution and residual for reduced amplification; it changes the effective inverse.")

Correct training mean / leaked combined mean: [13.   0.5] [58.666667  5.333333]
Centred / uncentred leading line projectors: [[0.9892748924933744, 0.10300524052492092], [0.10300524052492092, 0.010725107506625487]] [[0.998380120869607, 0.04021511061772047], [0.04021511061772047, 0.0016198791303925755]]
After feature two multiplied by 100, leading projector: [[0.000401121600379904, 0.02002400314226935], [0.02002400314226935, 0.9995988783996206]]
Scaled leading explained ratio: 0.998403833102307
Relative cutoff / base / changed weights / new residual: 1e-12 [1. 1.] [  1. 101.] 1.654361e-24
Relative cutoff / base / changed weights / new residual: 1e-08 [1. 0.] [1. 0.] 1.010000e-08
Dropping a weak singular direction trades resolution and residual for reduced amplification; it changes the effective inverse.

Deliverable: identify each changed objective or evaluation procedure, and report the inverse resolution trade-off with actual residuals. Distinguish deliberate weak-mode truncation from a claim that the exact nonzero direction disappeared.

Exercises 1–12 are required; 13–14 add 35 optional minutes. Use real factors and state rank, retained dimension and centring assumptions.

Use SVD, pseudoinverse and PCA selections in the authors’ Mathematics for Machine Learning companion. The official NumPy 1.26 SVD interface and pseudoinverse interface specify the numerical shapes and cutoff conventions. The original derivations above prove the geometry and approximation claims.

Without notes, derive positive singular directions from AᵀA, state full/reduced/compact shapes, and recover the four spaces. Prove one optimal rank-k error bound, write the pseudoinverse, and specify a training-only PCA pipeline with a label-information counterexample.

Exit task: diag(4,2,0) has rank two. Rank-one squared error is four and operator error two; exact pseudoinverse is diag(1/4,1/2,0). A target’s third coordinate remains residual. Explain why its reduced SVD still has three entries, why zero must not be inverted, and what a cutoff discarding two would change.

Ready to move on: you can distinguish exact map geometry, numerical resolution, reconstruction and task evaluation. The next block develops limits, derivatives and optimisation, using this linear-algebra foundation. Check the course overview for availability.