age sbp bmi
p1 64 138 27.4
p2 51 122 31.9
p3 73 155 24.1
p4 45 118 29.65 Matrices
A vector is a list of numbers. A matrix is a grid of them. That sounds like a small step, and the arithmetic in the first half of this chapter is genuinely mechanical.
The second half is where it stops being bookkeeping. A matrix is not really a grid of numbers at all — it is a function that takes a vector and returns a vector, and almost everything interesting about matrices follows from asking what that function does to space.
5.1 What a matrix is
A matrix is a rectangular array of numbers, arranged in rows and columns.
\[ \mathbf{A} = \begin{bmatrix} 2 & 1 \\ 0 & 3 \end{bmatrix} \]
A matrix with \(m\) rows and \(n\) columns is an \(m \times n\) matrix — rows first, always. We say it has shape or dimensions \(m \times n\), and write \(\mathbf{A} \in \mathbb{R}^{m \times n}\). The matrix above is \(2 \times 2\).
Individual entries get two subscripts, row then column: \(a_{ij}\) is the entry in row \(i\), column \(j\). So \(a_{12} = 1\) and \(a_{21} = 0\). Neither is bold, because each is a single number.
\[ \mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix} \]
A vector is the special case with one column: an \(n \times 1\) matrix. That is why column vectors are the default — it makes vectors and matrices the same kind of object rather than two kinds that need translating between.
5.2 A matrix as a table of data
The most common matrix in data science is the design matrix: one row per observation, one column per feature.
R fills a matrix column by column by default. Writing the numbers out in reading order and passing byrow = TRUE is almost always clearer — and forgetting it is one of the most common sources of a silently transposed matrix.
[,1] [,2]
[1,] 2 0
[2,] 1 3 [,1] [,2]
[1,] 2 1
[2,] 0 3Rows and columns mean different things here and are not interchangeable: a row is a patient, a column is a measurement. Most of the confusion people have with matrix orientation comes from losing track of which is which.
NoteIn machine learning
The convention is nearly universal: \(\mathbf{X}\) is \(n \times p\) with \(n\) observations and \(p\) features. Every model you fit — lm, glm, a random forest, a neural network — expects that orientation, and reports coefficients in the order of the columns.
5.3 Shape, indexing, and slicing
Index with [row, column]. Leaving one slot empty takes everything along that dimension.
A[1, 2] # single entry: row 1, column 2[1] 1A[1, ] # row 1, as a vector[1] 2 1A[, 2] # column 2, as a vector[1] 1 35.4 Special matrices
A handful of shapes come up constantly and each has a name worth knowing.
Square: same number of rows and columns. Only square matrices have determinants, inverses, or eigenvalues.
Identity, written \(\mathbf{I}\): square, ones on the diagonal, zeros elsewhere. It is the matrix that does nothing — the matrix equivalent of multiplying by 1.
\[ \mathbf{I}_3 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \]
diag(3) [,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1Diagonal: zero everywhere off the diagonal. Multiplying by one just rescales each coordinate independently, which is why diagonal matrices are the easy case of almost every problem.
Symmetric: \(\mathbf{A} = \mathbf{A}^\top\), so \(a_{ij} = a_{ji}\) — the matrix is its own mirror image across the diagonal. Covariance matrices, correlation matrices, and distance matrices are all symmetric, and symmetric matrices have properties so much better than general ones that Section 8.7 is devoted to them.
[,1] [,2] [,3]
[1,] 4 2 1
[2,] 2 5 3
[3,] 1 3 6[1] TRUEdraw_heatmap(S, show_values = TRUE)Triangular: all zeros above (lower triangular) or below (upper triangular) the diagonal. Systems involving them can be solved by substitution without any elimination work, which is the entire point of the decompositions in Chapter 9.
[,1] [,2] [,3]
[1,] 2 1 4
[2,] 0 3 5
[3,] 0 0 1Orthogonal: a square matrix whose columns are mutually orthogonal unit vectors, so \(\mathbf{Q}^\top\mathbf{Q} = \mathbf{I}\). These represent rotations and reflections — transformations that move things without stretching them.
5.5 Transpose
The transpose flips a matrix across its diagonal, turning rows into columns:
\[ (\mathbf{A}^\top)_{ij} = a_{ji} \]
An \(m \times n\) matrix transposes to an \(n \times m\) matrix.
Worked example.
\[ \mathbf{A} = \begin{bmatrix} 2 & 1 \\ 0 & 3\end{bmatrix} \qquad \mathbf{A}^\top = \begin{bmatrix} 2 & 0 \\ 1 & 3\end{bmatrix} \]
t(A) [,1] [,2]
[1,] 2 0
[2,] 1 3Two facts to keep:
\[ (\mathbf{A}^\top)^\top = \mathbf{A} \qquad\text{and}\qquad (\mathbf{A}\mathbf{B})^\top = \mathbf{B}^\top\mathbf{A}^\top \]
The second one surprises people: transposing a product reverses the order. It has to, or the shapes would not line up. We verify it numerically in Section 5.9.
5.6 Addition and scalar multiplication
Both are elementwise, exactly as for vectors, and addition requires identical shapes.
\[ (\mathbf{A} + \mathbf{B})_{ij} = a_{ij} + b_{ij} \qquad (c\mathbf{A})_{ij} = c\,a_{ij} \]
Worked example. With \(\mathbf{B} = \begin{bmatrix}1 & 4\\2 & 1\end{bmatrix}\):
\[ \mathbf{A} + \mathbf{B} = \begin{bmatrix} 2+1 & 1+4 \\ 0+2 & 3+1 \end{bmatrix} = \begin{bmatrix} 3 & 5 \\ 2 & 4 \end{bmatrix} \]
[,1] [,2]
[1,] 3 5
[2,] 2 42 * A [,1] [,2]
[1,] 4 2
[2,] 0 6Nothing surprising happens here. All the surprise is in multiplication.
5.7 Matrix–vector multiplication: two views
For \(\mathbf{A} \in \mathbb{R}^{m \times n}\) and \(\mathbf{x} \in \mathbb{R}^n\), the product \(\mathbf{A}\mathbf{x}\) is a vector in \(\mathbb{R}^m\):
\[ (\mathbf{A}\mathbf{x})_i = \sum_{j=1}^{n} a_{ij}x_j \tag{5.1}\]
The inner dimension has to match: \(\mathbf{A}\) has \(n\) columns, so \(\mathbf{x}\) needs \(n\) elements. There are two ways to read Equation 5.1, and you need both.
5.7.1 The row view
Each element of the output is the dot product of one row of \(\mathbf{A}\) with \(\mathbf{x}\).
Worked example. With \(\mathbf{x} = \begin{bmatrix}1\\2\end{bmatrix}\):
\[ \mathbf{A}\mathbf{x} = \begin{bmatrix} (2,1)\cdot(1,2) \\ (0,3)\cdot(1,2) \end{bmatrix} = \begin{bmatrix} 2 + 2 \\ 0 + 6 \end{bmatrix} = \begin{bmatrix} 4 \\ 6 \end{bmatrix} \]
This is the view to use when computing by hand, and the one that makes the shape rule obvious.
5.7.2 The column view
Each column of \(\mathbf{A}\) is a vector, and \(\mathbf{A}\mathbf{x}\) is the linear combination of those columns with coefficients taken from \(\mathbf{x}\):
\[ \mathbf{A}\mathbf{x} = x_1\mathbf{a}_1 + x_2\mathbf{a}_2 + \cdots + x_n\mathbf{a}_n \tag{5.2}\]
where \(\mathbf{a}_j\) is the \(j\)-th column. Same example:
\[ \mathbf{A}\mathbf{x} = 1\begin{bmatrix}2\\0\end{bmatrix} + 2\begin{bmatrix}1\\3\end{bmatrix} = \begin{bmatrix}2\\0\end{bmatrix} + \begin{bmatrix}2\\6\end{bmatrix} = \begin{bmatrix}4\\6\end{bmatrix} \]
[,1]
[1,] 4
[2,] 6A[, 1] * x[1] + A[, 2] * x[2] # the same thing, column view[1] 4 6Equation 5.2 is the more important of the two. It says that the set of vectors \(\mathbf{A}\) can possibly produce is exactly the set of linear combinations of its columns — a fact that decides when \(\mathbf{A}\mathbf{x} = \mathbf{b}\) has a solution (Section 6.5) and defines the column space (Section 7.9).
5.8 Matrix–matrix multiplication
Multiplying \(\mathbf{A} \in \mathbb{R}^{m \times n}\) by \(\mathbf{B} \in \mathbb{R}^{n \times p}\) gives \(\mathbf{C} \in \mathbb{R}^{m \times p}\), where
\[ c_{ij} = \sum_{k=1}^{n} a_{ik}b_{kj} \tag{5.3}\]
In words: entry \((i,j)\) of the product is the dot product of row \(i\) of \(\mathbf{A}\) with column \(j\) of \(\mathbf{B}\). Equivalently, column \(j\) of the product is \(\mathbf{A}\) applied to column \(j\) of \(\mathbf{B}\) — matrix multiplication is just matrix–vector multiplication done once per column.
The shapes must line up in the middle — the inner dimensions have to agree, and the outer ones survive into the result:
\[ \underset{m \times n}{\mathbf{A}} \;\; \underset{n \times p}{\mathbf{B}} \;\longrightarrow\; \underset{m \times p}{\mathbf{C}} \]
Worked example. Entry \((1,1)\) of \(\mathbf{A}\mathbf{B}\) is row 1 of \(\mathbf{A}\) dotted with column 1 of \(\mathbf{B}\): \((2,1)\cdot(1,2) = 2 + 2 = 4\). Doing all four:
\[ \mathbf{A}\mathbf{B} = \begin{bmatrix} 2 & 1 \\ 0 & 3\end{bmatrix} \begin{bmatrix} 1 & 4 \\ 2 & 1\end{bmatrix} = \begin{bmatrix} 4 & 9 \\ 6 & 3 \end{bmatrix} \]
A %*% B [,1] [,2]
[1,] 4 9
[2,] 6 3
WarningWatch out
* is elementwise multiplication in R, not matrix multiplication. It is a valid operation with a different meaning (the Hadamard product), and it fails silently when the shapes happen to match:
A * B # elementwise — almost never what you want [,1] [,2]
[1,] 2 4
[2,] 0 3A %*% B # matrix product [,1] [,2]
[1,] 4 9
[2,] 6 3Non-square shapes make the rule concrete:
[1] 2 2[1] 3 3C %*% D [,1] [,2]
[1,] 2 7
[2,] 10 10D %*% C [,1] [,2] [,3]
[1,] 5 5 4
[2,] 9 3 12
[3,] 4 3 4Both products exist and they are not even the same size.
5.9 Properties and non-properties
Matrix multiplication keeps most of the algebra you expect:
\[ \begin{aligned} (\mathbf{A}\mathbf{B})\mathbf{C} &= \mathbf{A}(\mathbf{B}\mathbf{C}) &&\text{associative} \\ \mathbf{A}(\mathbf{B} + \mathbf{C}) &= \mathbf{A}\mathbf{B} + \mathbf{A}\mathbf{C} &&\text{distributive} \\ \mathbf{A}\mathbf{I} = \mathbf{I}\mathbf{A} &= \mathbf{A} &&\text{identity} \end{aligned} \]
But it loses the one you rely on most without noticing.
Matrix multiplication is not commutative. In general \(\mathbf{A}\mathbf{B} \neq \mathbf{B}\mathbf{A}\):
Different matrices entirely. This is not a quirk of these particular numbers — it is the normal case, and it makes sense the moment you think of matrices as functions: rotating then stretching is not the same as stretching then rotating.
Two more traps:
- \(\mathbf{A}\mathbf{B} = \mathbf{0}\) does not imply \(\mathbf{A} = \mathbf{0}\) or \(\mathbf{B} = \mathbf{0}\).
- \(\mathbf{A}\mathbf{B} = \mathbf{A}\mathbf{C}\) does not imply \(\mathbf{B} = \mathbf{C}\). You cannot cancel a matrix unless you know it is invertible.
And the transpose rule promised earlier:
5.10 A matrix as a linear map
Here is the shift that makes the rest of linear algebra work.
Stop reading \(\mathbf{A}\) as a grid of numbers. Read it as a function: feed it a vector \(\mathbf{x}\), get back a vector \(\mathbf{A}\mathbf{x}\). The function is linear, meaning it respects addition and scaling:
\[ \mathbf{A}(\mathbf{x} + \mathbf{y}) = \mathbf{A}\mathbf{x} + \mathbf{A}\mathbf{y} \qquad \mathbf{A}(c\mathbf{x}) = c(\mathbf{A}\mathbf{x}) \]
Those two conditions have a strong geometric consequence: straight lines stay straight, parallel lines stay parallel, and the origin stays put. A linear map can rotate, stretch, shear, reflect, or flatten space — but it cannot bend or translate it.
The clearest way to see what a particular matrix does is to watch what it does to the unit square.
Where did the corners go? The corner at \((1,0)\) is \(\mathbf{e}_1\), the first standard basis vector, and \(\mathbf{A}\mathbf{e}_1\) picks out the first column of \(\mathbf{A}\). Likewise \(\mathbf{A}\mathbf{e}_2\) is the second column. So:
The columns of \(\mathbf{A}\) are the images of the basis vectors. To know what a matrix does, look at where it sends \(\mathbf{e}_1\) and \(\mathbf{e}_2\) — that is literally what the columns are.
Once matrices are functions, matrix multiplication stops being an arbitrary rule. \(\mathbf{A}\mathbf{B}\) is function composition: apply \(\mathbf{B}\) first, then \(\mathbf{A}\). That is why the order matters, and why \((\mathbf{A}\mathbf{B})^\top\) reverses.
Two special cases are worth seeing.
That second one is the geometric meaning of singular: the map is not reversible, because different inputs land on the same output. Its columns, \((1,2)\) and \((2,4)\), lie on the same line — the second is twice the first — so every linear combination of them lies on that line too.
NoteIn machine learning
Every layer of a neural network is a linear map followed by a non-linear function. The linear map is a matrix multiply; the non-linearity (ReLU, sigmoid) is what lets the network bend space rather than only rotating and stretching it. Stack linear maps without non-linearities and you get one linear map — a network of any depth collapses into a single matrix.
5.11 Block matrices
A matrix can be partitioned into sub-matrices, or blocks, and multiplied blockwise as though the blocks were single numbers — provided the partitions are compatible.
\[ \begin{bmatrix} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{21} & \mathbf{A}_{22} \end{bmatrix} \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \end{bmatrix} = \begin{bmatrix} \mathbf{A}_{11}\mathbf{x}_1 + \mathbf{A}_{12}\mathbf{x}_2 \\ \mathbf{A}_{21}\mathbf{x}_1 + \mathbf{A}_{22}\mathbf{x}_2 \end{bmatrix} \]
This is mostly a notational convenience, but it is how you keep track of intercepts, grouped features, and multi-output models without drowning in subscripts. In R, rbind() and cbind() build blocks:
cbind(intercept = 1, patients) # prepend a column of ones intercept age sbp bmi
p1 1 64 138 27.4
p2 1 51 122 31.9
p3 1 73 155 24.1
p4 1 45 118 29.65.12 The inverse
The inverse of a square matrix \(\mathbf{A}\), written \(\mathbf{A}^{-1}\), is the matrix that undoes it:
\[ \mathbf{A}\mathbf{A}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I} \]
Not every matrix has one. A matrix that does is invertible or non-singular; one that does not is singular, and geometrically that means it collapses space (Figure 5.6) so there is nothing to undo.
For a \(2 \times 2\) matrix there is a closed form worth memorizing:
\[ \begin{bmatrix} a & b \\ c & d \end{bmatrix}^{-1} = \frac{1}{ad - bc}\begin{bmatrix} d & -b \\ -c & a \end{bmatrix} \tag{5.4}\]
Worked example. For \(\mathbf{A} = \begin{bmatrix}2 & 1\\0 & 3\end{bmatrix}\), \(ad - bc = 6 - 0 = 6\), so
\[ \mathbf{A}^{-1} = \frac{1}{6}\begin{bmatrix}3 & -1\\0 & 2\end{bmatrix} = \begin{bmatrix}0.5 & -0.1\overline{6}\\0 & 0.\overline{3}\end{bmatrix} \]
Note the quantity \(ad - bc\) in the denominator: when it is zero there is no inverse. That quantity is the determinant, and Section 5.13 explains why it decides invertibility.
WarningWatch out
Do not compute an inverse to solve a linear system. To solve \(\mathbf{A}\mathbf{x} = \mathbf{b}\), write solve(A, b) — not solve(A) %*% b. The two-argument form is faster and numerically far more stable, because it factors the matrix rather than inverting it.
Same answer here. On an ill-conditioned matrix it will not be, and we return to why in Section 6.10.
5.13 Trace and determinant
Two numbers summarize a square matrix.
5.13.1 Trace
The trace is the sum of the diagonal entries:
\[ \operatorname{tr}(\mathbf{A}) = \sum_{i=1}^{n} a_{ii} \]
Its useful property is that it is invariant to the order of a product, even though the product itself is not:
\[ \operatorname{tr}(\mathbf{A}\mathbf{B}) = \operatorname{tr}(\mathbf{B}\mathbf{A}) \]
5.13.2 Determinant
The determinant \(\det(\mathbf{A})\), sometimes written \(|\mathbf{A}|\), is a single number measuring how much the map scales area (in 2D) or volume (in higher dimensions). For \(2 \times 2\):
\[ \det\begin{bmatrix} a & b \\ c & d \end{bmatrix} = ad - bc \]
Worked example. \(\det(\mathbf{A}) = (2)(3) - (1)(0) = 6\).
det(A)[1] 6The unit square has area 1, and its image under \(\mathbf{A}\) has area 6 — the determinant, exactly.
Three consequences follow directly from reading it as a scaling factor:
| \(\det(\mathbf{A})\) | Meaning |
|---|---|
| \(0\) | area collapses; the matrix is singular and has no inverse |
| negative | orientation flips (a reflection) |
| \(\|\det\| > 1\) | the map expands; \(< 1\), it shrinks |
The determinant of the singular matrix from Figure 5.6 is zero, as its flattened image requires:
det(Sing)[1] 0And because applying two maps in sequence scales area twice,
\[ \det(\mathbf{A}\mathbf{B}) = \det(\mathbf{A})\det(\mathbf{B}) \]
TipGoing deeper: why determinants are rare in practice
Determinants are wonderful for reasoning and terrible for computing. The textbook cofactor expansion costs \(O(n!)\) operations — hopeless beyond about \(n = 10\) — and even the practical \(O(n^3)\) route via factorization overflows or underflows the moment \(n\) gets large, because you are multiplying \(n\) numbers together.
Software that needs to know whether a matrix is invertible checks its condition number or its rank instead, never det(A) == 0. Use the determinant to understand; use a factorization to compute.
5.14 Computational cost
Multiplying an \(m \times n\) by an \(n \times p\) matrix requires \(mnp\) multiply-add operations. For two \(n \times n\) matrices that is \(O(n^3)\): doubling the size takes eight times as many operations.
sizes <- c(200, 400, 800)
times <- sapply(sizes, function(n) {
M <- matrix(rnorm(n * n), nrow = n)
unname(system.time(M %*% M)["elapsed"])
})
data.frame(
n = sizes,
seconds = times,
ratio = c(NA, round(times[-1] / times[-3], 1))
) n seconds ratio
1 200 0.005 NA
2 400 0.033 6.6
3 800 0.250 7.6draw_line(
x = sizes,
y = times,
xlab = "n",
ylab = "seconds"
)
WarningWatch out
Wall-clock time and operation count are related, but they are not the same thing, and the gap between them is easy to trip over.
The ratio column above should land somewhere near 8 — and it does partly because these sizes were chosen to make it. Start at \(n = 100\) instead and the first ratio collapses, because the whole multiply takes about as long as the timer can resolve and allocating the result dominates the arithmetic. At the other end, R hands large multiplies to an optimized BLAS that splits the work across CPU cores and reuses data already sitting in cache, so ratios can fall well below 8 — until the matrices stop fitting in cache, at which point they jump back above it.
\(O(n^3)\) is a statement about the number of operations. It correctly tells you that an \(n = 10{,}000\) multiply is out of reach while \(n = 1{,}000\) is routine. It is not a prediction of seconds on a particular machine, and reading it as one will mislead you in both directions.
Two practical consequences:
- Order your products. \((\mathbf{A}\mathbf{B})\mathbf{C}\) and \(\mathbf{A}(\mathbf{B}\mathbf{C})\) give the same answer at wildly different costs when the shapes are uneven. Exercise 6 makes this concrete.
- Never form a matrix you only need to apply. If you need \(\mathbf{A}^{-1}\mathbf{b}\), solve; if you need the diagonal of a product, compute only the diagonal.
5.15 Summary
| Operation | Notation | R | Shape |
|---|---|---|---|
| Transpose | \(\mathbf{A}^\top\) | t(A) |
\(n \times m\) |
| Addition | \(\mathbf{A} + \mathbf{B}\) | A + B |
same |
| Scaling | \(c\mathbf{A}\) | c * A |
same |
| Matrix–vector | \(\mathbf{A}\mathbf{x}\) | A %*% x |
\(m \times 1\) |
| Matrix–matrix | \(\mathbf{A}\mathbf{B}\) | A %*% B |
\(m \times p\) |
| Elementwise | \(\mathbf{A} \odot \mathbf{B}\) | A * B |
same |
| Inverse | \(\mathbf{A}^{-1}\) | solve(A) |
\(n \times n\) |
| Solve | \(\mathbf{A}^{-1}\mathbf{b}\) | solve(A, b) |
\(n \times 1\) |
| Trace | \(\operatorname{tr}(\mathbf{A})\) | sum(diag(A)) |
scalar |
| Determinant | \(\det(\mathbf{A})\) | det(A) |
scalar |
| \(\mathbf{A}^\top\mathbf{A}\) | \(\mathbf{A}^\top\mathbf{A}\) | crossprod(A) |
\(n \times n\) |
5.16 Exercises
1. Let \(\mathbf{M} = \begin{bmatrix}1 & 2\\3 & 4\end{bmatrix}\) and \(\mathbf{N} = \begin{bmatrix}0 & 1\\1 & 0\end{bmatrix}\). Compute \(\mathbf{M}\mathbf{N}\) and \(\mathbf{N}\mathbf{M}\) by hand. What does \(\mathbf{N}\) do?
TipSolution
\(\mathbf{M}\mathbf{N} = \begin{bmatrix}2 & 1\\4 & 3\end{bmatrix}\) — the columns of \(\mathbf{M}\) are swapped. \(\mathbf{N}\mathbf{M} = \begin{bmatrix}3 & 4\\1 & 2\end{bmatrix}\) — the rows are swapped.
\(\mathbf{N}\) is a permutation matrix. Multiplying on the right permutes columns; on the left, rows. A clean demonstration that order matters.
2. For \(\mathbf{A} = \begin{bmatrix}2 & 1\\0 & 3\end{bmatrix}\), find \(\mathbf{A}^{-1}\) by hand using Equation 5.4 and verify that \(\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}\).
TipSolution
3. Which of these are invertible? Answer from the determinant, then say what each does geometrically.
\[ \mathbf{P} = \begin{bmatrix}3 & 0\\0 & 3\end{bmatrix} \quad \mathbf{R} = \begin{bmatrix}1 & 2\\2 & 4\end{bmatrix} \quad \mathbf{T} = \begin{bmatrix}0 & 1\\-1 & 0\end{bmatrix} \]
TipSolution
[1] 9 0 1\(\det(\mathbf{P}) = 9\): invertible, scales everything by 3 in both directions, so area grows ninefold. \(\det(\mathbf{R}) = 0\): singular — its second column is twice the first, so it flattens the plane onto a line. \(\det(\mathbf{T}) = 1\): invertible, a rotation, which preserves area exactly.
4. Show that \(\mathbf{A}^\top\mathbf{A}\) is symmetric for any matrix \(\mathbf{A}\), square or not.
TipSolution
A matrix is symmetric when it equals its own transpose. Apply the reversal rule twice:
\[ (\mathbf{A}^\top\mathbf{A})^\top = \mathbf{A}^\top(\mathbf{A}^\top)^\top = \mathbf{A}^\top\mathbf{A} \]
G <- crossprod(C) # C is 2x3, so this is 3x3
G [,1] [,2] [,3]
[1,] 10 5 12
[2,] 5 5 4
[3,] 12 4 16[1] TRUEThis matters more than it looks: \(\mathbf{X}^\top\mathbf{X}\) is at the center of the normal equations, and its symmetry is what makes least squares tractable.
5. The design matrix patients has one row per patient. Compute the mean of each column two ways: with colMeans(), and as a matrix–vector product.
TipSolution
Multiplying on the left by a row vector of \(1/n\) averages the columns.
Summary statistics are linear maps. That is why they compose so cleanly with the rest of the machinery — centering a data matrix, for example, is a single matrix subtraction.
6. You need \(\mathbf{A}\mathbf{B}\mathbf{C}\) where \(\mathbf{A}\) is \(1000 \times 2\), \(\mathbf{B}\) is \(2 \times 1000\), and \(\mathbf{C}\) is \(1000 \times 2\). Count the multiply-adds for \((\mathbf{A}\mathbf{B})\mathbf{C}\) versus \(\mathbf{A}(\mathbf{B}\mathbf{C})\).
TipSolution
Multiplying \(m \times n\) by \(n \times p\) costs \(mnp\).
- \((\mathbf{A}\mathbf{B})\mathbf{C}\): \(\mathbf{A}\mathbf{B}\) costs \(1000 \times 2 \times 1000 = 2{,}000{,}000\) and produces a \(1000 \times 1000\) matrix; multiplying that by \(\mathbf{C}\) costs \(1000 \times 1000 \times 2 = 2{,}000{,}000\). Total 4,000,000.
- \(\mathbf{A}(\mathbf{B}\mathbf{C})\): \(\mathbf{B}\mathbf{C}\) costs \(2 \times 1000 \times 2 = 4{,}000\) and produces a \(2 \times 2\) matrix; multiplying \(\mathbf{A}\) by that costs \(1000 \times 2 \times 2 = 4{,}000\). Total 8,000.
set.seed(7)
Am <- matrix(rnorm(1000 * 2), nrow = 1000)
Bm <- matrix(rnorm(2 * 1000), nrow = 2)
Cm <- matrix(rnorm(1000 * 2), nrow = 1000)
# 100 repetitions: one pass of the good grouping is
# too fast to measure reliably
c(
left = unname(system.time(
for (i in 1:100) (Am %*% Bm) %*% Cm
)["elapsed"]),
right = unname(system.time(
for (i in 1:100) Am %*% (Bm %*% Cm)
)["elapsed"])
) left right
0.271 0.003 Same answer, 500 times fewer operations, and the bad grouping also allocates a \(1000 \times 1000\) intermediate. Associativity is free to exploit — but you have to choose to.