Math for Data Science · MDS-2
Vectors, matrices, and the geometry that powers machine learning — built up to SVD and PCA.
Linear algebra is the language of data: every dataset is a matrix, every model a transformation. This course builds intuition and computation side by side — from dot products and elimination to eigenvalues, least squares, and the SVD/PCA techniques that drive dimensionality reduction in real machine learning pipelines. You will compute everything by hand and verify it with NumPy.
Course Outline
Every lesson is self-contained: a full explanation, worked step-by-step examples (in Python and R/Julia), common mistakes, a try-it-yourself, and an interactive quiz. Jump to a unit:
A vector is one object seen two ways: an arrow in space and a row of data. Both views drive everything in machine learning.
A vector is an ordered list of numbers, written like v = [3, 4]. Geometrically it is an arrow from the origin to the point (3, 4), carrying a direction and a length. In data science the same list is a record: a customer might be [age, income, visits] = [34, 52000, 7]. The dimension is just how many numbers it holds. This dual view is the heart of linear algebra for ML: every data point becomes a point in high-dimensional space, so geometric ideas like distance and angle become tools for comparing records. NumPy stores a vector as a 1-D array. Treating rows of a dataset as vectors lets us measure similarity, cluster, and project — the operations behind recommendation engines and embeddings.
Worked Example 1
Problem. Write the data row 'age 30, height 175cm' as a vector and state its dimension.
Answer. v = [30, 175], a 2-dimensional vector.
Worked Example 2
Problem. Plot the vector v = [3, 4] geometrically and identify the point it reaches.
Answer. An arrow from (0,0) to the point (3, 4).
Worked Example 3
Problem. A grayscale 2x2 image has pixels 10, 20, 30, 40. Represent it as a vector.
Answer. [10, 20, 30, 40], a vector in R^4.
Problem. A sensor reads temperature 22, humidity 60, pressure 1013. Write it as a vector and give its dimension.
Solution. Order the readings: [22, 60, 1013]. Counting three entries, it is a 3-dimensional vector in R^3. In NumPy: np.array([22, 60, 1013]).shape -> (3,).
Adding vectors combines movements; scaling stretches them. Together they generate every linear combination.
Two operations define a vector space. Addition is element-wise: [1,2] + [3,4] = [4,6], geometrically placing one arrow at the tip of the other (the parallelogram rule). Scaling multiplies every entry by a scalar: 3·[1,2] = [3,6], stretching the arrow without changing its line. Combining these gives a linear combination: a·v + b·w for scalars a, b. This single idea underlies nearly all of ML: a neural-network layer forms linear combinations of inputs, and a regression prediction is a linear combination of features weighted by coefficients. The set of all linear combinations of some vectors is their span — the space you can reach. Understanding linear combinations is the prerequisite for bases, rank, and dimensionality.
Worked Example 1
Problem. Compute [2, 5] + [3, -1].
Answer. [5, 4]
Worked Example 2
Problem. Scale v = [4, -2] by -1.5.
Answer. [-6, 3]
Worked Example 3
Problem. Form the linear combination 2·[1,0] + 3·[0,1] - 1·[1,1].
Answer. [1, 2]
Problem. Compute 3·[1,2] - 2·[2,1].
Solution. Scale: 3·[1,2]=[3,6] and 2·[2,1]=[4,2]. Subtract: [3-4, 6-2] = [-1, 4]. NumPy: 3*np.array([1,2]) - 2*np.array([2,1]) -> array([-1, 4]).
The dot product turns two vectors into a single number measuring how much they point the same way.
The dot product of u and v is the sum of element-wise products: u·v = u1·v1 + u2·v2 + .... It connects algebra to geometry through u·v = ||u|| ||v|| cos(θ), where θ is the angle between them. So a large positive dot product means the vectors align, zero means they are perpendicular, and negative means they oppose. In ML the dot product is everywhere: it scores how well an input matches a weight vector in a neuron, computes attention between query and key vectors in transformers, and measures document similarity. It is also the building block of matrix multiplication — every entry of a matrix product is a dot product of a row and a column. Mastering it makes the rest of linear algebra click.
Worked Example 1
Problem. Compute [1, 2, 3]·[4, 5, 6].
Answer. 32
Worked Example 2
Problem. Show [1, 0] and [0, 1] are perpendicular using the dot product.
Answer. Dot product is 0, so they are orthogonal.
Worked Example 3
Problem. Use the dot product to find the angle between u=[1,1] and v=[1,0].
Answer. 45 degrees
Problem. Compute [2, -1, 3]·[1, 4, 0] and state whether the angle is acute, right, or obtuse.
Solution. Dot = 2·1 + (-1)·4 + 3·0 = 2 - 4 + 0 = -2. Since the dot product is negative, the angle is obtuse (greater than 90 degrees). NumPy: np.dot([2,-1,3],[1,4,0]) -> -2.
A norm measures a vector's size; L1 and L2 are the two most used, and dividing by a norm gives a unit vector.
A norm assigns a length to a vector. The L2 (Euclidean) norm is sqrt of the sum of squares: ||v||_2 = sqrt(v1^2 + v2^2 + ...), the straight-line distance from the origin. The L1 norm sums absolute values: ||v||_1 = |v1| + |v2| + ..., the 'taxicab' distance. These choices matter in ML: L2 regularization (ridge) shrinks weights smoothly, while L1 (lasso) drives some weights exactly to zero, producing sparse models. A unit vector is any vector divided by its norm, giving length 1 while preserving direction — used to express pure direction, normalize embeddings, and define cosine similarity. Knowing which norm to use is a real modeling decision, not just notation.
Worked Example 1
Problem. Find the L2 norm of v = [3, 4].
Answer. 5
Worked Example 2
Problem. Find the L1 norm of v = [3, -4, 1].
Answer. 8
Worked Example 3
Problem. Turn v = [3, 4] into a unit vector.
Answer. [0.6, 0.8]
Problem. Find the L2 norm of [1, 2, 2] and then its unit vector.
Solution. L2 norm = sqrt(1+4+4) = sqrt(9) = 3. Unit vector = [1/3, 2/3, 2/3] ≈ [0.333, 0.667, 0.667]. NumPy: v=np.array([1,2,2]); v/np.linalg.norm(v) -> array([0.333, 0.667, 0.667]).
Cosine similarity uses the dot product and norms to compare direction, ignoring magnitude — the standard for comparing data vectors.
Rearranging u·v = ||u|| ||v|| cos(θ) gives cosine similarity: cos(θ) = (u·v) / (||u|| ||v||). It ranges from -1 (opposite) through 0 (orthogonal) to 1 (identical direction). Because it divides out magnitude, it compares only direction — perfect for text and embeddings where a long document and a short one on the same topic should score as similar. Orthogonality (cos = 0, dot = 0) means the vectors carry independent information, which is why orthogonal bases and uncorrelated features are prized. Cosine similarity powers recommendation systems, semantic search, and clustering. The key insight: in high dimensions, direction often matters more than length, and cosine similarity captures exactly that.
Worked Example 1
Problem. Compute the cosine similarity of u=[1,0] and v=[1,1].
Answer. ≈ 0.707 (45 degrees apart)
Worked Example 2
Problem. Show that scaling does not change cosine similarity: compare u=[1,0] with v=[1,1] and with 5v=[5,5].
Answer. Both give ≈ 0.707; cosine ignores magnitude.
Worked Example 3
Problem. Are document vectors d1=[2,0,1] and d2=[0,3,0] similar? Use cosine similarity.
Answer. Cosine similarity is 0 — the documents are orthogonal (share no words).
Problem. Compute the cosine similarity of u=[3,4] and v=[4,3].
Solution. Dot = 3·4 + 4·3 = 24. Norms: ||u||=5, ||v||=5. cos = 24/(5·5) = 24/25 = 0.96. The vectors point in very similar directions. NumPy: np.dot([3,4],[4,3])/(5*5) -> 0.96.
NumPy arrays are how vectors live in code; shape and broadcasting rules govern every operation.
NumPy represents a vector as a 1-D ndarray created with np.array([...]). Its shape attribute reports dimensions: a length-3 vector has shape (3,). Broadcasting is NumPy's rule for combining arrays of different but compatible shapes without explicit loops — e.g. adding a scalar to a vector applies it to every entry, and adding a (3,) vector to a (2,3) matrix adds it to each row. This vectorized style is both concise and fast because it runs in optimized C, which is essential when datasets have millions of rows. Understanding shape prevents the most common bugs: a (3,) vector behaves differently from a (3,1) column or (1,3) row. Checking .shape before operating is a habit that saves hours of debugging in real pipelines.
Worked Example 1
Problem. Create the vector [10, 20, 30] in NumPy and report its shape.
Answer. shape = (3,)
Worked Example 2
Problem. Use broadcasting to add 5 to every entry of [1, 2, 3].
Answer. [6, 7, 8]
Worked Example 3
Problem. Add the row vector [1, 0, -1] to each row of the matrix [[1,1,1],[2,2,2]] via broadcasting.
Answer. [[2, 1, 0], [3, 2, 1]]
Problem. Given m = np.array([[10,20],[30,40]]) and v = np.array([1,2]), what is m + v?
Solution. v broadcasts across each row: row 1 [10,20]+[1,2]=[11,22]; row 2 [30,40]+[1,2]=[31,42]. Result: array([[11,22],[31,42]]). The (2,) vector aligns with the last axis (length 2) of the (2,2) matrix.
Represent three short documents as word-count vectors over a small shared vocabulary. Compute each pairwise dot product, L2 norm, and cosine similarity by hand for one pair, then reproduce all pairs with NumPy using np.dot and np.linalg.norm. Identify which two documents are most similar and explain why cosine similarity is preferred over raw dot product here.
Deliverable · A notebook showing the hand calculation, the NumPy code, the cosine-similarity matrix, and a short written conclusion.
1. Two vectors are orthogonal when their dot product equals:
Answer C. Orthogonal vectors meet at 90 degrees, and cos(90°) = 0, so their dot product is 0.
2. The L2 norm of the vector [3, 4] is:
Answer B. The L2 norm is sqrt(3^2 + 4^2) = sqrt(9 + 16) = sqrt(25) = 5.
3. Cosine similarity differs from the raw dot product because it:
Answer B. Cosine similarity divides by both norms, so it measures direction alignment independent of magnitude.
4. A linear combination of vectors v and w is any expression of the form:
Answer B. A linear combination scales each vector by a scalar and adds the results.
5. To turn a nonzero vector into a unit vector you:
Answer C. Dividing a vector by its own norm yields a vector pointing the same way with length 1.
I can compute the dot product, length, and angle between two vectors by hand and in NumPy.
I can use cosine similarity to measure how alike two data points are.
I can explain when two vectors are orthogonal and why it matters.
A matrix is two things at once: a grid of data and a machine that transforms vectors. Both views matter in ML.
A matrix is a rectangular array of numbers with m rows and n columns, written as an m×n matrix. The data view treats each row as a record and each column as a feature, so an entire dataset is one matrix X. The transformation view treats a matrix A as a function that takes a vector v and returns Av, a new vector — it rotates, scales, shears, or projects space in a linear way (straight lines stay straight, the origin stays fixed). These two views unify ML: a layer of a neural network multiplies inputs by a weight matrix, simultaneously a data operation and a geometric transformation. Recognizing that the same grid of numbers is both a table and an operator is the conceptual leap of this unit.
Worked Example 1
Problem. A dataset has 3 customers, each with [age, spend]. Write it as a matrix and give its shape.
Answer. A 3×2 matrix.
Worked Example 2
Problem. Apply the scaling matrix A = [[2,0],[0,3]] to the vector v = [1, 1].
Answer. [2, 3]
Worked Example 3
Problem. Apply the 90-degree rotation matrix R = [[0,-1],[1,0]] to v = [1, 0] and explain the result.
Answer. [0, 1] — a 90-degree counterclockwise rotation.
Problem. Apply A = [[1,1],[0,1]] (a shear) to v = [2, 3].
Solution. Row 1 dot v: 1·2 + 1·3 = 5. Row 2 dot v: 0·2 + 1·3 = 3. So Av = [5, 3]. The x-coordinate is sheared by adding y, while y is unchanged. NumPy: np.array([[1,1],[0,1]]) @ np.array([2,3]) -> array([5, 3]).
The simplest matrix operations are element-wise; the transpose flips a matrix across its diagonal and appears everywhere in ML formulas.
Matrix addition is element-wise and requires identical shapes: (A+B)[i][j] = A[i][j] + B[i][j]. Scalar multiplication scales every entry: (cA)[i][j] = c·A[i][j]. The transpose A^T swaps rows and columns, so an m×n matrix becomes n×m and entry (i,j) moves to (j,i). The transpose is central in data science: the Gram matrix A^T A summarizes feature correlations and appears in least squares, covariance is computed from centered-data transposes, and many identities like (AB)^T = B^T A^T rely on it. A symmetric matrix satisfies A = A^T. These operations are cheap but foundational — every higher technique (normal equations, covariance, kernels) is written using transposes and sums.
Worked Example 1
Problem. Add A = [[1,2],[3,4]] and B = [[5,6],[7,8]].
Answer. [[6, 8], [10, 12]]
Worked Example 2
Problem. Compute the transpose of A = [[1,2,3],[4,5,6]].
Answer. [[1, 4], [2, 5], [3, 6]]
Worked Example 3
Problem. Show A = [[2,1],[1,3]] is symmetric, and scale it by 0.5.
Answer. A is symmetric; 0.5A = [[1, 0.5], [0.5, 1.5]].
Problem. Given A = [[1,0,2],[3,1,0]], compute A^T and state its shape.
Solution. A is 2×3 so A^T is 3×2. Swap rows and columns: A^T = [[1,3],[0,1],[2,0]]. NumPy: np.array([[1,0,2],[3,1,0]]).T -> array([[1,3],[0,1],[2,0]]), shape (3,2).
A matrix times a vector is best understood as a weighted sum of the matrix's columns — the key intuition behind span and column space.
The product Av can be computed row-by-row (each output entry is the dot product of a row of A with v), but the deeper view is column-based: Av is a linear combination of A's columns, where the entries of v are the weights. If A has columns c1, c2 and v = [a, b], then Av = a·c1 + b·c2. This reframes solving Ax = b as asking which combination of columns produces b — directly motivating column space, rank, and span. In ML this is exactly how a linear model forms predictions: features (columns) are weighted by learned coefficients (the vector) and summed. Seeing Av as mixing columns, not just crunching rows, unlocks the geometry of linear systems.
Worked Example 1
Problem. Compute A v with A = [[1,2],[3,4]] and v = [5, 6] using the row view.
Answer. [17, 39]
Worked Example 2
Problem. Recompute the same product as a linear combination of columns.
Answer. [17, 39] (5·c1 + 6·c2)
Worked Example 3
Problem. For A = [[1,0,2],[0,1,1]] and v = [3,4,1], compute Av and name its shape.
Answer. [5, 5], a 2-dimensional vector.
Problem. Compute Av for A = [[2,1],[0,3],[1,1]] and v = [1, 2].
Solution. A is 3×2, v length 2, so Av is length 3. Row 1: 2·1+1·2=4; row 2: 0·1+3·2=6; row 3: 1·1+1·2=3. Av = [4, 6, 3]. As columns: 1·[2,0,1] + 2·[1,3,1] = [4,6,3]. NumPy: np.array([[2,1],[0,3],[1,1]]) @ np.array([1,2]) -> array([4,6,3]).
Multiplying matrices composes transformations; each output entry is a row-times-column dot product, and AB usually differs from BA.
To multiply A (m×k) by B (k×n), the inner dimensions must match; the result is m×n. Entry (i,j) of AB is the dot product of row i of A with column j of B. Conceptually, AB is the composition of two transformations: first apply B, then A, to a vector — (AB)v = A(Bv). Because order of operations matters (rotating then scaling differs from scaling then rotating), matrix multiplication is generally non-commutative: AB ≠ BA. This composition view explains why neural networks stack weight matrices and why the chain of transformations in a pipeline must be ordered correctly. Mastering shape rules and the row-column recipe is essential before determinants, inverses, and decompositions.
Worked Example 1
Problem. Multiply A = [[1,2],[3,4]] by B = [[5,6],[7,8]].
Answer. [[19, 22], [43, 50]]
Worked Example 2
Problem. Show AB ≠ BA using the matrices above.
Answer. BA = [[23,34],[31,46]] ≠ AB, so multiplication is non-commutative.
Worked Example 3
Problem. Multiply A (2×3) = [[1,0,1],[2,1,0]] by B (3×2) = [[1,2],[0,1],[3,0]] and give the shape.
Answer. [[4, 2], [2, 5]], a 2×2 matrix.
Problem. Compute AB for A = [[2,0],[1,3]] and B = [[1,4],[2,1]].
Solution. (1,1)=2·1+0·2=2; (1,2)=2·4+0·1=8; (2,1)=1·1+3·2=7; (2,2)=1·4+3·1=7. AB = [[2,8],[7,7]]. NumPy: np.array([[2,0],[1,3]]) @ np.array([[1,4],[2,1]]) -> array([[2,8],[7,7]]).
A few special matrix shapes appear constantly: the identity does nothing, diagonals scale axes independently, and symmetric matrices have beautiful eigen-structure.
The identity matrix I has 1s on the main diagonal and 0s elsewhere; it is the multiplicative neutral element, so IA = AI = A and Iv = v. A diagonal matrix has nonzero entries only on the main diagonal; multiplying by it scales each coordinate independently, which makes powers and inverses trivial (just power or reciprocate each diagonal entry). A symmetric matrix satisfies A = A^T; these are pivotal in data science because covariance matrices and Gram matrices A^T A are always symmetric, and the spectral theorem guarantees them real eigenvalues and orthogonal eigenvectors — the foundation of PCA. Recognizing these structures lets you shortcut computation and predict an algorithm's behavior.
Worked Example 1
Problem. Verify that the 2×2 identity leaves v = [7, -2] unchanged.
Answer. Iv = [7, -2], unchanged.
Worked Example 2
Problem. Multiply the diagonal D = diag(2, 5) by v = [3, 4].
Answer. [6, 20]
Worked Example 3
Problem. Show that for any A = [[1,2],[3,4]], the product A^T A is symmetric.
Answer. A^T A = [[10,14],[14,20]], which is symmetric.
Problem. Compute D^2 for the diagonal matrix D = diag(3, -2).
Solution. For diagonal matrices, squaring just squares each diagonal entry: D^2 = diag(9, 4) = [[9,0],[0,4]]. NumPy: D = np.diag([3,-2]); D @ D -> array([[9,0],[0,4]]).
In NumPy, @ does true matrix multiplication while * does element-wise — confusing them is one of the most common and silent bugs in data code.
NumPy stores a matrix as a 2-D ndarray; .shape reports (rows, columns). The @ operator (equivalently np.matmul or np.dot) performs true matrix multiplication using the row-column rule and respects inner-dimension matching. The * operator performs element-wise (Hadamard) multiplication and requires broadcast-compatible shapes — it will silently produce a wrong but valid-looking result if you meant matrix multiplication. This single distinction causes a large share of real bugs in ML pipelines. Always check shapes with .shape before multiplying, and use @ for linear-algebra products. Knowing this keeps transformations, weight multiplications, and Gram matrices correct, and it carries over to deep-learning libraries that follow the same convention.
Worked Example 1
Problem. Given A = [[1,2],[3,4]] and B = [[1,0],[0,1]], compute A @ B in NumPy.
Answer. [[1, 2], [3, 4]] (unchanged, since B is identity)
Worked Example 2
Problem. Contrast A @ B with A * B for A = [[1,2],[3,4]], B = [[5,6],[7,8]].
Answer. @ gives [[19,22],[43,50]]; * gives [[5,12],[21,32]].
Worked Example 3
Problem. Predict and verify the shape of A @ v where A is (2,3) and v has shape (3,).
Answer. Shape (2,); the product is [6, 7].
Problem. For A = [[2,0],[0,2]] and B = [[1,2],[3,4]], compute both A @ B and A * B.
Solution. A @ B (A doubles everything): [[2,4],[6,8]]. A * B element-wise: [[2·1,0·2],[0·3,2·4]] = [[2,0],[0,8]]. They differ. NumPy: A@B -> array([[2,4],[6,8]]); A*B -> array([[2,0],[0,8]]).
Construct rotation and scaling matrices for 2D space. By hand, apply your rotation matrix to a single point, then in NumPy apply it to a small set of points stored as columns of a matrix using the @ operator. Verify that (AB) applied to a vector equals A applied to (B applied to the vector), and show that A * B (element-wise) gives a different, wrong result.
Deliverable · A notebook with the hand calculation, the transformed point set plotted, and a short note on why @ and * differ.
1. If A is 3x2 and B is 2x4, the product AB has shape:
Answer A. Inner dimensions (2 and 2) match, and the result takes the outer dimensions: 3x4.
2. Matrix multiplication is generally:
Answer B. In general AB ≠ BA; order matters because transformations compose in sequence.
3. Multiplying any vector by the identity matrix I gives:
Answer C. The identity matrix is the multiplicative neutral element: Iv = v.
4. The transpose of a matrix product (AB)^T equals:
Answer B. Transposing a product reverses the order: (AB)^T = B^T A^T.
5. In NumPy, the correct operator for true matrix multiplication is:
Answer C. The @ operator (or np.matmul) does matrix multiplication; * does element-wise multiplication.
I can multiply matrices by hand and predict the shape of the result.
I can interpret a matrix-vector product as a linear transformation of space.
I can distinguish element-wise from true matrix multiplication in NumPy.
Every system of linear equations compresses into one matrix equation Ax = b, the gateway to solving with linear algebra.
A linear system is a set of equations like 2x + y = 5, x - y = 1, all sharing the unknowns. We package the coefficients into a matrix A, the unknowns into a vector x, and the right-hand sides into a vector b, giving the compact form Ax = b. Reading it through the column view, Ax is a linear combination of A's columns, so solving asks which weights x produce b. This single equation generalizes to thousands of variables and is how regression, circuit analysis, and equilibrium models are stated. The augmented matrix [A | b] keeps everything in one grid for elimination. Mastering this translation from words to Ax = b is the first and most important step in solving real systems.
Worked Example 1
Problem. Write the system 2x + y = 5, x - y = 1 in the form Ax = b.
Answer. A = [[2,1],[1,-1]], x = [x,y], b = [5,1].
Worked Example 2
Problem. Form the augmented matrix for x + 2y + z = 4, 3y - z = 2, 2x + y = 0.
Answer. [[1,2,1,4],[0,3,-1,2],[2,1,0,0]]
Worked Example 3
Problem. Check that x = [2, 1] solves Ax = b for A = [[2,1],[1,-1]], b = [5,1].
Answer. Yes; Ax = [5,1] = b.
Problem. Write 3x - 2y = 4, x + y = 7 as Ax = b and verify x = [2, 5].
Solution. A = [[3,-2],[1,1]], b = [4,7]. Check: 3·2 - 2·5 = 6-10 = -4, not 4. So x = [2,5] is NOT a solution. (The actual solution is x=18/5, y=17/5.) NumPy: np.array([[3,-2],[1,1]]) @ np.array([2,5]) -> array([-4, 7]).
Gaussian elimination uses three reversible row operations to simplify a system into an easily solved triangular form.
Gaussian elimination systematically zeroes out entries below each pivot to reach an upper-triangular (echelon) form, from which back-substitution recovers the solution. It uses exactly three elementary row operations, each preserving the solution set: swap two rows, multiply a row by a nonzero scalar, and add a multiple of one row to another. Applied to the augmented matrix [A | b], these operations mirror the legal algebraic manipulations of the equations. This is the workhorse algorithm behind solving systems, computing rank, and finding inverses, and it is what np.linalg.solve does internally (with pivoting for stability). Doing it by hand builds the intuition for pivots, rank, and the three solution cases that follow.
Worked Example 1
Problem. Eliminate to solve x + y = 3, 2x + 3y = 8.
Answer. x = 1, y = 2
Worked Example 2
Problem. Use a row swap to handle a zero pivot: solve 0x + y = 2, x + y = 5.
Answer. x = 3, y = 2
Worked Example 3
Problem. Forward-eliminate the 3×3 system with augmented matrix [[1,1,1,6],[0,2,1,8],[2,1,3,14]] (first row already pivoted).
Answer. x = 0, y = 2, z = 4
Problem. Solve 2x + 4y = 10, 3x + 7y = 16 by elimination.
Solution. Augmented [[2,4,10],[3,7,16]]. R1 -> R1/2 = [1,2,5]. R2 -> R2 - 3·R1 = [0,1,1]. So y = 1; x + 2·1 = 5 gives x = 3. Solution x=3, y=1. NumPy: np.linalg.solve([[2,4],[3,7]],[10,16]) -> array([3., 1.]).
Echelon forms are the standardized 'staircase' results of elimination; RREF is the unique, fully simplified version you can read solutions off directly.
A matrix is in row echelon form (REF) when each leading nonzero entry (pivot) sits to the right of the pivot above it and all-zero rows are at the bottom — a downward staircase. Reduced row echelon form (RREF) goes further: every pivot is 1 and is the only nonzero entry in its column. RREF is unique for a given matrix, so it is the canonical fingerprint used to read off solutions, rank, and free variables without ambiguity. Where REF needs back-substitution, RREF gives the solution almost for free. SymPy's rref() produces it exactly with rational arithmetic, avoiding floating-point error. Understanding the staircase structure makes pivots, rank, and the solution cases immediately visible.
Worked Example 1
Problem. Is [[1,2,3],[0,1,4],[0,0,1]] in row echelon form?
Answer. Yes, it is in row echelon form.
Worked Example 2
Problem. Reduce [[1,2,5],[0,1,2]] (REF) to RREF.
Answer. [[1,0,1],[0,1,2]] (x = 1, y = 2)
Worked Example 3
Problem. Find the RREF of A = [[1,2,1],[2,4,3]] and identify the pivot columns.
Answer. RREF = [[1,2,0],[0,0,1]]; pivot columns are 1 and 3.
Problem. Reduce [[2,4,6],[1,1,2]] to RREF.
Solution. R1 -> R1/2 = [1,2,3]. R2 -> R2 - R1 = [0,-1,-1]. R2 -> -R2 = [0,1,1]. R1 -> R1 - 2·R2 = [1,0,1]. RREF = [[1,0,1],[0,1,1]]. SymPy: Matrix([[2,4,6],[1,1,2]]).rref() -> ([[1,0,1],[0,1,1]], (0,1)).
Pivots count the genuinely independent equations; their number is the rank, and the variables without pivots are free, signaling infinitely many solutions.
After elimination, each pivot marks an independent constraint. The rank of a matrix is the number of pivots — equivalently the number of independent rows or columns. Columns that contain a pivot correspond to pivot (basic) variables; columns without a pivot correspond to free variables, which can take any value while the basic variables adjust to satisfy the equations. The count is governed by rank-nullity: (number of variables) = (pivot variables) + (free variables). Rank is one of the most informative numbers in data science: it tells you the true dimensionality of your features, exposes redundancy (collinear columns), and predicts whether a system or regression has a unique solution. Reading rank and free variables straight off the echelon form is a core skill.
Worked Example 1
Problem. Find the rank of A = [[1,2],[2,4]].
Answer. rank = 1
Worked Example 2
Problem. For RREF [[1,0,3],[0,1,2]] of a 2×3 system, identify pivot and free variables.
Answer. x1, x2 are basic; x3 is free (rank 2, 1 free variable).
Worked Example 3
Problem. Determine the rank of A = [[1,1,1],[2,2,2],[3,3,3]].
Answer. rank = 1
Problem. Find the rank of A = [[1,2,3],[0,1,4],[0,0,0]].
Solution. It is already in echelon form. Pivots sit in rows 1 and 2 (positions (1,1) and (2,2)); row 3 is all zeros. So there are 2 pivots and rank = 2. NumPy: np.linalg.matrix_rank([[1,2,3],[0,1,4],[0,0,0]]) -> 2.
Every linear system falls into exactly one of three outcomes, and the echelon form tells you which by comparing rank to variables and checking for contradictions.
A linear system Ax = b has exactly one of three fates. It is inconsistent (no solution) when elimination produces a row like [0 0 0 | c] with c ≠ 0, a literal contradiction 0 = c. It has a unique solution when it is consistent and every variable is a pivot variable (rank equals the number of unknowns, no free variables). It has infinitely many solutions when it is consistent but has at least one free variable (rank less than the number of unknowns). Compactly: compare rank(A) to rank([A|b]) — unequal means no solution; equal-and-full means one; equal-but-deficient means infinitely many. This classification underlies whether a model is identifiable and whether data uniquely determines parameters.
Worked Example 1
Problem. Classify x + y = 2, 2x + 2y = 5.
Answer. No solution (inconsistent).
Worked Example 2
Problem. Classify x + y = 3, x - y = 1.
Answer. Exactly one solution: x = 2, y = 1.
Worked Example 3
Problem. Classify x + 2y = 4, 2x + 4y = 8.
Answer. Infinitely many solutions: x = 4 - 2t, y = t.
Problem. Classify the system x + y + z = 1, 2x + 2y + 2z = 3.
Solution. Augmented [[1,1,1,1],[2,2,2,3]]. R2 -> R2 - 2·R1 = [0,0,0,1], i.e. 0 = 1, a contradiction. The system is inconsistent and has no solution (rank(A)=1 < rank([A|b])=2).
NumPy solves numerically and fast; SymPy solves symbolically and exactly — each suits different needs in practice.
np.linalg.solve(A, b) returns the unique solution of a square, full-rank system using fast LU decomposition with partial pivoting; it raises an error if A is singular. For non-square or rank-deficient systems, np.linalg.lstsq gives a least-squares answer. SymPy works symbolically: Matrix(A).rref() returns the exact reduced row echelon form with rational arithmetic, and linsolve handles parametrized solution sets including free variables, expressing infinitely many solutions explicitly. Use NumPy when you have numeric data and need speed; use SymPy when you need exact answers, want to see the RREF, or must describe a solution family. Together they let you both compute and understand systems, and they confirm by hand work without rounding error.
Worked Example 1
Problem. Solve 3x + y = 9, x + 2y = 8 with np.linalg.solve.
Answer. x = 2, y = 3
Worked Example 2
Problem. Use SymPy rref to solve the augmented system [[1,1,5],[2,-1,1]].
Answer. x = 2, y = 3
Worked Example 3
Problem. Detect a singular system: try np.linalg.solve on A = [[1,2],[2,4]], b = [3, 6].
Answer. solve raises a Singular matrix error; lstsq gives [0.6, 1.2].
Problem. Solve 4x - y = 1, x + y = 4 using np.linalg.solve, and verify by hand.
Solution. A = [[4,-1],[1,1]], b = [1,4]. Adding equations conceptually: from x+y=4, y=4-x; sub: 4x-(4-x)=1 -> 5x=5 -> x=1, y=3. NumPy: np.linalg.solve([[4,-1],[1,1]],[1,4]) -> array([1., 3.]). Check: 4·1 - 3 = 1, 1 + 3 = 4. Correct.
Write a 3-equation, 3-unknown system as an augmented matrix and perform Gaussian elimination by hand to RREF, recording each row operation. Then use SymPy's rref() to confirm your result and np.linalg.solve to find the solution. Determine the rank and state whether the system has zero, one, or infinitely many solutions, justifying with the pivot count.
Deliverable · A notebook showing hand elimination steps, the SymPy/NumPy verification, the rank, and the solution classification.
1. The rank of a matrix equals the number of:
Answer C. Rank is the count of pivots, which equals the number of independent rows (or columns).
2. A square system Ax = b has a unique solution exactly when A is:
Answer B. A full-rank square matrix is invertible, giving the single solution x = A^{-1}b.
3. If elimination produces a row like [0 0 0 | 5], the system is:
Answer B. That row says 0 = 5, a contradiction, so the system has no solution.
4. A free variable in the solution indicates:
Answer C. Each free variable can take any value, producing infinitely many solutions.
5. Which is NOT a valid elementary row operation?
Answer D. Deleting rows is not an elementary row operation; the three valid ones preserve the solution set.
I can solve a linear system by Gaussian elimination and read off the rank.
I can determine whether a system has zero, one, or infinitely many solutions.
I can solve and verify systems with np.linalg.solve and SymPy's rref.
The determinant is one number that says how much a matrix stretches space and whether it flips orientation — a geometric summary of a transformation.
The determinant of a square matrix A, written det(A), is the factor by which A scales area (in 2D) or volume (in higher dimensions). A determinant of 3 means the transformation triples areas; a determinant of 1 preserves area. The sign encodes orientation: a negative determinant means the transformation flips space (like a mirror). Most importantly, det(A) = 0 means A squashes space into a lower dimension — area collapses to zero — which signals that A is singular and not invertible. In data science the determinant flags collinearity and degeneracy (a zero determinant means redundant features), and it appears in the multivariate Gaussian density and change-of-variables formulas. Seeing it as a volume scaler makes its algebraic rules intuitive.
Worked Example 1
Problem. Interpret det([[2,0],[0,3]]).
Answer. det = 6 (areas multiplied by 6).
Worked Example 2
Problem. What does det([[1,2],[2,4]]) = 0 tell you geometrically?
Answer. det = 0: space collapses to a line; A is singular.
Worked Example 3
Problem. Interpret the negative determinant of the reflection [[0,1],[1,0]].
Answer. det = -1: area preserved but orientation reversed.
Problem. Without full computation, what is det([[5,0],[0,0]]) and what does it mean?
Solution. It is diagonal with entries 5 and 0; the product is 5·0 = 0. Geometrically the y-axis is collapsed to a point, so space flattens to a line and the matrix is singular. NumPy: np.linalg.det([[5,0],[0,0]]) -> 0.0.
Small determinants have direct formulas: the 2×2 cross-difference and the 3×3 cofactor expansion.
For a 2×2 matrix [[a,b],[c,d]], det = ad - bc, the main-diagonal product minus the anti-diagonal product. For a 3×3 matrix, cofactor expansion along the first row gives det = a(ei - fh) - b(di - fg) + c(dh - eg), where each term multiplies an entry by the 2×2 determinant of the submatrix left after deleting its row and column, with alternating signs (+, -, +). This recursive minor-expansion generalizes to any size, though it grows expensive; in practice large determinants are computed via LU decomposition (what NumPy uses). Mastering the 2×2 and 3×3 by hand builds the foundation for understanding invertibility, Cramer's rule, and eigenvalue characteristic polynomials, all of which lean on these formulas.
Worked Example 1
Problem. Compute det([[4,3],[6,3]]).
Answer. -6
Worked Example 2
Problem. Compute det([[1,2,3],[0,1,4],[5,6,0]]) by cofactor expansion along row 1.
Answer. 1
Worked Example 3
Problem. Compute det([[2,0,1],[1,3,2],[1,0,2]]) using the first column (it has two useful zeros via row 2,3? expand row with zeros).
Answer. 9
Problem. Compute det([[3,1,0],[2,4,1],[0,0,5]]).
Solution. Expand along the bottom row (two zeros): only the (3,3) entry contributes: +5·det([[3,1],[2,4]]) = 5·(3·4 - 1·2) = 5·10 = 50. NumPy: round(np.linalg.det([[3,1,0],[2,4,1],[0,0,5]])) -> 50.
A handful of determinant rules let you reason about transformations without grinding through the full computation.
Key properties: det(I) = 1; swapping two rows flips the sign; multiplying a row by a scalar c multiplies det by c; adding a multiple of one row to another leaves det unchanged (the operation used in elimination). The determinant is multiplicative: det(AB) = det(A)·det(B), and det(A^T) = det(A). For a triangular matrix, det is just the product of the diagonal entries — which is why elimination to triangular form computes determinants efficiently. The cornerstone fact: det(A) = 0 if and only if A is singular (not invertible, rank-deficient, columns linearly dependent). In data science a zero or near-zero determinant warns of collinear features and numerically unstable inverses, making these properties practical diagnostics.
Worked Example 1
Problem. Use the product rule to find det(AB) where det(A)=3 and det(B)=-2.
Answer. -6
Worked Example 2
Problem. Find the determinant of the triangular matrix [[2,7,1],[0,3,5],[0,0,4]].
Answer. 24
Worked Example 3
Problem. A row operation R2 -> R2 - 5·R1 is applied to A. How does det change, and what if we then swap two rows?
Answer. The elimination step keeps det the same; the swap negates it, giving -det(A).
Problem. If det(A) = 4 for a 3×3 matrix A, what is det(2A)?
Solution. Scaling all 3 rows by 2 multiplies the determinant by 2^3 = 8. So det(2A) = 8·det(A) = 8·4 = 32. NumPy: np.linalg.det(2*A) equals 8*np.linalg.det(A).
The inverse undoes a matrix transformation; it exists exactly when the matrix is square and has a nonzero determinant.
The inverse A^{-1} of a square matrix A is the unique matrix with A·A^{-1} = A^{-1}·A = I. Multiplying by A^{-1} reverses the transformation A performed, which is how Ax = b is solved as x = A^{-1}b (conceptually — in practice we solve directly for stability). An inverse exists if and only if A is square and det(A) ≠ 0 (equivalently, full rank, independent columns). For a 2×2 matrix there is a closed form: A^{-1} = (1/det)·[[d,-b],[-c,a]]. Inverses underlie the least-squares normal equations and covariance whitening, but in numerical practice we avoid forming them explicitly because they are slower and less stable than solving the system directly. Knowing when an inverse exists is essential before relying on it.
Worked Example 1
Problem. Find the inverse of A = [[4,7],[2,6]].
Answer. [[0.6, -0.7], [-0.2, 0.4]]
Worked Example 2
Problem. Verify the inverse from Example 1 satisfies A·A^{-1} = I.
Answer. A·A^{-1} = I, so the inverse is correct.
Worked Example 3
Problem. Does A = [[1,2],[2,4]] have an inverse?
Answer. No; det = 0 so A is not invertible.
Problem. Find the inverse of A = [[2,1],[1,1]].
Solution. det = 2·1 - 1·1 = 1. A^{-1} = (1/1)·[[1,-1],[-1,2]] = [[1,-1],[-1,2]]. Check (1,1): 2·1 + 1·(-1) = 1. NumPy: np.linalg.inv([[2,1],[1,1]]) -> array([[1.,-1.],[-1.,2.]]).
For matrices larger than 2×2, the systematic way to invert by hand is to row-reduce [A | I] until the left side becomes I.
Gauss-Jordan elimination finds an inverse by augmenting A with the identity to form [A | I], then applying row operations to drive the left block to reduced row echelon form. When the left side becomes I, the right side has become A^{-1}, because the same sequence of operations that turns A into I turns I into A^{-1}. If the left side cannot reach I (a zero row appears), then A is singular and no inverse exists — the method detects non-invertibility automatically. This generalizes the 2×2 formula to any size and reveals why invertibility is tied to full rank. It is the conceptual basis for how software computes inverses, and doing it by hand cements the link between elimination, rank, and inversion.
Worked Example 1
Problem. Invert A = [[2,0],[0,4]] by Gauss-Jordan.
Answer. [[0.5, 0], [0, 0.25]]
Worked Example 2
Problem. Invert A = [[1,1],[0,1]] by Gauss-Jordan.
Answer. [[1, -1], [0, 1]]
Worked Example 3
Problem. Begin inverting A = [[1,2,0],[0,1,0],[0,0,2]] by Gauss-Jordan and finish.
Answer. [[1, -2, 0], [0, 1, 0], [0, 0, 0.5]]
Problem. Invert A = [[1,3],[2,7]] using Gauss-Jordan (or the 2×2 shortcut to verify).
Solution. det = 1·7 - 3·2 = 1. Augment [[1,3,1,0],[2,7,0,1]]. R2 -> R2 - 2·R1 = [0,1,-2,1]. R1 -> R1 - 3·R2 = [1,0,7,-3]. So A^{-1} = [[7,-3],[-2,1]]. NumPy: np.linalg.inv([[1,3],[2,7]]) -> array([[7.,-3.],[-2.,1.]]).
NumPy gives fast numeric determinants and inverses; SymPy gives exact symbolic ones, and the condition number warns when inversion is unreliable.
np.linalg.det(A) returns the determinant numerically and np.linalg.inv(A) the inverse, both via LU decomposition; inv raises LinAlgError on a singular matrix. SymPy's Matrix(A).det() and .inv() give exact rational results, ideal for verifying hand work without rounding. A crucial numerical caveat is the condition number, np.linalg.cond(A): when it is large, A is nearly singular and its inverse amplifies tiny errors enormously, so results are untrustworthy. Best practice in data science is to avoid explicit inverses where possible — use np.linalg.solve or least-squares routines — and to inspect the condition number before trusting an inversion. Knowing both libraries lets you compute quickly and check exactly.
Worked Example 1
Problem. Compute the determinant and inverse of A = [[3,1],[1,2]] in NumPy.
Answer. det = 5; inverse = [[0.4, -0.2], [-0.2, 0.6]].
Worked Example 2
Problem. Verify the same inverse exactly with SymPy.
Answer. Exact inverse [[2/5,-1/5],[-1/5,3/5]], det = 5.
Worked Example 3
Problem. Inspect a near-singular matrix A = [[1,1],[1,1.0001]] with the condition number.
Answer. det ≈ 0.0001 and cond ≈ 40000 — near-singular, inverse unreliable.
Problem. Use NumPy to find det and inverse of A = [[2,5],[1,3]], then verify A @ A_inv ≈ I.
Solution. det = 2·3 - 5·1 = 1. A^{-1} = [[3,-5],[-1,2]]. NumPy: np.linalg.det([[2,5],[1,3]]) -> 1.0; Ai = np.linalg.inv([[2,5],[1,3]]) -> array([[3.,-5.],[-1.,2.]]); np.array([[2,5],[1,3]]) @ Ai -> array([[1.,0.],[0.,1.]]).
For a chosen 3x3 matrix, compute its determinant by cofactor expansion by hand, then find its inverse by Gauss-Jordan elimination. Verify both with np.linalg.det and np.linalg.inv, and confirm A @ A_inv is the identity (allowing for floating-point error). Finally, build a near-singular matrix, report its determinant and condition number, and comment on why inverting it numerically is risky.
Deliverable · A notebook with the hand determinant and inverse, the NumPy verification, and a short note on the near-singular case.
1. A square matrix is invertible if and only if its determinant is:
Answer B. Any nonzero determinant means full rank and an existing inverse; zero means singular.
2. The determinant of the 2x2 matrix [[a, b], [c, d]] is:
Answer C. The 2x2 determinant is the main-diagonal product minus the anti-diagonal product: ad - bc.
3. A determinant of 0 geometrically means the transformation:
Answer C. Zero volume scaling means the transformation flattens space, losing a dimension.
4. For invertible matrices, (AB)^{-1} equals:
Answer B. Inverting a product reverses the order, like undoing operations: (AB)^{-1} = B^{-1} A^{-1}.
5. To find an inverse by Gauss-Jordan, you reduce the augmented matrix:
Answer B. Row-reducing [A | I] until the left block becomes I leaves A^{-1} on the right.
I can compute a determinant and explain what a zero determinant implies.
I can find a matrix inverse by hand and verify A·A^{-1} = I.
I can use np.linalg.det and np.linalg.inv and recognize when a matrix is singular.
Eigenvectors are the special directions a matrix only stretches without rotating; the eigenvalue is the stretch factor.
An eigenvector of a square matrix A is a nonzero vector v whose direction is unchanged when A acts on it: Av = λv, meaning the output is just the input scaled by the eigenvalue λ. Geometrically, while A rotates and shears most vectors, eigenvectors lie along the axes that A merely stretches (or flips, if λ < 0) by the factor λ. These directions reveal a transformation's intrinsic structure: the long axis of a covariance matrix (the direction of most variance) is its top eigenvector, which is exactly how PCA finds principal components. Eigenvalues also govern the long-run behavior of repeated transformations (Markov chains, power iteration). Understanding Av = λv as 'A only scales v' is the conceptual core of the entire unit.
Worked Example 1
Problem. Verify that v = [1, 0] is an eigenvector of A = [[3,0],[0,5]] and find its eigenvalue.
Answer. Yes; eigenvalue λ = 3.
Worked Example 2
Problem. Check whether v = [1, 1] is an eigenvector of A = [[2,1],[1,2]].
Answer. Yes; eigenvalue λ = 3.
Worked Example 3
Problem. Show v = [1, 1] is NOT an eigenvector of A = [[1,2],[0,3]].
Answer. v = [1,2] is not an eigenvector (output not parallel to input).
Problem. Is v = [2, -2] an eigenvector of A = [[0,1],[1,0]]? If so, give λ.
Solution. Av: row1 = 0·2 + 1·(-2) = -2; row2 = 1·2 + 0·(-2) = 2, so Av = [-2, 2] = -1·[2,-2] = -v. Yes, it is an eigenvector with λ = -1. NumPy: np.array([[0,1],[1,0]]) @ np.array([2,-2]) -> array([-2, 2]).
Eigenvalues are the roots of det(A - λI) = 0, the characteristic polynomial that encodes a matrix's stretch factors.
Av = λv rearranges to (A - λI)v = 0. For a nonzero v to exist, A - λI must be singular, so its determinant must be zero: det(A - λI) = 0. Expanding this determinant gives the characteristic polynomial in λ, whose roots are precisely the eigenvalues. For a 2×2 matrix the polynomial is λ² - (trace)λ + det = 0, a quick shortcut. Each eigenvalue's multiplicity is how many times it repeats as a root. Two useful checks fall out: the eigenvalues sum to the trace and multiply to the determinant. This procedure converts the geometric idea of 'directions A only stretches' into a concrete algebraic equation you can solve, and it underlies how stability and variance directions are computed.
Worked Example 1
Problem. Find the eigenvalues of A = [[2,0],[0,3]].
Answer. λ = 2 and λ = 3
Worked Example 2
Problem. Find the eigenvalues of A = [[2,1],[1,2]] using trace and determinant.
Answer. λ = 1 and λ = 3
Worked Example 3
Problem. Find the eigenvalues of A = [[0,-1],[1,0]] (a 90-degree rotation).
Answer. λ = i and λ = -i (complex)
Problem. Find the eigenvalues of A = [[4,1],[2,3]].
Solution. trace = 4+3 = 7; det = 4·3 - 1·2 = 10. Characteristic: λ² - 7λ + 10 = 0 -> (λ-2)(λ-5) = 0, so λ = 2 and λ = 5. NumPy: np.linalg.eigvals([[4,1],[2,3]]) -> array([2., 5.]).
Once you have an eigenvalue, its eigenvectors are the nonzero solutions of (A - λI)v = 0 — the null space of A - λI.
For each eigenvalue λ, substitute it into A - λI and solve the homogeneous system (A - λI)v = 0. Because A - λI is singular by construction, this system has nonzero solutions; they form the eigenspace for λ (all solutions, plus the zero vector). In practice you reduce A - λI, find the free variable, and express the eigenvector direction. Eigenvectors are only determined up to scale, so any nonzero multiple is equally valid — software like NumPy normalizes them to length 1. Repeated eigenvalues may yield a multi-dimensional eigenspace with several independent eigenvectors. Finding these directions is what lets us diagonalize and is exactly the step PCA performs on a covariance matrix to extract principal directions.
Worked Example 1
Problem. Find an eigenvector of A = [[2,1],[1,2]] for λ = 3.
Answer. v = [1, 1] (or its normalization [0.707, 0.707])
Worked Example 2
Problem. Find an eigenvector of A = [[2,1],[1,2]] for λ = 1.
Answer. v = [1, -1]
Worked Example 3
Problem. Find an eigenvector of A = [[4,1],[2,3]] for λ = 5.
Answer. v = [1, 1]
Problem. Find an eigenvector of A = [[4,1],[2,3]] for λ = 2.
Solution. A - 2I = [[2,1],[2,1]]. The equation is 2v1 + v2 = 0, so v2 = -2v1. Take v1 = 1: v = [1, -2]. Check: Av = [4·1+1·(-2), 2·1+3·(-2)] = [2, -4] = 2·[1,-2]. Correct.
If a matrix has enough independent eigenvectors, it can be written as PDP^{-1}, turning hard matrix powers into easy diagonal powers.
A square matrix A is diagonalizable if it has n linearly independent eigenvectors. Stack them as the columns of P and put the matching eigenvalues on the diagonal of D; then A = PDP^{-1}. This factorization decouples the transformation: in the eigenvector coordinate system, A acts as simple independent scalings. Its biggest payoff is matrix powers: A^n = P D^n P^{-1}, and raising the diagonal D to the n-th power just raises each eigenvalue to the n-th power — vastly cheaper than multiplying A by itself n times. This is how long-run dynamics (Markov chains, recurrence relations like Fibonacci) are solved in closed form, and it is the discrete cousin of the SVD used throughout data science.
Worked Example 1
Problem. Diagonalize A = [[2,0],[0,3]].
Answer. P = I, D = diag(2, 3).
Worked Example 2
Problem. Diagonalize A = [[2,1],[1,2]] using its eigen-data.
Answer. P = [[1,1],[1,-1]], D = diag(3, 1).
Worked Example 3
Problem. Use diagonalization to compute A^10 for A = [[2,0],[0,3]].
Answer. A^10 = [[1024, 0], [0, 59049]]
Problem. Given A = PDP^{-1} with D = diag(1, 0.5), what is A^3 in terms of D?
Solution. A^3 = P D^3 P^{-1}, and D^3 = diag(1^3, 0.5^3) = diag(1, 0.125). So the eigenvalue 0.5 shrinks to 0.125 while 1 stays fixed. NumPy: np.diag([1,0.5])**3 conceptually gives diag(1, 0.125).
Symmetric matrices are the best-behaved: always real eigenvalues and a full set of orthogonal eigenvectors — the spectral theorem.
The spectral theorem says any real symmetric matrix (A = A^T) has only real eigenvalues and a complete set of mutually orthogonal eigenvectors. Normalizing those eigenvectors to length 1 gives an orthonormal matrix Q, so A = Q D Q^T — the eigenvectors form a perfect rotated coordinate system in which A is pure scaling. This is gold for data science because covariance matrices and Gram matrices A^T A are always symmetric and positive semidefinite (nonnegative eigenvalues), guaranteeing real, orthogonal principal directions. PCA is exactly the spectral decomposition of a covariance matrix: the eigenvectors are the principal components and the eigenvalues are the variances along them. The orthogonality means the components are uncorrelated, which is why PCA produces clean, independent axes.
Worked Example 1
Problem. Confirm the eigenvectors of the symmetric A = [[2,1],[1,2]] are orthogonal.
Answer. Dot product 0 — eigenvectors are orthogonal.
Worked Example 2
Problem. Show the symmetric A = [[5,4],[4,5]] has real eigenvalues.
Answer. λ = 1 and λ = 9, both real.
Worked Example 3
Problem. For the covariance-like matrix A = [[3,0],[0,1]], identify principal directions and variances.
Answer. Principal directions [1,0] (variance 3) and [0,1] (variance 1).
Problem. Find the eigenvalues of the symmetric A = [[4,2],[2,4]] and confirm they are real.
Solution. trace = 8, det = 16 - 4 = 12. Characteristic: λ² - 8λ + 12 = 0 -> (λ-2)(λ-6) = 0, so λ = 2 and λ = 6, both real. NumPy: np.linalg.eigvalsh([[4,2],[2,4]]) -> array([2., 6.]).
NumPy computes eigenvalues and eigenvectors directly; knowing its output conventions avoids subtle bugs.
np.linalg.eig(A) returns a tuple (w, V): w is a 1-D array of eigenvalues and V is a matrix whose COLUMNS are the corresponding unit-length eigenvectors — V[:, i] pairs with w[i]. The eigenvectors are normalized to length 1 and may carry an arbitrary sign, so they can differ from hand calculations by a scalar; that is fine since direction is what matters. For symmetric matrices, prefer np.linalg.eigh, which is faster, returns real sorted eigenvalues, and yields a genuinely orthonormal eigenvector matrix. To verify, reconstruct A as V @ diag(w) @ inv(V) (or V @ diag(w) @ V.T for symmetric). Understanding that eigenvectors are columns, normalized, and sign-ambiguous prevents the most common eigen-decomposition mistakes in practice.
Worked Example 1
Problem. Use np.linalg.eig on A = [[2,0],[0,3]] and read off the pairs.
Answer. Eigenvalues [2, 3] with eigenvectors [1,0] and [0,1].
Worked Example 2
Problem. Reconstruct A = [[2,1],[1,2]] from its eigendecomposition.
Answer. Reconstruction recovers A = [[2,1],[1,2]].
Worked Example 3
Problem. Compare np.linalg.eig and np.linalg.eigh on symmetric A = [[2,1],[1,2]].
Answer. Both give eigenvalues {1, 3}; eigh sorts them and guarantees orthonormal vectors.
Problem. After w, V = np.linalg.eig(A), how do you get the eigenvector for the largest eigenvalue?
Solution. Find the index of the max eigenvalue with i = np.argmax(w), then take the column V[:, i]. For example if w = [1., 3.], argmax gives i = 1, so V[:, 1] is the eigenvector for λ = 3. Remember eigenvectors are columns and are unit length.
For a 2x2 matrix, find the eigenvalues from the characteristic polynomial and the matching eigenvectors by hand. Confirm with np.linalg.eig, noting that NumPy returns unit-length eigenvectors. Then diagonalize the matrix as PDP^{-1} and use it to compute A raised to the 10th power efficiently, checking against np.linalg.matrix_power.
Deliverable · A notebook with the hand eigen-computation, the NumPy verification, the diagonalization, and the fast matrix power.
1. An eigenvector v of A satisfies:
Answer B. Eigenvectors keep their direction under A, only scaling by the eigenvalue λ.
2. Eigenvalues are found by solving:
Answer B. Nontrivial eigenvectors exist only when A - λI is singular, i.e. det(A - λI) = 0.
3. A symmetric real matrix is guaranteed to have:
Answer C. The spectral theorem guarantees real eigenvalues and an orthogonal set of eigenvectors.
4. The sum of a matrix's eigenvalues equals its:
Answer C. The trace (sum of diagonal entries) equals the sum of the eigenvalues.
5. Diagonalizing A as PDP^{-1} makes computing A^n easy because:
Answer B. A^n = P D^n P^{-1}, and raising a diagonal D to a power just powers each diagonal entry.
I can find eigenvalues and eigenvectors of a small matrix by hand.
I can diagonalize a matrix and use it to compute matrix powers quickly.
I can interpret eigenvectors as directions a transformation only stretches.
Orthogonal vectors carry independent information; making them unit length too gives an orthonormal set, the ideal coordinate system.
A set of vectors is orthogonal if every pair has dot product zero — they meet at right angles and share no overlapping direction. It is orthonormal if, in addition, each vector has length 1. Orthonormal sets are the perfect basis: coordinates are found by simple dot products (no solving systems), the basis matrix Q satisfies Q^T Q = I, and Q^T is its own inverse. This is why orthogonality is prized in data science: uncorrelated features behave independently, and orthonormal bases (as in PCA and the SVD) make projections and reconstructions trivial and numerically stable. Recognizing and constructing orthonormal sets is the foundation for projections, Gram-Schmidt, QR, and ultimately the decompositions that compress data.
Worked Example 1
Problem. Are u = [1, 0] and v = [0, 1] orthonormal?
Answer. Yes, orthonormal.
Worked Example 2
Problem. Are u = [1, 1] and v = [1, -1] orthogonal? Are they orthonormal?
Answer. Orthogonal, but not orthonormal (norms are sqrt(2)).
Worked Example 3
Problem. Normalize u = [1, 1] and v = [1, -1] into an orthonormal set.
Answer. q1 = [0.707, 0.707], q2 = [0.707, -0.707].
Problem. Are u = [2, 1] and v = [-1, 2] orthogonal? Normalize them if so.
Solution. Dot = 2·(-1) + 1·2 = 0, so orthogonal. Each has norm sqrt(5). Orthonormal versions: u/sqrt(5) = [0.894, 0.447], v/sqrt(5) = [-0.447, 0.894]. NumPy: np.dot([2,1],[-1,2]) -> 0.
Projection finds the closest point in a subspace to a vector; the leftover error (residual) is orthogonal to that subspace.
Projecting a vector b onto a line spanned by a gives the point proj = ((a·b)/(a·a))·a — the shadow of b along a. The residual b - proj is orthogonal to a, capturing the part of b that the line cannot represent. Projecting onto a subspace spanned by the columns of a matrix A generalizes this: proj = A(A^T A)^{-1} A^T b, where P = A(A^T A)^{-1} A^T is the projection matrix. Projection is the geometric heart of least squares — fitting a model is projecting the data onto the space your model can reach, and the residual measures the unavoidable error. This idea drives regression, dimensionality reduction, and signal denoising throughout data science.
Worked Example 1
Problem. Project b = [3, 4] onto the line spanned by a = [1, 0].
Answer. Projection [3, 0]; residual [0, 4].
Worked Example 2
Problem. Project b = [2, 2] onto the line spanned by a = [1, 1].
Answer. Projection [2, 2]; residual [0, 0].
Worked Example 3
Problem. Project b = [1, 2] onto the line spanned by a = [3, 4].
Answer. Projection [1.32, 1.76]; residual [-0.32, 0.24].
Problem. Project b = [4, 0] onto the line spanned by a = [1, 1].
Solution. a·b = 1·4 + 1·0 = 4; a·a = 2. scalar = 4/2 = 2. proj = 2·[1,1] = [2, 2]. Residual = [4,0] - [2,2] = [2,-2], orthogonal to [1,1] (dot = 2-2 = 0). NumPy: (np.dot([1,1],[4,0])/2)*np.array([1,1]) -> array([2,2]).
Gram-Schmidt converts any basis into an orthonormal one; bundling the result gives the QR factorization used for stable least squares.
Gram-Schmidt takes a set of independent vectors and produces an orthonormal set spanning the same space. It works greedily: keep the first vector (normalized), then for each next vector subtract its projections onto all previously chosen directions (removing overlap) and normalize the remainder. Collecting the orthonormal vectors as columns of Q and recording the coefficients in an upper-triangular R yields the QR factorization A = QR, where Q has orthonormal columns (Q^T Q = I) and R is upper triangular. QR is the numerically stable workhorse for least squares: instead of forming A^T A (which squares the condition number), you solve Rx = Q^T b. Understanding Gram-Schmidt explains where orthonormal bases come from and why QR is preferred in practice.
Worked Example 1
Problem. Apply Gram-Schmidt to a1 = [1, 0] and a2 = [1, 1].
Answer. Orthonormal basis q1 = [1,0], q2 = [0,1].
Worked Example 2
Problem. Orthogonalize a1 = [3, 0] and a2 = [2, 2] (normalize a1 only conceptually).
Answer. q1 = [1, 0], q2 = [0, 1].
Worked Example 3
Problem. Factor A = [[1,1],[1,0]] as QR with NumPy and confirm Q^T Q = I.
Answer. A = QR with orthonormal Q (Q^T Q = I) and upper-triangular R.
Problem. Apply Gram-Schmidt to a1 = [0, 3] and a2 = [4, 1].
Solution. q1 = [0,3]/3 = [0,1]. proj of a2 onto q1 = (a2·q1)q1 = (1)[0,1] = [0,1]; w = [4,1] - [0,1] = [4,0]. Normalize: q2 = [1,0]. Orthonormal basis: q1 = [0,1], q2 = [1,0].
When there are more equations than unknowns, b usually lies outside the column space, so no exact solution exists — we seek the best approximation instead.
An overdetermined system has more equations than unknowns (a tall, thin A). Its solvability depends on whether b lies in the column space of A — the set of all vectors A can produce. With real, noisy data, b almost never lands exactly in that lower-dimensional space, so Ax = b has no exact solution. Rather than give up, we ask for the x that gets Ax as close to b as possible, minimizing the residual ||Ax - b||. Geometrically, the best Ax is the projection of b onto the column space, and the residual is orthogonal to that space. This shift from 'solve exactly' to 'minimize error' is the conceptual leap into least squares — the standard tool for fitting models to imperfect data.
Worked Example 1
Problem. Explain why fitting a line through 3 non-collinear points (x,y) = (0,1),(1,3),(2,2) has no exact solution.
Answer. No line passes through all three points; the system is inconsistent (overdetermined).
Worked Example 2
Problem. For A = [[1],[1],[1]] and b = [2, 4, 6], is Ax = b solvable exactly?
Answer. Not solvable exactly; the least-squares answer x = 4 (the mean).
Worked Example 3
Problem. Verify with NumPy that A = [[1,0],[0,1],[1,1]], b = [1,1,3] has no exact solution but a least-squares one.
Answer. Inconsistent exactly; least-squares solution ≈ [1.333, 1.333].
Problem. Does A = [[1],[1]] with b = [3, 5] have an exact solution? What is the best single value x?
Solution. A's column space is multiples of [1,1]; b = [3,5] is not such a multiple, so no exact x. The least-squares x is the mean of b: (3+5)/2 = 4, giving Ax = [4,4]. NumPy: np.linalg.lstsq([[1],[1]], [3,5], rcond=None)[0] -> array([4.]).
Setting the residual orthogonal to the column space yields the normal equations, the algebraic recipe for the best-fit solution.
Least squares seeks x minimizing ||Ax - b||. The minimizer makes the residual r = b - Ax orthogonal to the column space of A, i.e. A^T r = 0, which rearranges to the normal equations A^T A x = A^T b. When A has independent columns, A^T A is square, symmetric, and invertible, so x = (A^T A)^{-1} A^T b. These equations are how linear regression coefficients are computed in closed form: A is the design matrix (features), x the coefficients, b the targets. The matrix (A^T A)^{-1} A^T is the pseudoinverse. While conceptually central, forming A^T A squares the condition number, so numerically QR or SVD is preferred — but the normal equations remain the clearest explanation of why the fit is best.
Worked Example 1
Problem. Fit a constant to data b = [1, 2, 3] using A = [[1],[1],[1]] via the normal equations.
Answer. x = 2 (the mean).
Worked Example 2
Problem. Fit y = mx through the origin to points (1,2),(2,2),(3,4) via normal equations.
Answer. m ≈ 1.286
Worked Example 3
Problem. Fit y = mx + c to (0,1),(1,3),(2,2) via the normal equations.
Answer. Best-fit line y = 0.5x + 1.5.
Problem. Fit y = mx + c to (0,0),(1,1),(2,3) using the normal equations.
Solution. A = [[0,1],[1,1],[2,1]], b = [0,1,3]. A^T A = [[5,3],[3,3]]; A^T b = [0+1+6, 0+1+3] = [7, 4]. Solve [[5,3],[3,3]][m,c]=[7,4]: subtract -> 2m = 3 so m = 1.5; then 3·1.5+3c = 4 -> c = -1/6 ≈ -0.167. Line y = 1.5x - 0.167. NumPy: np.linalg.solve(A.T@A, A.T@b) -> array([1.5, -0.167]).
np.linalg.lstsq solves least squares directly and stably, returning the best-fit coefficients plus diagnostics like residuals and rank.
np.linalg.lstsq(A, b, rcond=None) computes the least-squares solution to Ax ≈ b using a stable SVD-based routine, avoiding the explicit A^T A that can be ill-conditioned. It returns four items: the solution x, the sum of squared residuals, the rank of A, and the singular values. To fit a line you build the design matrix with a column for x and a column of ones for the intercept, stack the targets in b, and call lstsq. This is the practical, production-grade way to fit linear models and is what underlies tools like NumPy's polyfit and the linear core of scikit-learn. Knowing how to assemble the design matrix and read lstsq's outputs is a direct, employable data-science skill.
Worked Example 1
Problem. Use lstsq to fit y = mx + c to (0,1),(1,3),(2,2) and confirm it matches the normal-equation answer.
Answer. y = 0.5x + 1.5 (matches the normal equations).
Worked Example 2
Problem. Fit a line to (1,1),(2,2),(3,2),(4,3) and report the slope and intercept.
Answer. y = 0.6x + 0.5
Worked Example 3
Problem. Verify the residual from a lstsq fit is orthogonal to the columns of A for (0,1),(1,3),(2,2).
Answer. A^T r = [0, 0]: the residual is orthogonal to A's columns, confirming the fit.
Problem. Fit y = mx + c to (1,2),(2,4),(3,5) with np.linalg.lstsq.
Solution. A = np.column_stack([[1,2,3], np.ones(3)]); b = [2,4,5]. np.linalg.lstsq(A, b, rcond=None)[0] -> array([1.5, 0.667]). So y = 1.5x + 0.667. (By hand: A^T A = [[14,6],[6,3]], A^T b = [25,11]; solving gives m = 1.5, c = 2/3.)
Generate or take a small set of (x, y) data points that roughly follow a line plus noise. Build the design matrix A and set up the normal equations A^T A x = A^T b by hand for the slope and intercept. Solve them with NumPy, then compare to np.linalg.lstsq, and confirm the residual vector is orthogonal to the columns of A. Plot the data and the fitted line.
Deliverable · A notebook with the normal-equation setup, the solved coefficients, the lstsq comparison, the orthogonality check, and the plot.
1. The least squares solution minimizes:
Answer B. Least squares finds the x making Ax as close as possible to b, minimizing the residual length.
2. The normal equations for least squares are:
Answer B. Setting the residual orthogonal to the column space yields A^T A x = A^T b.
3. The residual vector at the least squares solution is:
Answer C. The best approximation leaves an error perpendicular to the subspace it was projected onto.
4. Gram-Schmidt converts a basis into one that is:
Answer B. Gram-Schmidt produces mutually orthogonal, unit-length vectors — an orthonormal basis.
5. In the QR factorization A = QR, the matrix Q has:
Answer A. Q has orthonormal columns; R is upper triangular, making least squares numerically stable.
I can project a vector onto a subspace and compute the residual.
I can set up and solve the least squares normal equations for a best-fit line.
I can fit a model to noisy data with np.linalg.lstsq and interpret the result.
The SVD factors any matrix into a rotation, a scaling, and another rotation — the most general and useful decomposition in data science.
The singular value decomposition writes any m×n matrix A as A = U Σ V^T. Here U (m×m) and V (n×n) are orthogonal (their columns are orthonormal), and Σ is m×n diagonal holding the nonnegative singular values σ1 ≥ σ2 ≥ ... in decreasing order. Geometrically, V^T rotates the input, Σ scales along the new axes, and U rotates into the output space — every linear map is a rotate-stretch-rotate. Unlike eigendecomposition, the SVD exists for every matrix (rectangular, singular, anything). The columns of V are the right singular vectors (input directions), the columns of U the left singular vectors, and the singular values measure each direction's importance. The SVD underpins PCA, low-rank compression, recommender systems, and the stable pseudoinverse — it is the central tool of this unit.
Worked Example 1
Problem. Interpret the SVD of the diagonal matrix A = [[3,0],[0,2]].
Answer. U = V = I, singular values 3 and 2.
Worked Example 2
Problem. Compute the SVD of A = [[0,2],[0,0]] and read off the singular values.
Answer. Singular values 2 and 0 (rank 1).
Worked Example 3
Problem. Reconstruct A = [[1,0],[0,1],[1,1]] from its SVD in NumPy.
Answer. U @ diag(S) @ Vt reproduces A; singular values ≈ [1.732, 1.0].
Problem. What are the singular values of A = [[5,0],[0,0],[0,1]]?
Solution. The nonzero rows scale by 5 and by 1 along orthogonal directions, so the singular values are 5 and 1 (sorted: 5, 1). NumPy: np.linalg.svd([[5,0],[0,0],[0,1]])[1] -> array([5., 1.]). The matrix has rank 2.
Singular values rank directions by importance; keeping only the largest gives the best low-rank approximation, the basis of compression.
The number of nonzero singular values equals the rank of A — they count its genuinely independent directions. Because Σ orders them from largest to smallest, the top singular values carry the most of the matrix's 'energy' (its Frobenius norm squared equals the sum of squared singular values). The Eckart-Young theorem says the best rank-k approximation of A, in least-squares sense, is obtained by keeping only the top k singular values and their singular vectors: A_k = Σ_{i=1}^k σi ui vi^T. Dropping small singular values discards the least important structure, which is exactly how images, recommender matrices, and embeddings are compressed and denoised. This connects directly to PCA: keeping top components is keeping top singular directions.
Worked Example 1
Problem. What is the rank of a matrix whose singular values are 5, 3, 0, 0?
Answer. rank = 2
Worked Example 2
Problem. A matrix has singular values [10, 1, 0.1]. What fraction of energy does the top component hold?
Answer. ≈ 99.0% of the energy is in the top singular direction.
Worked Example 3
Problem. Build a rank-1 approximation of A = [[2,2],[2,2]] from its SVD.
Answer. A_1 = [[2,2],[2,2]] — the rank-1 approximation equals A exactly.
Problem. Singular values are [8, 6, 0]. Give the rank and the energy fraction captured by the top two.
Solution. Two nonzero values, so rank = 2. Total energy = 8² + 6² + 0² = 64 + 36 = 100. Top two capture (64+36)/100 = 1.0 = 100% (the third is zero). NumPy: S=np.array([8,6,0]); np.sum(S[:2]**2)/np.sum(S**2) -> 1.0.
The covariance matrix summarizes how features vary together; its eigenvectors are the principal components and its eigenvalues the variances.
For centered data (each feature minus its mean), the covariance matrix C = (1/(n-1)) X^T X is symmetric and captures how every pair of features co-varies: diagonal entries are individual variances, off-diagonals are covariances. By the spectral theorem, C has real eigenvalues and orthogonal eigenvectors. The eigenvectors are the principal components — orthogonal directions in feature space — and each eigenvalue is the variance of the data along its component. The largest eigenvalue points along the direction of greatest spread. PCA is exactly this eigendecomposition: it finds the axes that explain the most variance and are mutually uncorrelated. Because C is built from X^T X, its eigenvectors coincide with the right singular vectors of X, which is why SVD computes PCA directly and more stably.
Worked Example 1
Problem. Compute the covariance matrix of centered data X = [[1,1],[-1,-1]] (already mean-zero per column).
Answer. C = [[2, 2], [2, 2]]
Worked Example 2
Problem. Find the principal components of C = [[2,0],[0,1]].
Answer. PC1 = [1,0] (variance 2), PC2 = [0,1] (variance 1).
Worked Example 3
Problem. Find the top principal component of C = [[2,2],[2,2]].
Answer. PC1 = [0.707, 0.707] with variance 4 (data lies on one line).
Problem. Given C = [[5,0],[0,2]], list the principal components and the variance each explains.
Solution. C is diagonal: eigenvalues 5 and 2 with eigenvectors [1,0] and [0,1]. PC1 = [1,0] explains variance 5; PC2 = [0,1] explains variance 2. Total variance 7, so PC1 explains 5/7 ≈ 71.4%. NumPy: np.linalg.eigh([[5,0],[0,2]]) -> eigenvalues array([2., 5.]).
PCA reduces dimension by projecting centered data onto its top principal components, keeping the directions of greatest variance.
Given centered data X, PCA projects each point onto the top k principal components (the top right singular vectors of X, equivalently the top eigenvectors of the covariance matrix). The projection coordinates are the new low-dimensional representation, called scores: if V_k holds the top k components as columns, the scores are X V_k. This keeps the directions where data varies most and discards the rest, compressing many correlated features into a few uncorrelated ones while preserving structure. Because the components are orthonormal, the projection is just dot products — fast and stable. PCA is the workhorse of dimensionality reduction: it powers visualization (project to 2D), noise reduction, and preprocessing before clustering or regression, all by projecting onto the most informative axes.
Worked Example 1
Problem. Project the centered point [3, 1] onto the principal component v = [0.707, 0.707].
Answer. Score ≈ 2.828 (its coordinate along PC1).
Worked Example 2
Problem. Project two centered points X = [[2,0],[0,2]] onto PC1 = [0.707, 0.707].
Answer. Scores [1.414, 1.414].
Worked Example 3
Problem. Reduce centered X = [[1,1],[2,2],[3,3]] to 1D using its top component.
Answer. 1D scores [1.414, 2.828, 4.243], losing no information.
Problem. Project centered points [[4,0],[0,4]] onto PC1 = [0.707, 0.707].
Solution. Point [4,0]: 4·0.707 + 0 = 2.828. Point [0,4]: 0 + 4·0.707 = 2.828. Scores along PC1: [2.828, 2.828]. NumPy: np.array([[4,0],[0,4]]) @ np.array([0.707,0.707]) -> array([2.828, 2.828]).
The explained variance ratio tells you what share of total variance each component keeps, guiding how many to retain.
Each principal component's eigenvalue (or singular value squared) is the variance it captures. The explained variance ratio for component i is its eigenvalue divided by the sum of all eigenvalues — the fraction of total variance it accounts for. Summing these in order gives the cumulative explained variance, and you typically keep just enough components to reach a target (say 90% or 95%). A scree plot of the ratios shows where they level off (the 'elbow'), suggesting a natural cutoff. This is the practical decision in dimensionality reduction: too few components lose signal, too many keep noise. Choosing by explained variance balances compression against fidelity and is a standard, defensible criterion in real PCA pipelines.
Worked Example 1
Problem. Eigenvalues are [6, 3, 1]. Compute each component's explained variance ratio.
Answer. [0.6, 0.3, 0.1] (60%, 30%, 10%).
Worked Example 2
Problem. Using eigenvalues [6, 3, 1], how many components are needed to retain at least 85% of variance?
Answer. 2 components (cumulative 90%).
Worked Example 3
Problem. Singular values are [4, 2, 1]. Compute the explained variance ratios (use squared singular values).
Answer. ≈ [0.762, 0.190, 0.048]; two components reach ~95%.
Problem. Eigenvalues [10, 5, 4, 1]. How many components retain at least 75% of variance?
Solution. Total = 20. Ratios: 0.5, 0.25, 0.2, 0.05. Cumulative: 0.5, 0.75, 0.95, 1.0. The 75% target is met at 2 components (cumulative exactly 0.75). NumPy: ev=np.array([10,5,4,1]); np.cumsum(ev/ev.sum()) -> array([0.5,0.75,0.95,1.0]).
NumPy's svd gives the raw decomposition; scikit-learn's PCA wraps centering, projection, and explained variance into one tested tool.
np.linalg.svd(X, full_matrices=False) returns U, S (singular values), and Vt for a data matrix X; after centering X, the rows of Vt are the principal components and S²/(n-1) are the component variances. scikit-learn's PCA automates this: PCA(n_components=k).fit(X) centers the data internally, computes the SVD, and exposes .components_ (the principal directions), .explained_variance_ratio_ (the fractions), and .transform(X) (the projected scores). Use raw SVD when you want full control or to understand the mechanics; use sklearn PCA in production for its convenience, correctness, and integration with pipelines. Both give the same answer up to sign conventions on the singular vectors. Knowing both lets you verify by hand and ship reliably.
Worked Example 1
Problem. Get the principal components of centered X = [[1,0],[-1,0],[0,1],[0,-1]] via np.linalg.svd.
Answer. Components along [1,0] and [0,1] with equal variance.
Worked Example 2
Problem. Use scikit-learn PCA to reduce X = [[2,0],[0,2],[-2,0],[0,-2]] to 1 component.
Answer. 1 component captures 50% of variance; transform gives the 1-D scores.
Worked Example 3
Problem. Confirm np.linalg.svd and sklearn PCA agree on explained variance for centered X.
Answer. Both report identical explained variance ratios.
Problem. After U,S,Vt = np.linalg.svd(X_centered, full_matrices=False), how do you get the explained variance ratio?
Solution. Square the singular values and normalize: ratios = (S**2) / np.sum(S**2). For example if S = [4, 3], then S**2 = [16, 9], total 25, so ratios = [0.64, 0.36]. This matches sklearn's PCA(...).explained_variance_ratio_.
Take a small numeric dataset (e.g. a handful of features), center it by subtracting the column means, and compute its SVD with np.linalg.svd. Project the data onto the top two principal components and report the explained variance ratio for each. Reconstruct a low-rank approximation using only the top components and compare it to the original, then verify your result against scikit-learn's PCA.
Deliverable · A notebook with the centering step, the SVD, the 2D projection plot, the explained-variance report, and the sklearn cross-check.
1. In the SVD A = U Σ V^T, the matrix Σ contains:
Answer B. Σ is diagonal with the nonnegative singular values, ordered from largest to smallest.
2. Before running PCA you should first:
Answer B. PCA describes variance about the mean, so each feature must be centered first.
3. Principal components are directions of:
Answer B. PCA finds orthogonal axes capturing the most variance, ranked by singular value.
4. A low-rank approximation from the SVD keeps:
Answer C. Dropping small singular values and keeping the largest gives the best low-rank approximation.
5. Explained variance ratio tells you:
Answer B. It reports the fraction of total variance each principal component accounts for, guiding how many to keep.
I can compute and interpret the SVD of a data matrix.
I can use PCA to reduce dimensionality and explain the variance retained.
I can build a low-rank approximation and connect it to compression.
Where this leads
Part of the open-source Crunch Academy tech-education path. Pair it with the data-science and ML tracks when you're ready.
All Math for Data Science courses Crunch Academy home