Linear Algebra for ML

You're investigating a production incident during a high-traffic release. Pull three customer support tickets into an incident table. Each row tracks two binary signals: did the user report high latency, and did they request an emergency rollback? Ticket t2 flags both issues at once. Intuitively, the joint latency-plus-rollback pattern dominates your dataset, while the contrast between latency-only and rollback-only looks like a secondary detail.

What happens if you compress this table to save memory or speed up vector retrieval? If you discard the weaker contrast direction, a latency-only ticket and a rollback-only ticket collapse into identical representations. Your search or routing pipeline suddenly can't tell an API timeout from a broken deployment.

Singular value decomposition (SVD) finds perpendicular input directions and measures the output length produced by each unit direction. For this uncentered table, those strengths describe magnitude relative to zero, rather than variation around the average ticket. The decision is which of those directions to keep, and the cost is a number we can calculate from this table.

In the vectors, matrices, and tensors lesson, we worked with rows, columns, and dot products. Keep those operations close: three tickets and two features give us everything we need to compute, visualize, and diagnose a full SVD by hand.

A matrix can hide repeated directions

Write those flags as a 3×2 matrix. Each row is one ticket; each column records a latency or rollback signal:

Look at t2 and predict its score along two diagonals. Because it contains both signals, it should reinforce [latency + rollback]; its contribution to [latency - rollback] should cancel.

First compare the raw patterns [1, 1] and [1, -1]. Ticket t2 = [1, 1] scores 1 + 1 = 2 on the first and 1 - 1 = 0 on the second. Both pattern vectors have length $\sqrt{1^2 + 1^2} = \sqrt{2}$. Divide by that length to compare unit vectors, whose length is 1:

Multiply A by either direction. Each entry is one ticket's score for that pattern. The common scores are $[1/\sqrt{2}, 1/\sqrt{2}, 2/\sqrt{2}]$, so their squared lengths add to $1/2 + 1/2 + 2 = 3$. For contrast, the scores are $[1/\sqrt{2}, -1/\sqrt{2}, 0]$ and the squared lengths add to $1$. Taking square roots gives strengths $\sqrt{3}$ and $1$. NumPy checks each step:

import numpy as np

A = np.array([
    [1.0, 0.0],  # latency
    [0.0, 1.0],  # rollback
    [1.0, 1.0],  # latency and rollback
])
common = np.array([1.0, 1.0]) / np.sqrt(2)
contrast = np.array([1.0, -1.0]) / np.sqrt(2)

common_scores = A @ common
contrast_scores = A @ contrast

print("A_shape", A.shape)
print("common_scores", common_scores.round(4).tolist())
print("contrast_scores", contrast_scores.round(4).tolist())
print("strengths", round(np.linalg.norm(common_scores), 4), round(np.linalg.norm(contrast_scores), 4))

A_shape (3, 2)
common_scores [0.7071, 0.7071, 1.4142]
contrast_scores [0.7071, -0.7071, 0.0]
strengths 1.7321 1.0

The common direction has strength $\sqrt{3}$; the contrast direction has strength $1$. Nothing has been compressed yet. This is a measurement of which feature mix the rows support most strongly.

View larger

Coordinates come from a basis

The common and contrast directions are linear combinations of the original feature axes: scale each axis, then add. The span is every vector reachable by those combinations.

These directions are orthogonal: their dot product is $(1 - 1)/2 = 0$, so they meet at a right angle. Neither is a multiple of the other, which makes them linearly independent. Together they form a basis: a set of independent vectors that spans the entire two-dimensional feature space. Because they're also unit vectors, this basis is orthonormal.

To practice changing coordinates, allow mention counts instead of yes-or-no flags for a new row $x = [3, 1]$. Its common score is $(3 + 1)/\sqrt{2} = 2\sqrt{2}$; its contrast score is $(3 - 1)/\sqrt{2} = \sqrt{2}$.

Because this basis is orthonormal, its coordinates are dot products with the two unit directions:

Each coefficient is an orthogonal projection coordinate. Multiply a coefficient by its unit direction to recover the part of x carried by that pattern:

So $x = (x \cdot v_{\text{common}})v_{\text{common}} + (x \cdot v_{\text{contrast}})v_{\text{contrast}}$. The A @ v scores above measure every ticket along one direction; these coordinates rebuild one ticket from its directions. Later, projecting onto fewer SVD directions will keep selected coordinates and deliberately discard the rest.

What if you use b = [1, 1] for the common direction without making it a unit vector? The line is the same, but its coordinate scale changes. For any nonzero real direction b, the projection is

The denominator matters. To find the nearest point c * b, require its residual to be perpendicular to the line: $(x - cb)\cdot b = 0$. Rearranging gives $c = (x\cdot b)/(b\cdot b)$. For our row, that's 4 / 2 = 2, so the projection is still [2, 2]. With a unit direction, b · b = 1, which recovers the earlier formula.^{[1]Reference 1Numerical Linear Algebra.https://people.maths.ox.ac.uk/trefethen/books.html}

import numpy as np

x = np.array([3.0, 1.0])
b = np.array([1.0, 1.0])
projection = (x @ b) / (b @ b) * b
residual = x - projection
wrong_projection = (x @ b) * b

print("projection", projection.tolist())
print("residual", residual.tolist())
print("residual_dot_direction", float(residual @ b))
print("without_denominator", wrong_projection.tolist())
assert np.allclose(residual @ b, 0)
assert not np.allclose((x - wrong_projection) @ b, 0)

projection [2.0, 2.0]
residual [1.0, -1.0]
residual_dot_direction 0.0
without_denominator [4.0, 4.0]

The perpendicular residual is a useful check: the discarded [1, -1] contrast can't move you any closer along the common line. Using [4, 4] instead overshoots it.

We chose this basis by inspection. For real data, you need the matrix to reveal it. A.T @ A makes that choice visible without changing the feature space.

The Gram matrix reveals those directions by hand

For a data matrix A, compare feature columns with one another by multiplying on the feature side:

The first diagonal entry is the latency column's dot product with itself: $1^2 + 0^2 + 1^2 = 2$. The off-diagonal entry is the columns' overlap: $1\cdot0 + 0\cdot1 + 1\cdot1 = 1$. G has shape [2, 2] because its rows and columns both index features.

A nonzero vector v satisfying $Gv = \lambda v$ is an eigenvector. Multiplying by G keeps it on the same line through the origin. Here the eigenvalue $\lambda$ is nonnegative because G is a Gram matrix, so it stretches or shrinks the vector without reversing it. For a unit eigenvector, $\lambda$ is squared strength because

when v has length 1.

For this matrix, the common direction is scaled by 3 and the contrast direction by 1:

View larger

import numpy as np

A = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
G = A.T @ A
common = np.array([1.0, 1.0]) / np.sqrt(2)
contrast = np.array([1.0, -1.0]) / np.sqrt(2)

print("G", G.tolist())
print("G_common", (G @ common).round(4).tolist())
print("three_common", (3 * common).round(4).tolist())
print("G_contrast_equals_contrast", np.allclose(G @ contrast, contrast))

G [[2.0, 1.0], [1.0, 2.0]]
G_common [2.1213, 2.1213]
three_common [2.1213, 2.1213]
G_contrast_equals_contrast True

SVD names the directions and strengths

Any real matrix, including this 3×2 incident table, has a singular value decomposition (SVD):^{[1]Reference 1Numerical Linear Algebra.https://people.maths.ox.ac.uk/trefethen/books.html}

The incident matrix A has shape [3, 2]. Its reduced SVD keeps min(3, 2) = 2 directions. Reduced SVD still includes zero singular values if the matrix has dependent columns; it hasn't made a compression decision yet.^{[2]Reference 2numpy.linalg.svdhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html}

Here, the right singular vectors are the normalized common and contrast directions. Their singular values tell us how strongly A stretches those directions:

For each positive singular value, the corresponding left singular vector is $u_i = A v_i / \sigma_i$. For the common pattern, dividing its three scores by $\sqrt{3}$ gives $u_1 = [1, 1, 2]^T/\sqrt{6}$. For contrast, $u_2 = [1, -1, 0]^T/\sqrt{2}$. Columns of U have unit length; the scaled columns U * S contain the actual ticket scores. Put all three factors back together and they reconstruct every original entry.

For computation, call the SVD routine directly. It returns U, the one-dimensional singular-value array S, and Vt, which is $V^T$. Before running it, predict the shapes: reduced U should have one row per ticket, and Vt should have one row per retained feature direction.

import numpy as np

A = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
reconstructed = U @ np.diag(S) @ Vt

print("shapes", U.shape, S.shape, Vt.shape)
print("singular_values", S.round(4).tolist())
print("reconstruction_error", float(np.max(np.abs(A - reconstructed))))
print("reconstructs_A", np.allclose(A, reconstructed))

shapes (3, 2) (2,) (2, 2)
singular_values [1.7321, 1.0]
reconstruction_error 4.440892098500626e-16
reconstructs_A True

The sign of a singular vector may differ between correct implementations. If matching directions in U and V flip sign together, their product stays the same. Repeated singular values also permit rotations within the tied subspace. Compare singular values and reconstruction, not raw vectors. The tiny printed reconstruction error is floating-point roundoff, not discarded information.

Rank counts independent patterns

The rank of a matrix is the number of non-zero singular values. It counts how many independent directions the matrix really needs, not how many columns happen to be present.

Now make the failure concrete. Suppose every rollback feature exactly repeats the latency feature:

The second column adds no new information. B still has two columns, but only one independent direction, so its rank is 1. If you had to predict a new row from one column, the duplicate couldn't improve that prediction.

import numpy as np

B = np.array([
    [1.0, 1.0],
    [0.0, 0.0],
    [1.0, 1.0],
])
singular_values = np.linalg.svd(B, compute_uv=False)

print("shape", B.shape)
print("singular_values", singular_values.round(6).tolist())
print("rank", int(np.linalg.matrix_rank(B)))

shape (3, 2)
singular_values [2.0, 0.0]
rank 1

A feature table with rank below its column count deserves inspection. It may contain a deliberate encoding, or two features may duplicate each other and waste compute. Rank exposes the structural fact; domain meaning tells you whether that fact is a bug.

In general, rank is the dimension of a matrix's column space, the span of its columns, and also the dimension of its row space.

Exact rank is a mathematical definition. Floating-point calculations needn't produce a perfect zero for a dependent direction, so np.linalg.matrix_rank treats sufficiently small singular values as zero using a tolerance. Its default accounts for roundoff and matrix dimensions; a measurement-noise threshold may need to be larger. Numerical rank depends on that threshold.^{[3]Reference 3numpy.linalg.matrix_rankhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.matrix_rank.html}

Truncation trades detail for a smaller matrix

Return to full-rank incident matrix A. Keeping both singular directions reconstructs it exactly. What would one direction lose? Keeping only the strongest creates a rank-one approximation:

Among all rank-one matrices, this approximation has the smallest Frobenius reconstruction error: the square root of the sum of squared entry errors.^{[4]Reference 4The approximation of one matrix by another of lower rankhttps://doi.org/10.1007/BF02288367} Equivalently, it minimizes that squared sum. Notice what the guarantee covers: reconstructing A. It doesn't promise that a downstream classifier or retriever keeps its quality.

import numpy as np

A = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
rank_one = U[:, :1] @ np.diag(S[:1]) @ Vt[:1, :]

error = np.linalg.norm(A - rank_one, ord="fro")

print("rank_one")
print(rank_one.round(3))
print("frobenius_error", round(float(error), 4))
print("discarded_singular_value", round(float(S[1]), 4))

rank_one
[[0.5 0.5]
 [0.5 0.5]
 [1.  1. ]]
frobenius_error 1.0
discarded_singular_value 1.0

The discarded direction has strength 1. Look at the top two rows of A₁: a latency-only ticket and a rollback-only ticket both become [0.5, 0.5]. Their contrast is gone. Four entries each have absolute error 0.5, so the Frobenius error is $\sqrt{4\cdot0.5^2} = 1$.

A₁ is still a 3×2 matrix. The smaller representation is its one coordinate per ticket, A @ Vt[:1].T, plus the saved feature direction needed for reconstruction. Rebuilding the full matrix doesn't save storage. This distinction matters when compressing embeddings.

View larger

PCA is SVD after centering each feature

So far, zero meant "no signal." Principal component analysis (PCA) changes the question: which directions describe how tickets vary around their average row?

That question requires centering first. Given a matrix X of tickets by numeric features, subtract each column's average:

For our three tickets, both column means are $2/3$. Subtraction changes t0 from [1, 0] to [1/3, -2/3]: it has more latency and less rollback than the average ticket. The centered rows are [1/3, -2/3], [-2/3, 1/3], and [1/3, 1/3].

Run SVD on those rows and the right singular vectors give principal-component directions. They're also eigenvectors of the sample covariance matrix centered.T @ centered / (n - 1), where n is the number of rows. Its entries measure how feature columns vary together. For this example, centering makes the contrast direction strongest. Check which two tickets differ most from each other before running the code.^{[5]Reference 5sklearn.decomposition.PCAhttps://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html}

import numpy as np

X = np.array([
    [1.0, 0.0],
    [0.0, 1.0],
    [1.0, 1.0],
])
centered = X - X.mean(axis=0)
_, S, Vt = np.linalg.svd(centered, full_matrices=False)
pc1 = Vt[0]
if pc1[0] < 0:
    pc1 = -pc1

retained = S[0] ** 2 / np.sum(S ** 2)

print("column_means", X.mean(axis=0).round(3).tolist())
print("centered_means_are_zero", np.allclose(centered.mean(axis=0), 0))
print("pc1", pc1.round(4).tolist())
print("pc1_variance_fraction", round(float(retained), 4))

column_means [0.667, 0.667]
centered_means_are_zero True
pc1 [0.7071, -0.7071]
pc1_variance_fraction 0.75

Uncentered SVD favored common signal because the rows sit mostly in the positive diagonal direction from the origin. Centered PCA favors contrast because the latency-only and rollback-only tickets are the most separated around their mean. Both results are correct: we asked different questions of the same rows.

PCA finds variance, not relevance. It also centers without automatically rescaling feature units. If one feature is measured in milliseconds and another in seconds, choose the scaling deliberately before fitting. For a new row, reuse the training column means and the fitted directions; don't fit a fresh PCA basis per query.^{[5]Reference 5sklearn.decomposition.PCAhttps://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html}

Avoid materializing a centered sparse matrix

PCA centers every feature mathematically. In a sparse term-document matrix, subtracting a nonzero column mean turns each implicit zero into $-\mu_j$. Explicitly building that centered matrix destroys sparsity. This doesn't make PCA impossible on sparse input: scikit-learn supports sparse inputs with selected solvers, including arpack. The implementation can handle centering without constructing the full centered data table.^{[5]Reference 5sklearn.decomposition.PCAhttps://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html}

A 100,000 by 10,000 matrix has one billion cells. At 1% density, a sparse format stores about 10 million values plus indices. The dense values alone take about 4 GB as float32, or about 8 GB as float64. Those are order-of-magnitude values; index widths and format add their own cost.

If directions relative to the origin are what you want, TruncatedSVD works directly on sparse counts or TF-IDF values (term counts reweighted by how common a term is across documents). Unlike PCA, it doesn't center; this is the usual latent semantic analysis setup.^{[6]Reference 6TruncatedSVDhttps://scikit-learn.org/stable/modules/generated/sklearn.decomposition.TruncatedSVD.html} The choice changes the question, not just memory use. Count the nonzero entries before and after centering a small dense array with roughly 1% nonzero values:

import numpy as np

rng = np.random.default_rng(42)
shape = (1000, 100)  # 100,000 possible entries
X = (rng.random(shape) < 0.01).astype(float)
centered = X - X.mean(axis=0)

print("explicit_values_before", int(np.count_nonzero(X)))
print("explicit_values_after", int(np.count_nonzero(centered)))
print("column_means_all_nonzero", bool(np.all(X.mean(axis=0) != 0)))

assert np.count_nonzero(X) < 1500
assert np.count_nonzero(centered) == 100000

explicit_values_before 1047
explicit_values_after 100000
column_means_all_nonzero True

The example is deliberately small. For a large sparse matrix, don't materialize X or X - mean; use a sparse-aware decomposition path.

Least squares projects a target onto a feature span

SVD also supports the fitting problem behind linear regression. Suppose the three tickets have target severities

The two columns of A are the available latency and rollback features. The first two rows demand weights 1 and 2, but those weights predict 3 for the last row rather than its target 2. No combination matches all three targets. Choose weights $w$ that make the squared residual as small as possible:

This is least squares.^{[1]Reference 1Numerical Linear Algebra.https://people.maths.ox.ac.uk/trefethen/books.html} The fitted vector $Aw$ is the orthogonal projection of y onto the span of A's columns. In other words, choose the point in that feature span closest to the target. For this small example, the hand calculation uses the Gram matrix already computed:

That gives $Aw = [2/3, 5/3, 7/3]$ and residual $y - Aw = [1/3, 1/3, -1/3]$. The residual is perpendicular to both feature columns. Move anywhere else inside their span and the squared error can only stay the same or grow.

NumPy's least-squares solver checks the arithmetic. It finds the weights that minimize the Euclidean residual norm; that also minimizes its square. Inspect the prediction and residual, not just a single loss number:^{[7]Reference 7numpy.linalg.lstsqhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html}

import numpy as np

A = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
y = np.array([1.0, 2.0, 2.0])

weights, _, rank, singular_values = np.linalg.lstsq(A, y, rcond=None)
prediction = A @ weights
residual = y - prediction

print("weights", weights.round(4).tolist())
print("prediction", prediction.round(4).tolist())
print("residual", residual.round(4).tolist())
print("squared_error", round(float(np.sum(residual ** 2)), 4))
print("residual_orthogonal", bool(np.allclose(A.T @ residual, 0)))
print("rank", int(rank), "singular_values", singular_values.round(4).tolist())

weights [0.6667, 1.6667]
prediction [0.6667, 1.6667, 2.3333]
residual [0.3333, 0.3333, -0.3333]
squared_error 0.3333
residual_orthogonal True
rank 2 singular_values [1.7321, 1.0]

Now reuse the duplicate-feature matrix B from the rank example. It maps [1, -1] to zero: increasing one weight and decreasing the other by the same amount leaves every prediction unchanged. Such a direction is in the nullspace, the set of input vectors the matrix maps to zero.^{[1]Reference 1Numerical Linear Algebra.https://people.maths.ox.ac.uk/trefethen/books.html} The closest fitted vector is still unique, but the coefficients that produce it aren't.

For the same target y, the best sum of the two weights is 1.5. NumPy's lstsq chooses the solution with the smallest Euclidean weight norm among equally good fits, splitting that sum into [0.75, 0.75].^{[7]Reference 7numpy.linalg.lstsqhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html} Check a second solution:

import numpy as np

B = np.array([[1.0, 1.0], [0.0, 0.0], [1.0, 1.0]])
y = np.array([1.0, 2.0, 2.0])
weights, _, rank, _ = np.linalg.lstsq(B, y, rcond=None)
alternative = weights + 2.0 * np.array([1.0, -1.0])

print("rank", int(rank))
print("minimum_norm_weights", weights.round(4).tolist())
print("alternative_weights", alternative.round(4).tolist())
print("prediction", (B @ weights).round(4).tolist())
print("same_predictions", bool(np.allclose(B @ weights, B @ alternative)))
print("minimum_norm_is_smaller", bool(np.linalg.norm(weights) < np.linalg.norm(alternative)))
assert np.allclose(B @ np.array([1.0, -1.0]), 0)
assert np.allclose(B @ weights, B @ alternative)

rank 1
minimum_norm_weights [0.75, 0.75]
alternative_weights [2.75, -1.25]
prediction [1.5, 0.0, 1.5]
same_predictions True
minimum_norm_is_smaller True

Both solutions predict [1.5, 0, 1.5]. This table can identify the sum of the feature weights, but can't identify their separate contributions. Interpreting either coefficient as an independently measured feature effect would go beyond the data.

The normal equations make the geometry easy to see, but forming $A^T A$ can square a condition number. For a large or ill-conditioned design matrix, use a direct least-squares or SVD-based routine instead of building the Gram matrix yourself. That numerical warning is next.

Condition number warns about solving through weak directions

SVD also reveals a failure mode. If a matrix stretches one direction strongly and almost erases another, a solve that tries to reverse that transformation must amplify whatever landed in the weak direction.

For example, a diagonal matrix with entries 1 and 0.0001 maps [1, 1] to [1, 0.0001]. Reversing it divides the second output by 0.0001. An output error of only 0.0001 therefore adds 1 to the recovered second coordinate.

The ratio between the largest and smallest stretch measures this imbalance. It's the 2-norm condition number:

when the smallest singular value is nonzero. For an invertible square matrix, a condition number near 1 means directions have comparable scales; it doesn't mean the matrix is close to the identity. A large condition number warns about sensitivity of inverse-like operations. If a singular value is zero, a square matrix has no inverse and its condition number is infinite. Whitening (rescaling projected components toward unit variance) also divides by their strengths, so weak directions need care.^{[8]Reference 8numpy.linalg.condhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.cond.html} The next example isolates that sensitivity:

import numpy as np

A = np.diag([1.0, 0.0001])
x_true = np.array([1.0, 1.0])
y = A @ x_true
y_with_noise = y + np.array([0.0, 0.0001])

x_clean = np.linalg.solve(A, y)
x_noisy = np.linalg.solve(A, y_with_noise)

print("condition_number", float(np.linalg.cond(A)))
print("clean_solution", x_clean.tolist())
print("noisy_solution", x_noisy.tolist())
print("small_output_noise_doubled_second_coordinate", x_noisy[1] == 2.0)

condition_number 10000.0
clean_solution [1.0, 1.0]
noisy_solution [1.0, 2.0]
small_output_noise_doubled_second_coordinate True

Storing vectors in a matrix with a weak direction doesn't automatically make cosine or dot-product retrieval unstable. Amplification appears when a pipeline performs an inverse-like operation, such as a solve or whitening transform, along that weak direction.

Forming AᵀA can square the condition number

A.T @ A made our hand calculation easy, but it isn't the safe general implementation of SVD. For a full-column-rank matrix, forming it squares the 2-norm condition number:

A condition number of 10,000 becomes 100,000,000 in the Gram matrix. The full-column-rank assumption matters: a wide matrix can have a finite ratio of singular values while A.T @ A is singular. Its eigenvalues are the squared singular values, with extra zeros when there are more columns than rows.^{[1]Reference 1Numerical Linear Algebra.https://people.maths.ox.ac.uk/trefethen/books.html} Weak directions can lose useful precision before the decomposition begins. Use a direct SVD implementation in numerical code; keep A.T @ A as a teaching device for the relationship between eigenvectors and singular vectors.

See that precision loss in a nearly duplicated version of B. Change the second row from [0, 0] to [0, 1e-8]. The columns are now independent, and the original matrix still stores their tiny difference. In exact arithmetic, the bottom-right Gram entry would be 2 + 1e-16. In float64, that rounds to 2: the Gram matrix loses the difference and becomes singular. Predict what happens to a solve through that rounded matrix before running the second half:

import numpy as np

A = np.diag([1.0, 0.0001])
G = A.T @ A
print("cond_A", float(np.linalg.cond(A)))
print("cond_AtA", float(np.linalg.cond(G)))
print("ratio", float(np.linalg.cond(G) / np.linalg.cond(A)))

B = np.array([[1.0, 1.0], [0.0, 1e-8], [1.0, 1.0]], dtype=np.float64)
y = B @ np.array([1.0, 1.0])
gram = B.T @ B
direct_weights, _, direct_rank, _ = np.linalg.lstsq(B, y, rcond=None)

print("direct_rank", int(direct_rank))
print("direct_weights", direct_weights.round(4).tolist())
print("rounded_gram", gram.tolist())
try:
    np.linalg.solve(gram, B.T @ y)
except np.linalg.LinAlgError as error:
    print("gram_solve", type(error).__name__ + ": " + str(error))
else:
    raise AssertionError("The rounded Gram matrix should be singular")
assert direct_rank == 2
assert np.allclose(direct_weights, [1.0, 1.0], atol=1e-7, rtol=0)
assert np.array_equal(gram, np.full((2, 2), 2.0))

cond_A 10000.0
cond_AtA 100000000.0
ratio 10000.0
direct_rank 2
direct_weights [1.0, 1.0]
rounded_gram [[2.0, 2.0], [2.0, 2.0]]
gram_solve LinAlgError: Singular matrix

The direct solver recovers [1, 1] for this exact target; the Gram solve fails with LinAlgError: Singular matrix. Using a direct solver avoids this extra loss from forming B.T @ B. It doesn't remove the original problem's sensitivity to noisy observations along the weak direction.

Choosing a cutoff uses squared singular values

For centered data, squared singular values are proportional to the variation explained by each principal direction. For an uncentered count or embedding matrix, they still measure Frobenius energy captured by each direction. In either case, the cumulative ratio gives a reconstruction-based cutoff, where r is the number of singular values. Predict the smallest k that should pass a 95% target from the spectrum below, then check it:

import numpy as np

singular_values = np.array([12.4, 8.7, 3.1, 0.9, 0.3])
fractions = np.cumsum(singular_values ** 2) / np.sum(singular_values ** 2)
target = 0.95
k = int(np.searchsorted(fractions, target) + 1)

print("retained_by_component", fractions.round(4).tolist())
print("smallest_k_for_95_percent", k)
print("retained_at_k", round(float(fractions[k - 1]), 4))

retained_by_component [0.6408, 0.9562, 0.9962, 0.9996, 1.0]
smallest_k_for_95_percent 2
retained_at_k 0.9562

This rule answers a reconstruction question. For a retrieval system, follow it with ranking metrics on held-out queries. Matrix energy isn't a proxy for relevance by itself.

Compress a retrieval matrix and check the ranking

Represent four incident tickets with three term features: latency, rollback, and trace. The query is mostly about a latency trace with a smaller rollback signal. We'll project tickets and query onto the top two right-singular directions, then compare rankings before and after compression with cosine similarity.

Cosine similarity divides each dot product by both vector lengths, so it compares direction rather than raw magnitude. Before running the fixture, predict whether dropping only the weakest direction should change the top ticket.

import numpy as np

tickets = np.array([
    [3.0, 0.0, 1.0],  # latency trace
    [0.0, 3.0, 1.0],  # rollback trace
    [2.0, 1.0, 0.5],  # latency with rollback mention
    [1.0, 2.0, 1.5],  # rollback with latency mention
])
query = np.array([2.5, 0.5, 1.0])

def cosine_scores(matrix: np.ndarray, vector: np.ndarray) -> np.ndarray:
    row_norms = np.linalg.norm(matrix, axis=1)
    vector_norm = np.linalg.norm(vector)
    if np.any(row_norms == 0) or vector_norm == 0:
        raise ValueError("cosine similarity requires non-zero vectors")
    return (matrix @ vector) / (row_norms * vector_norm)

_, singular_values, Vt = np.linalg.svd(tickets, full_matrices=False)
basis = Vt[:2].T
retained_energy = np.sum(singular_values[:2] ** 2) / np.sum(singular_values ** 2)

original_scores = cosine_scores(tickets, query)
compressed_scores = cosine_scores(tickets @ basis, query @ basis)
original_order = np.argsort(original_scores)[::-1]
compressed_order = np.argsort(compressed_scores)[::-1]

print("singular_values", singular_values.round(4).tolist())
print("retained_energy", round(float(retained_energy), 4))
print("original_order", original_order.tolist())
print("compressed_order", compressed_order.tolist())
print("top_ticket_preserved", bool(original_order[0] == compressed_order[0]))
print("max_score_change", round(float(np.max(np.abs(original_scores - compressed_scores))), 4))

try:
    cosine_scores(tickets, np.zeros(3))
except ValueError as error:
    print("zero_vector_guard", str(error))

singular_values [4.7011, 3.1677, 0.6044]
retained_energy 0.9888
original_order [0, 2, 3, 1]
compressed_order [0, 2, 3, 1]
top_ticket_preserved True
max_score_change 0.0221
zero_vector_guard cosine similarity requires non-zero vectors

For this tiny fixture, dropping the tail direction changes scores a little and keeps the same order. The explicit norm guard matters too: cosine similarity is undefined for a zero-length query or ticket vector. A serving pipeline should reject, quarantine, or replace that vector instead of returning NaN scores.

The fixture is a sanity check, not a credible evaluation. On a real retrieval set, measure Recall@k (the fraction of relevant tickets recovered in the top k returned candidates) and inspect failures before accepting reduced dimensions.

Low-rank factors in LoRA and embedding retrieval

The same rank question returns once a model has to adapt or retrieve. The arithmetic above is the part you can check by hand. Run that same question on the task's own data before you trust the result there.

LoRA doesn't compute the SVD of a full trained update during fine-tuning. It learns thin factors directly. SVD gives you the language for "low rank"; later training lessons explain how those factors are optimized.

Practice with measured answers

Give each direction a name and a strength, and say what a cutoff discards. Use the same fixture to test each idea under a small change. Start with the symmetric matrix, then stress the interpretation rather than the arithmetic:

1. Start with the symmetric matrix $C = \begin{bmatrix}3 & 1 \\ 1 & 3\end{bmatrix}$: compute $C^T C$, its two eigenvalues, singular values, and condition number.
2. In the retrieval exercise, keep only one singular direction. Did the first result stay the same? How large did the maximum score change become?
3. Add a large constant offset [100, 100] to every row of the PCA example. Compare SVD before and after centering. Which version describes differences between rows?
4. Fit y = [1, 2, 2] with the columns of A using least squares. What are the weights, prediction, and residual?

Check your reasoning against these answers:

1. $C^T C = \begin{bmatrix}10 & 6 \\ 6 & 10\end{bmatrix}$, eigenvalues are 16 and 4, singular values are 4 and 2, and $\kappa_2(C) = 2$.
2. In one dimension, nonzero scalars with the same sign have cosine similarity 1. All four tickets and this query project to the same sign, so all four scores are 1 up to roundoff. The order of those ties isn't evidence of relevance and can vary with numerical precision or sorting behavior. The maximum score change is approximately 0.7113.
3. Uncentered SVD is pulled toward [1, 1] by the offset. Centering removes it and recovers the contrast direction [1, -1]/sqrt(2), up to sign, just as before.
4. The weights come out to [2/3, 5/3], yielding prediction [2/3, 5/3, 7/3] with residual [1/3, 1/3, -1/3]. That residual is orthogonal to both columns of A.

Our three-ticket matrix leaves one concrete warning: its weaker contrast direction is exactly what distinguishes a latency-only ticket from a rollback-only ticket. SVD tells you how much reconstruction error discarding it costs. Only the downstream task can tell you whether losing that distinction is acceptable.