Linear Algebra for Deep Learning

The dot product, the L2 norm, unit vectors, and cosine similarity are built up from scratch in Linear Algebra for Machine Learning in the companion series; this article assumes them and only recalls what it needs. What is new here is the tensor vocabulary, the difference between the two kinds of multiplication, the shape rule for applying one weight matrix to a batch of inputs, and why transposes appear in backpropagation.

Dense layers combine weighted inputs, add biases, and pass their outputs to activations. Following these calculations requires distinguishing matrix multiplication from elementwise operations and keeping track of what each array axis represents. We will use those ideas to calculate a layer output and understand the transposes in its gradients. Reductions such as sum(axis=...) are introduced in Python and NumPy Basics.

Scalars, vectors, matrices, tensors

A scalar is one number. A vector is an ordered list, written \(x\in\mathbb{R}^{n}\) and treated as a column by convention. A matrix is a rectangle of numbers, \(W\in\mathbb{R}^{m\times n}\) with \(m\) rows and \(n\) columns. A tensor, as the word is used in this series and in the frameworks, is an array with any number of axes: a scalar is a 0-axis tensor, a vector a 1-axis tensor, a matrix a 2-axis tensor, and a batch of RGB images a 4-axis tensor whose axes in the layout used here are batch, channels, height, and width. Other layouts place channels last.

The notation \(\mathbb{R}^{m\times n}\) means an \(m\)-by-\(n\) grid of real numbers. In code, give those axes a meaning as well as a size: an axis might index features, examples, or image channels.

import numpy as np

s = 3.0                                  # scalar
x = np.array([1.0, 2.0, 3.0])            # vector,  shape (3,)
W = np.array([[1.0, 2.0, 3.0],
              [4.0, 5.0, 6.0]])          # matrix,  shape (2, 3)
T = np.zeros((32, 3, 224, 224))          # tensor,  (batch, channels, height, width)

print(x.shape, W.shape, T.shape)         # (3,) (2, 3) (32, 3, 224, 224)
print(W.ndim, T.ndim)                    # 2 4

A one-dimensional NumPy array has shape (3,). Transposing it leaves that shape unchanged: there is no second axis to exchange with the first. Use reshape(3, 1) to represent a column explicitly; its transpose then has shape (1, 3).

print(x.T.shape)
# (3,)
x_col = x.reshape(3, 1)
print(x_col.shape, x_col.T.shape)
# (3, 1) (1, 3)

The two ways to multiply

Elementwise multiplication (the Hadamard product, \(A\odot B\)) multiplies matching positions. The mathematical operation is defined for two operands of the same shape, and the result has that shape too. In code it is * — but NumPy’s * also accepts broadcast-compatible shapes, so multiplying arrays with shapes (2,1) and (1,3) returns a (2,3) array. The vectorization article covers those rules; the affine layer later in this article already relies on them.

Matrix multiplication (\(AB\)) combines rows and columns: entry \((i,j)\) of the result is the dot product of row \(i\) of \(A\) with column \(j\) of \(B\). In code it is @.

\[C_{ij}=\sum_{k=1}^{n}A_{ik}B_{kj}\]

A = np.array([[1.0, 2.0],
              [3.0, 4.0]])
B = np.array([[5.0, 6.0],
              [7.0, 8.0]])

print(A * B)      # elementwise
# [[ 5. 12.]
#  [21. 32.]]
print(A @ B)      # matrix product
# [[19. 22.]
#  [43. 50.]]

The top-left entry of the matrix product is \(19 = 1\cdot5 + 2\cdot7\): multiply the first row of \(A\) by the first column of \(B\), then add the two products.

The shape rule

For \(AB\) to exist, the inner dimensions must match:

\[(m\times \mathbf{n})\;(\mathbf{n}\times p)\;\longrightarrow\;(m\times p)\]

The sum runs over the shared \(n\) entries, leaving one result for each row of the first matrix and column of the second. Writing these shapes beside the code helps check whether the intended inputs and outputs line up.

W = np.random.randn(4, 3)       # (4, 3)
x = np.random.randn(3, 1)       # (3, 1)
print((W @ x).shape)            # (4, 1) # sum over the 3 input features

X = np.random.randn(3, 100)     # 100 examples, each of size 3
print((W @ X).shape)            # (4, 100) # one matmul handles all 100 examples

try:
    x @ W                       # (3,1) @ (4,3): inner 1 vs 4
except ValueError:
    print("ValueError: inner dimensions 1 and 4 do not match")
# ValueError: inner dimensions 1 and 4 do not match

The weighted sum for one input is \(Wx\); applied to a hundred inputs stored as columns it is \(WX\), with no loop. This series stores one example per column, so \(X\) is (features, batch). The companion machine learning series stores one example per row, and PyTorch’s nn.Linear expects features on the last axis: (batch, features) for a two-dimensional batch, or more generally (..., features) with any number of leading axes. For a row-based batch with the same weight matrix, the corresponding product is X_rows @ W.T.

Matrix multiplication is associative, \((AB)C=A(BC)\), and distributive, but not commutative: \(AB\ne BA\) in general, and \(BA\) may not even be a valid shape. Order matters everywhere in this series.

Transpose

The transpose \(A^T\) exchanges rows and columns, turning shape \((m,n)\) into \((n,m)\). To see why it appears in backpropagation, let \(Z=WX\) and let \(L\) be a scalar loss that depends on this layer through \(Z\). Write \(dZ=\partial L/\partial Z\) for the gradient supplied by later computations: each entry tells us how the loss changes when the corresponding entry of \(Z\) changes. The calculus article introduces the chain rule used here.

For one output entry, \(Z_{ij}=\sum_k W_{ik}X_{kj}\). Changing \(W_{ik}\) affects row \(i\) of the output; its contribution to output \(Z_{ij}\) is scaled by \(X_{kj}\). Summing the loss contributions over examples gives:

\[\frac{\partial L}{\partial W_{ik}}=\sum_j dZ_{ij}X_{kj}=(dZ\,X^T)_{ik}.\]

The transpose puts those entries of \(X\) in the column needed for the dot product. Similarly, changing \(X_{kj}\) affects output column \(j\), and summing over output units gives \(\partial L/\partial X_{kj}=\sum_i W_{ik}dZ_{ij}\). Together these give \(dW=dZ\,X^T\) and \(dX=W^T dZ\). Shapes provide a useful check, but cannot establish the derivatives: multiplying either expression by 2 would preserve its shape and generally change its value.

The random dZ below stands in for a gradient supplied by later computations. This example checks the output shapes; it does not compute a loss or obtain its gradient.

W = np.random.randn(4, 3)
X = np.random.randn(3, 100)
Z = W @ X                             # the forward pass
dZ = np.random.randn(4, 100)          # dL/dZ, same shape as Z

print((W.T @ dZ).shape)               # (3, 100) # shape of X
print((dZ @ X.T).shape)               # (4, 3) # shape of W

Two useful identities are \((A^T)^T=A\) and \((AB)^T=B^T A^T\). Transposing a product reverses the order of its factors.

Where dot products and norms show up here

The dot product \(x^{T}y\), its geometric reading \(\lVert x\rVert\lVert y\rVert\cos\theta\), the L2 norm, unit vectors, and cosine similarity are developed from the ground up in Linear Algebra for Machine Learning. This section only notes where they turn up in the rest of this series.

Each output of a dense layer contains a dot product between a weight row and an input vector. When a transformer computes \(QK^{T}\), each row of \(Q\) holds one query vector and each row of \(K\) holds one key vector. This row-based layout differs from the column-based input layout used above. The product collects a dot product between every query and every key. Each score reflects both the direction and the magnitude of the two vectors; scaled dot-product attention then divides by \(\sqrt{d_k}\), where \(d_k\) is the number of components in a query or key vector, applies a softmax across keys to produce nonnegative weights that sum to one for each query, and uses those weights to average the value vectors. Transformers from Scratch develops the attention calculation.

An L2 regularization penalty uses the squared L2 norm to discourage large weights. Norm-based gradient clipping uses the norm of the gradient to limit its magnitude. Dividing a nonzero vector by its L2 norm gives a unit vector, and the dot product of two unit vectors equals their cosine similarity — which is how the embedding articles compare representations when magnitude should not count. The L1 norm \(\lVert x\rVert_1=\sum_i|x_i|\) sums absolute values; using it as a penalty can encourage sparse weights.

What you can safely skip

Determinants, matrix inverses, and hand-solving \(Ax=b\) are not prerequisites for the layer calculations here. The affine layers and elementwise activations used in these introductory examples can be evaluated and differentiated without forming a matrix inverse.

Eigenvalues help describe repeated application of a fixed square matrix, which provides a useful linear model of recurrence. In a nonlinear recurrent network, activation derivatives also enter the Jacobian at each time step. A Jacobian is the matrix of derivatives of one vector output with respect to its vector input. Both recurrent and feedforward gradients can involve products of different Jacobians; the eigenvalues of individual factors alone do not determine finite-depth amplification. Singular values describe how much a linear map can stretch or shrink a vector and appear in principal component analysis and spectral normalization. These concepts are introduced where they are needed.

The array operations assumed here are covered in Python and NumPy Basics.

Checking shapes at layer boundaries

An affine layer adds a learned bias to the linear map: \(Z=WX+b\). With one example per column, \(Z_{ij}=\sum_k W_{ik}X_{kj}+b_i\). Each output unit has its own bias, shared across the examples. A bias stored as \((n_{\text{out}},1)\) broadcasts across the columns; it can make the output nonzero even when the input is zero.

W_small = np.array([[1., 2., 0.], [0., -1., 1.]])
X_small = np.array([[1., 2.], [3., 0.], [2., 4.]])
b_small = np.array([[0.5], [-1.]])
print(W_small @ X_small + b_small)
# [[ 7.5  2.5]
#  [-2.   3. ]]

The first column is the output for input \((1,3,2)^T\): its first entry is \(1\cdot1+2\cdot3+0\cdot2+0.5=7.5\), and its second is \(0\cdot1-1\cdot3+1\cdot2-1=-2\). The same weights and biases produce the second column from the second input. An activation, such as ReLU applied elementwise, would be a separate operation after this affine calculation.

Annotate shapes in comments as you write, and check them at function boundaries. These checks can catch incompatible dimensions and unintended broadcasting, including cases where NumPy returns an array without raising an error.

def layer_forward(W, X, b):
    """W: (n_out, n_in) | X: (n_in, m) | b: (n_out, 1) -> (n_out, m)"""
    assert W.ndim == 2 and X.ndim == 2, "W and X must both be 2-D"
    assert W.shape[1] == X.shape[0], f"inner dim {W.shape[1]} vs {X.shape[0]}"
    assert b.shape == (W.shape[0], 1), f"b must be {(W.shape[0], 1)}, got {b.shape}"
    Z = W @ X + b                       # b broadcasts across the m columns
    assert Z.shape == (W.shape[0], X.shape[1])
    return Z

print(layer_forward(np.random.randn(4, 3), np.random.randn(3, 10),
                    np.zeros((4, 1))).shape)      # (4, 10)

The weight and input gradients derived above also apply when this bias is added. For the bias itself, changing \(b_i\) changes every \(Z_{ij}\) in its row by the same amount, so \(db_i=\sum_j dZ_{ij}\). In NumPy, db = dZ.sum(axis=1, keepdims=True) adds the contributions across examples and keeps shape (n_out, 1). If the loss averages over examples, that scaling is already included in dZ; do not divide by the batch size a second time.

The three input checks serve different purposes. Without the check on b.shape, a bias of shape (10,) passes the other checks and still returns a (4,10) result — but it has broadcast across the wrong axis, adding a per-example offset instead of a per-neuron bias, and the final check on Z.shape does not notice. The ndim check comes first so that the later .shape[1] lookups are safe.

The first network built from these pieces is in Neural Networks for Beginners: From Architecture to Forward Propagation.