Linear Algebra
you can see.

Linear algebra is the language of neural networks

Imagine an online store that estimates delivery time from distance, traffic, and number of stops. Three numbers describe one order. A neural network takes those numbers and transforms them into a prediction. The transformation is linear algebra.

Unlike a chapter that simply lists operations, this chapter shows what gets multiplied by what, how the geometry changes, and how the same calculations appear in production model layers.

Scalars, vectors, coordinates, and norms

What is a vector?

A vector is an ordered collection of numbers. Geometrically, a 2D vector can be drawn as an arrow. In ML, its entries are feature values, often with different units and meanings.

The Euclidean (L2) norm measures vector length. It comes from the Pythagorean theorem:

Unit vectors and normalization

Dividing a nonzero vector by its length produces a unit vector, pointing in the same direction with length one. This is useful when direction matters more than magnitude—for example, comparing text embeddings.

import math
import numpy as np

x = [3.0, 4.0]

# WHY: squaring and summing follows the Pythagorean distance rule.
squared_length = sum(component ** 2 for component in x)
length = math.sqrt(squared_length)
print("Manual L2 norm:", length)  # 5.0

# WHY: NumPy's norm works directly on vectors and matrices.
x_np = np.array(x, dtype=np.float64)
print("NumPy L2 norm:", np.linalg.norm(x_np))

# WHY: avoid division by zero for a zero vector.
unit = x_np / length if length > 0 else np.zeros_like(x_np)
print("Unit vector:", unit)  # [0.6 0.8]

Other useful norms: L1 sums absolute values; L2 is geometric length; squared L2 often appears in loss functions. Regularization can penalize large parameter norms.

Dot product: similarity and weighted scores

Suppose a recommendation system gives three features importance scores. The dot product combines them into one value. We disassemble it into paired multiplications, then add.

Connection to an angle

For nonzero vectors, the dot product also measures directional alignment:

If vectors point the same way, cosine similarity is +1; if perpendicular, 0; if opposite, −1. This becomes important in embedding search and semantic retrieval.

import numpy as np
import torch

features = [5.0, 0.8, 3.0]
weights = [2.0, 4.0, 3.0]

# WHY: zip aligns each feature with the corresponding weight.
manual = sum(x * w for x, w in zip(features, weights))
print("Manual:", manual)  # 22.2

# WHY: @ computes the dot product for 1D NumPy arrays.
x_np = np.array(features)
w_np = np.array(weights)
print("NumPy:", x_np @ w_np)

# WHY: PyTorch tensors support the same mathematical operation.
x_t = torch.tensor(features)
w_t = torch.tensor(weights)
print("PyTorch:", torch.dot(x_t, w_t).item())

Matrices: organize many features and examples

A matrix is a rectangular arrangement of numbers. Each row can describe a separate delivery and each column a feature. If there are 4 deliveries and 3 features, X.shape == (4,3).

Matrix–vector multiplication: one prediction for every row

Shapes tell the story: (4,3) @ (3,) → (4,). This is just four dot products, one per row. Adding a scalar bias applies the same baseline to all four results.

import numpy as np

# WHY: one row per example; one column per input feature.
X = np.array([[5, .8, 3], [2, .3, 1],
              [8, .9, 4], [3, .5, 2]])

# WHY: one weight for each feature.
w = np.array([2.0, 4.0, 3.0])
b = 5.0  # WHY: a shared baseline prediction (intercept).

# Disassemble the batch computation into individual dot products.
manual = np.array([sum(row[j] * w[j] for j in range(3)) + b
                   for row in X])

# Reassemble: one optimized matrix-vector operation.
vectorized = X @ w + b
assert np.allclose(manual, vectorized)
print(vectorized)  # [27.2 13.2 36.6 19. ]

Transpose

The transpose switches rows and columns: a (4,3) matrix becomes (3,4). It is crucial for gradient calculations, batch data, and the way some frameworks store weights.

Trap: transpose does not mean inverse. Rotating the index layout and undoing a transformation are different operations.

Matrix multiplication: disassemble, compute, reassemble

A fully connected neural network with 3 input features and 2 output neurons needs a weight matrix with shape (3,2). Each output neuron takes its own weighted sum of the input features.

Unlike element-wise multiplication (X * W), matrix multiplication combines a row with a column. The inner dimensions must agree; otherwise the operation is undefined.

import numpy as np

X = [[1, 2, 3], [4, 5, 6]]  # WHY: batch of 2, each with 3 features.
W = [[1, 2], [0, -1], [3, 1]]  # WHY: map 3 inputs to 2 outputs.

# WHY: allocate one number per sample and output neuron.
Z = [[0] * 2 for _ in range(2)]

for i in range(2):          # choose sample / output row
    for j in range(2):      # choose output neuron / column
        for k in range(3):  # pair input feature with its weight
            Z[i][j] += X[i][k] * W[k][j]

print("Loop result:", Z)  # [[10, 3], [22, 9]]

# WHY: @ replaces the three Python loops with matrix multiplication.
fast = np.array(X) @ np.array(W)
assert np.array_equal(fast, np.array(Z))
print("NumPy result:\n", fast)

Linear transformations: watch space change

Matrix multiplication is also geometry. A 2×2 matrix transforms each point in a 2D plane: stretching, squashing, shearing or reflecting it. A transformation is linear when it preserves vector addition and scalar multiplication.

For a 2D matrix, absolute determinant tells how area scales; its sign tells whether orientation flips. A zero determinant collapses the plane into a line or point, losing information.

import numpy as np

# WHY: each row represents an (x, y) coordinate.
points = np.array([[1., 0.], [0., 1.], [1., 1.]])
A = np.array([[1.5, 0.5], [0.0, 1.0]])

# WHY: row-wise point storage uses A.T to apply y = A @ x.
transformed = points @ A.T

print("Transformed points:\n", transformed)
print("Determinant:", np.linalg.det(A))  # 1.5

# WHY: rank indicates how many independent directions remain.
print("Matrix rank:", np.linalg.matrix_rank(A))

Rank, span, independence, inverse, and solving systems

Imagine two sensors that always report the same information. Although your dataset has two columns, it effectively contains one independent signal. Rank captures the number of independent directions represented by a matrix.

Columns are linearly dependent if one can be made from the others. The span of columns is the set of all linear combinations. Full rank means no independent direction is redundant, relative to the smaller dimension.

Inverse is not transpose

A square matrix is invertible if it has full rank (equivalently, nonzero determinant). But practical ML code generally solves linear systems using dedicated solvers instead of computing an explicit inverse.

import numpy as np

A = np.array([[2., 1.], [1., 3.]])
b = np.array([5., 7.])

# WHY: solve Ax=b directly rather than calculating inv(A) @ b.
x = np.linalg.solve(A, b)
print("Solution:", x)  # [1.6 1.8]
assert np.allclose(A @ x, b)

# WHY: a redundant matrix has fewer independent directions.
redundant = np.array([[1., 2.], [2., 4.]])
print("Rank:", np.linalg.matrix_rank(redundant))  # 1

# WHY: use least squares if data has more equations than unknowns.
X = np.array([[1., 1.], [1., 2.], [1., 3.]])
y = np.array([2.1, 2.9, 4.2])
weights, residuals, rank, singular_values = np.linalg.lstsq(
    X, y, rcond=None
)
print("Least-squares weights:", weights)

Practical note: conditioning also matters. Even an invertible matrix can be difficult to solve accurately if its columns are nearly dependent.

Eigenvalues and eigenvectors: directions a matrix preserves

Most arrows change both length and direction when transformed by a matrix. But special directions, called eigenvectors, only get scaled (and possibly reversed).

For example, a diagonal scaling matrix stretches horizontal vectors twice and shrinks vertical vectors by half: the coordinate axes are eigenvector directions with eigenvalues 2 and 0.5.

Eigenvectors are especially useful for understanding symmetric covariance matrices and principal component analysis. Not all real matrices have real eigenvalues or a full set of independent eigenvectors.

import numpy as np

# WHY: symmetric matrices have real eigenvalues and orthogonal eigenvectors.
A = np.array([[2., 0.], [0., 0.5]])

# WHY: eigh is designed for real symmetric / Hermitian matrices.
values, vectors = np.linalg.eigh(A)

for i in range(len(values)):
    v = vectors[:, i]  # WHY: eigenvectors are returned as COLUMNS.
    assert np.allclose(A @ v, values[i] * v)
    print("Eigenvalue:", values[i], "Eigenvector:", v)

Singular value decomposition: rotate, stretch, rotate

Eigenvectors help explain some transformations, but singular value decomposition (SVD) works for any real matrix, including rectangular ones. Think of it as three coordinated steps: change the input basis, stretch along independent axes, then place the result in output space.

Why does this matter in AI? SVD underpins low-rank approximations, compression, dimensionality reduction, and understanding feature spaces. Low-rank structure also motivates parameter-efficient model adaptations, though those methods have additional engineering details.

import numpy as np

A = np.array([[3., 1.], [1., 3.], [2., 2.]])

# WHY: reduced SVD returns compact factors for a rectangular matrix.
U, singular_values, Vt = np.linalg.svd(A, full_matrices=False)

# WHY: multiply each column of U by its singular value.
reconstructed = (U * singular_values) @ Vt
assert np.allclose(A, reconstructed)

# WHY: k=1 keeps only the strongest singular direction.
k = 1
low_rank = (U[:, :k] * singular_values[:k]) @ Vt[:k, :]

print("Singular values:", singular_values)
print("Rank-one approximation:\n", low_rank)
print("Reconstruction error:", np.linalg.norm(A - low_rank))

When singular values decay quickly, relatively few directions may approximate the original matrix well. The approximation error is not generally zero if components are discarded.

Reassemble: the exact linear algebra inside deep learning

Let us now reconstruct a batch neural-network layer. A single row of X holds the features of one delivery; every column of W corresponds to a neuron. The matrix product computes all neuron weighted sums for all examples at once.

Example shapes: X: (4,3), W: (3,2), b: (2,), Z: (4,2). The bias broadcasts across the four examples.

import numpy as np

X = np.array([[5., .8, 3.], [2., .3, 1.],
              [8., .9, 4.], [3., .5, 2.]])
W = np.array([[2., -1.], [4., 3.], [3., 2.]])
b = np.array([5., 1.])

# WHY: @ calculates all weighted sums using row-column dot products.
Z = X @ W + b

# WHY: ReLU keeps positive activations and clips negatives to zero.
H = np.maximum(0, Z)

print("Linear output shape:", Z.shape)  # (4, 2)
print("Activated output shape:", H.shape)  # (4, 2)
print("Activations:\n", H)

import torch
from torch import nn

X = torch.tensor([[5., .8, 3.], [2., .3, 1.],
                  [8., .9, 4.], [3., .5, 2.]])
W = torch.tensor([[2., -1.], [4., 3.], [3., 2.]])
b = torch.tensor([5., 1.])

layer = nn.Linear(3, 2)  # WHY: map three features to two neuron outputs.

with torch.no_grad():  # WHY: assignment is setup, not an operation to differentiate.
    # WHY: PyTorch stores nn.Linear weights as (out_features, in_features).
    layer.weight.copy_(W.T)
    layer.bias.copy_(b)

# WHY: nn.Linear computes X @ layer.weight.T + layer.bias.
Z = layer(X)
H = torch.relu(Z)

# WHY: verify our framework and manual formula agree.
assert torch.allclose(Z, X @ W + b)
print(H)

Where these ideas appear later

Interview questions: reason, do not memorize

Try to answer each question aloud before expanding the model answer.

Exercises, self-checks, and official references

Official resources

• NumPy: matrix multiplication
• NumPy: linear algebra reference
• NumPy: singular value decomposition
• PyTorch: nn.Linear
• PyTorch: torch.linalg
• Dive into Deep Learning: linear algebra

Next: Chapter 04 · Calculus and chain rule. We will use vectors, dot products, and transposes to understand gradients and backpropagation.