R (4.6) Module 3 – Factors, matrices, lists and data frames
Matrices, arrays and linear algebra
Build R matrices and arrays, index them without losing dimensions, summarise rows and columns, multiply with %*% and solve least squares stably with QR.
What you will learn
- Build matrices and arrays, and index them by row, column and slice without dropping dimensions by accident
- Apply matrix algebra with t(), %*%, crossprod() and solve(), and name the main decompositions
- Summarise rows and columns with rowSums(), colMeans() and apply()
- Compute least-squares coefficients stably with QR instead of an explicit inverse
Before you start
What you’ll learn · Before you start
What you will learn
- Build matrices and arrays, and index them by row, column and slice without dropping dimensions by accident
- Apply matrix algebra with t(), %*%, crossprod() and solve(), and name the main decompositions
- Summarise rows and columns with rowSums(), colMeans() and apply()
- Compute least-squares coefficients stably with QR instead of an explicit inverse
Before you start
On this page
A function that turns the rows of a matrix into proportions works on 100 rows and breaks on 1, because m[1, ]
quietly returns a plain vector instead of a matrix. And a regression solved through the normal equations, on a data
set that ships with R, keeps only about 8 of the 15 to 16 significant digits a double can hold. Matrices are where
R’s vector rules meet linear algebra, and both problems come from rules worth knowing.
On runs of R 4.6.1 you will build matrices and arrays and index them without losing dimensions, summarise their rows and columns, multiply and solve with the right operators, and compute least-squares coefficients with QR, measured against NIST’s reference values.
A matrix is a vector with dimensions, filled column by column
matrix(1:6, nrow = 3) makes a matrix: A vector with two dimensions, set by its dim attribute, so that its values sit in rows and columns. Like every atomic vector it holds values of one type. of 3 rows, and from the six values R works out that it
needs 2 columns. The values go down the first column, then down the second, unless you ask for byrow = TRUE. Step
through both fills:
The drawing shows the whole scene; its 5 steps are listed below.
Every step as text
- Six cells. nrow = 3 and six values give 3 rows and 6 / 3 = 2 columns. R works out the number of columns itself.
- Column 1 first. The values go down the first column, 1, 2 and 3, because byrow is FALSE by default.
- Then column 2. The next three values, 4, 5 and 6, fill the second column, so the first row reads 1 4.
- byrow = TRUE. With byrow = TRUE the same values go across the rows instead, 1 2, then 3 4, then 5 6.
- Stored by column. Without byrow, R keeps the matrix as the vector 1:6 with a dim attribute of 3 and 2, so m[5], the fifth stored value, is the cell in row 2 of column 2.
The same matrix, built and taken apart at the console:
from session import run, version
version()
# R reads matrix_basics.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f matrix_basics.R
""")Output
# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f matrix_basics.R
> m <- matrix(1:6, nrow = 3)
> m
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
> matrix(1:6, nrow = 3, byrow = TRUE)
[,1] [,2]
[1,] 1 2
[2,] 3 4
[3,] 5 6
> dim(m)
[1] 3 2
> attributes(m)
$dim
[1] 3 2
> # A matrix is a vector with a dim attribute.
> v <- 1:6
> dim(v) <- c(3, 2)
> identical(v, m)
[1] TRUE
> length(m)
[1] 6
> # Five values do not fill two rows evenly.
> matrix(1:5, nrow = 2)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 1
Warning message:
In matrix(1:5, nrow = 2) :
data length [5] is not a sub-multiple or multiple of the number of rows [2]
> # Names for the rows and the columns.
> rownames(m) <- c("a", "b", "c")
> colnames(m) <- c("x", "y")
> m
x y
a 1 4
b 2 5
c 3 6
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 matrix_basics_session.py
attributes(m) holds a single entry, dim, with the values 3 and 2. A matrix is an ordinary vector with that
attribute: Extra information R stores beside a vector's values, such as its names, its dimensions or its class. attributes(x) lists them. More in “Atomic vectors, types and coercion”: giving the vector 1:6 the dimensions c(3, 2) produced a matrix identical to
matrix(1:6, nrow = 3), and length(m) is still 6. The values are stored as one
atomic vector: A vector whose elements all have one of R's six basic types: logical, integer, double, complex, character or raw. A single value such as 42 is an atomic vector of length 1. More in “Atomic vectors, types and coercion”, down the columns, in
column-major order: Storing a matrix column after column: all of column 1, then all of column 2. R stores matrices and arrays this way, so m[5] is the fifth value counted down the columns.. So a matrix holds one type, and a single character value
would turn every number into text.
The help page of matrix says that byrow = FALSE, the default, fills by columns. matrix(1:5, nrow = 2) shows what
happens when the values do not fit: R still builds a 2 × 3 matrix, starts the five values again to fill the sixth
cell, and warns that the data length is not a sub-multiple or multiple of the number of rows. Treat that warning as an
error in your own code; a matrix with a recycled value in it is rarely what you meant. rownames() and colnames()
name the rows and columns, and the names then work as indices.
Checkpoint
Check: Column-major storage
Indexing by row and column, and keeping both dimensions
m[i, j] takes rows i and columns j, and each index follows the rules of vector indexing: positions, negative
positions, logical vectors or names. Leave one index empty and you get every row, or every column. Step through what
each form selects on a matrix with named rows and columns:
The drawing shows the whole scene; its 6 steps are listed below.
Every step as text
- m. The matrix has rows a, b and c and columns x and y. Each box is a cell, named by its row and its column, with its value below.
- m[2, 1]. A row index and a column index pick one cell, row 2 and column 1, the value 2. Names work too, so m["b", "x"] is the same cell.
- m[2, ]. An empty column index means every column, so this is all of row 2. R returns a plain vector, named x and y, not a matrix.
- m[, "y"]. An empty row index means every row, so this is the whole column y, again as a plain vector, named a, b and c.
- m[m > 3]. A logical matrix as the index picks the cells where it is TRUE and returns their values as a vector, in column order.
- drop = FALSE. m[2, , drop = FALSE] picks the same row as before but keeps both dimensions, so the result is a matrix of 1 row and 2 columns.
The run of those calls, and of a filter that happens to keep one row:
from session import run, version
version()
# R reads indexing.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f indexing.R
""")Output
# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f indexing.R
> m <- matrix(1:6, nrow = 3, dimnames = list(c("a", "b", "c"), c("x", "y")))
> m[2, 1]
[1] 2
> m["b", "y"]
[1] 5
> m[2, ]
x y
2 5
> m[, "y"]
a b c
4 5 6
> m[m > 3]
[1] 4 5 6
> # One row or one column comes back as a plain vector.
> dim(m[2, ])
NULL
> m[2, , drop = FALSE]
x y
b 2 5
> dim(m[2, , drop = FALSE])
[1] 1 2
> # A filter that happens to keep one row.
> big <- m[m[, "x"] > 2, ]
> big
x y
3 6
> ncol(big)
NULL
> big <- m[m[, "x"] > 2, , drop = FALSE]
> ncol(big)
[1] 2
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 indexing_session.py
m[2, ] comes back as a plain vector with names, and dim(m[2, ]) is NULL. That is the drop argument at work:
the help page of Extract says that with drop = TRUE, the default for matrices, the result is coerced to the lowest
possible dimension, so one row or one column loses its matrix shape. Indexing with a logical matrix, m[m > 3],
returns the selected values as a vector in column order.
The danger is code that expects a matrix and usually gets one. m[m[, "x"] > 2, ] keeps the rows whose x is above
2; on this matrix only row c qualifies, so the result is a vector, and ncol() of it is NULL. A function written
this way passes every test with several matching rows and fails on the first day with one. Adding drop = FALSE
keeps a 1 × 2 matrix, and ncol() is 2 again.
Taking rows of a matrix inside a function without drop = FALSE
x[keep, ] returns a vector when keep selects one row, and the next ncol(x) or x %*% w misbehaves. Write
x[keep, , drop = FALSE] in any function that takes a matrix, and test it on a matrix of one row.
Checkpoint
Check: Indexing and drop
rowSums() and colMeans() first, apply() for the rest
apply(m, 1, f) calls f once for each row of m, and apply(m, 2, f) once for each column; the second argument
is the margin, the dimension the function runs along. Step through both margins with sum():
The drawing shows the whole scene; its 6 steps are listed below.
Every step as text
- The matrix. matrix(1:6, nrow = 3) has three rows and two columns. apply() calls a function once for each row, with MARGIN = 1, or once for each column, with MARGIN = 2.
- Row 1. With MARGIN = 1, sum() first receives row 1, the values 1 and 4.
- Row 2. Then row 2, the values 2 and 5.
- Row 3. Then row 3, the values 3 and 6. apply(m, 1, sum) returns the three sums in order, 5 7 9, which is what rowSums(m) returns.
- Column 1. With MARGIN = 2, sum() receives each column instead. Column 1 holds 1, 2 and 3.
- Column 2. Column 2 holds 4, 5 and 6, so apply(m, 2, sum) returns 6 15, the same as colSums(m).
The run adds the shortcuts and a function that returns more than one value:
from session import run, version
version()
# R reads margins.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f margins.R
""")Output
# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f margins.R
> m <- matrix(1:6, nrow = 3)
> apply(m, 1, sum)
[1] 5 7 9
> apply(m, 2, sum)
[1] 6 15
> rowSums(m)
[1] 5 7 9
> colMeans(m)
[1] 2 5
> # A function that returns a vector gives one column per row.
> apply(m, 1, function(r) r / sum(r))
[,1] [,2] [,3]
[1,] 0.2 0.2857143 0.3333333
[2,] 0.8 0.7142857 0.6666667
> # Divide each row by its sum and keep the shape.
> m / rowSums(m)
[,1] [,2]
[1,] 0.2000000 0.8000000
[2,] 0.2857143 0.7142857
[3,] 0.3333333 0.6666667
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 margins_session.py
For sums and means, rowSums(), colSums(), rowMeans() and colMeans() give the same answers as apply(); their
help page says they are equivalent to apply() with sum or mean but a lot faster, so use them first and keep
apply() for other functions, such as max or sd.
When the function returns a vector, apply() puts each result in a column. The proportions of each row came back as
2 rows and 3 columns, the transpose of what was meant, and code that reads the result by rows reads the wrong numbers.
Wrap it in t(), or avoid apply() here: m / rowSums(m) divides each row by its own sum, because the three sums
are recycling: R's rule for vectors of different lengths in an element-by-element operation: the shorter one is repeated from its start until it is as long as the longer one. R warns only when the lengths do not divide evenly. More in “Vector arithmetic, recycling and sequences” down each column of the column-major matrix, so the first sum always meets
row 1. The result keeps the shape of m, even when m has a single row.
Key idea
The margin is the dimension apply() runs along: 1 for rows, 2 for columns. Use rowSums() and colMeans() for
sums and means, and expect apply() to return one column per row when the function returns a vector.
Checkpoint
Check: Row and column summaries
An array adds a third index: a plate read three times
An array: A vector with any number of dimensions, set by its dim attribute. A matrix is the two-dimensional case, and arr[i, j, k] picks one value of a three-dimensional array. takes the same idea to three or more dimensions. A plate of 96 wells, 8 rows by 12 columns, read at 0, 2 and 4 hours, fits an array of dimensions 8, 12 and 3, with a name for each row, column and time. Step through what each slice of such an array returns:
The drawing shows the whole scene; its 5 steps are listed below.
Every step as text
- dim(plate). The array has three dimensions, 8 rows, 12 columns and 3 reading times, so every value has three indices, such as plate["A", "1", "2h"].
- plate[, , "2h"]. Fixing the third index to 2h and leaving the others empty takes the whole plate read at 2 hours. One dimension is fixed, so R drops it and returns a matrix.
- plate["A", , ]. Fixing the row to A takes row A of every plate, 12 wells at 3 times. The row dimension is dropped, so the result is a 12 x 3 matrix.
- drop = FALSE. plate["A", , , drop = FALSE] selects the same 36 values but keeps all three dimensions, so code written for arrays keeps working.
- colMeans(plate). colMeans() averages over the first dimension, the 8 rows, and keeps the others, so it returns the mean of each column at each time, a 12 x 3 matrix.
The run fills a plate with simulated readings, seeded with set.seed(1) so that every run gets the same numbers, then
slices it and averages it:
from session import run, version
version()
# R reads plate.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f plate.R
""")Output
# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f plate.R
> # Simulated readings of a plate of 96 wells (rows A to H,
> # columns 1 to 12), read at 0, 2 and 4 hours.
> set.seed(1)
> plate <- array(round(runif(8 * 12 * 3, 0.1, 1.2), 2),
+ dim = c(8, 12, 3),
+ dimnames = list(LETTERS[1:8], 1:12, c("0h", "2h", "4h")))
> dim(plate)
[1] 8 12 3
> plate["A", "1", "2h"]
[1] 0.6
> dim(plate[, , "2h"])
[1] 8 12
> dim(plate["A", , ])
[1] 12 3
> dim(plate["A", , , drop = FALSE])
[1] 1 12 3
> # Column means at each time: one value per column and time.
> col_means <- colMeans(plate)
> dim(col_means)
[1] 12 3
> round(col_means[1:3, ], 3)
0h 2h 4h
1 0.764 0.726 0.643
2 0.570 0.744 0.466
3 0.760 0.490 0.770
> all.equal(col_means, apply(plate, c(2, 3), mean))
[1] TRUE
> round(apply(plate, 3, mean), 3)
0h 2h 4h
0.667 0.673 0.593
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 plate_session.py
Every value now has three indices, as in plate["A", "1", "2h"], and the drop rule is the same as for matrices:
fix one index and that dimension disappears. plate[, , "2h"] is the 8 × 12 matrix of one reading time, and
plate["A", , ] is a 12 × 3 matrix, row A across all columns and times. drop = FALSE keeps the array shape,
1 × 12 × 3, which is what a function written for arrays needs.
Summaries work along any set of dimensions. colMeans(plate) averages over the first dimension, the 8 rows of the
plate, and keeps the other two, so it returns a 12 × 3 matrix, the mean of every column at every time.
apply(plate, c(2, 3), mean) says the same thing in words, “keep dimensions 2 and 3”, and all.equal() confirms that
the two results agree. apply(plate, 3, mean) keeps only the time and gives one mean per reading: 0.667, 0.673 and
0.593.
Checkpoint
Check: Arrays of three dimensions
A matrix product combines rows with columns
* on two matrices multiplies the values in the same places, as it does for vectors. The
matrix product: The product of linear algebra, A %*% B in R: entry [i, j] is the sum of row i of A times column j of B, so the number of columns of A must equal the number of rows of B. is %*%: each entry combines a row of the first matrix with a
column of the second. Step through one entry of A %*% B, then the same place in A * B:
The drawing shows the whole scene; its 5 steps are listed below.
Every step as text
- Row times column. Entry [1, 2] of A %*% B combines row 1 of A with column 2 of B. The two have the same length, 2, which is why the columns of A must match the rows of B.
- First pair. The first value of the row times the first value of the column, 2 × 4 = 8.
- Second pair. The second value of the row times the second value of the column, 1 × 2 = 2.
- The sum. The products add up to 8 + 2 = 10, the entry in row 1, column 2 of A %*% B, as the run prints.
- A * B instead. A * B multiplies only the values in the same place, A[1, 2] × B[1, 2] = 1 × 4 = 4, and nothing is added. That is element-wise arithmetic, as for vectors.
The run multiplies both ways, then solves a system of equations:
from session import run, version
version()
# R reads algebra.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f algebra.R
""")Output
# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f algebra.R
> A <- matrix(c(2, 1, 1, 3), nrow = 2)
> B <- matrix(c(1, 0, 4, 2), nrow = 2)
> A * B
[,1] [,2]
[1,] 2 4
[2,] 0 6
> A %*% B
[,1] [,2]
[1,] 2 10
[2,] 1 10
> t(A)
[,1] [,2]
[1,] 2 1
[2,] 1 3
> crossprod(A)
[,1] [,2]
[1,] 5 5
[2,] 5 10
> identical(crossprod(A), t(A) %*% A)
[1] TRUE
> # Solve A x = b without forming the inverse.
> b <- c(5, 10)
> solve(A, b)
[1] 1 3
> # The inverse is exact on paper, not in a double.
> A %*% solve(A)
[,1] [,2]
[1,] 1.000000e+00 0
[2,] -5.551115e-17 1
> isTRUE(all.equal(A %*% solve(A), diag(2)))
[1] TRUE
> basename(extSoftVersion()[["BLAS"]])
[1] "libopenblas.0.dylib"
> # A badly conditioned 10 x 10 system with the answer x = 1.
> H <- 1 / outer(1:10, 1:10, "+")
> b <- H %*% rep(1, 10)
> signif(max(abs(solve(H, b) - 1)), 2)
[1] 0.0024
> signif(max(abs(solve(H) %*% b - 1)), 2)
[1] 0.027
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 algebra_session.py
t(A) transposes. crossprod(A) is t(A) %*% A, which its help page calls formally equivalent but faster, and the
run shows the two are identical here. To solve A %*% x = b, call solve(A, b). An Introduction to R points out
that computing solve(A) %*% b instead, the inverse times b, is both inefficient and potentially unstable, and the
run measures it on a badly conditioned 10 × 10 system whose answer is ten 1s: solve(H, b) is off by at most 0.0024,
the inverse route by 0.027, eleven times more.
Even a small, well-behaved inverse is not exact in doubles. A %*% solve(A) should be the identity matrix, and it
printed -5.551115e-17 where a 0 belongs. Such last digits can differ between builds of R, which may use different
libraries of linear-algebra routines, the BLAS: Basic Linear Algebra Subprograms: the library of routines for matrix products and related work that R calls. Builds of R can use different BLAS libraries, which can change the last digits of results.; this run used OpenBLAS, and R’s installation
manual explains how R is built with an external BLAS and warns that some parallel ones are not even deterministic. So
compare matrices with isTRUE(all.equal(...)), never with == or identical().
Reading crossprod() as the cross product of vectors
crossprod(x, y) is t(x) %*% y, a matrix of dot products, not the cross product of two 3-D vectors from physics.
For that, write the formula yourself, or use a package that provides it.
R has the main matrix decompositions built in, each solving one kind of problem:
qr(X)factorsXinto an orthogonal and a triangular matrix, the basis of least squares (next section).chol(S)factors a symmetric positive-definite matrix into a triangular one and its transpose, a fast way to solve systems with such matrices.eigen(S)gives the eigenvalues and eigenvectors, whichprincomp()uses for principal components.svd(X)gives the singular value decomposition, from whichkappa(X, exact = TRUE)computes the condition number used in the next section.
Checkpoint
Check: Matrix products and solving
Least squares: solve with QR, not with the normal equations
Fitting a straight line, or any linear model, means finding the coefficients b that make X %*% b as close as
possible to y, where X is the design matrix, one column per term. Algebra gives a formula, the
normal equations: The system t(X) %*% X %*% b = t(X) %*% y, whose solution b is the least-squares fit of y on X. Solving it directly works with the square of the condition number of X., t(X) %*% X %*% b = t(X) %*% y, and it is tempting to solve
them directly. The other route factors X itself with a QR decomposition: Writing a matrix X as Q %*% R, with Q orthogonal and R triangular. Least squares solved this way, as lm() and qr.coef() do, works on X itself instead of t(X) %*% X.. Follow
both routes on a hard case:
The drawing shows the whole scene; its 5 steps are listed below.
Every step as text
- The design matrix. X holds a column of 1s and six economic series that move together. Its condition number, kappa(X, exact = TRUE), is about 23.8 million: errors in the last digits can grow that much.
- Forming t(X) %*% X. crossprod(X) squares the condition number, to about 5.7e14, so the system the normal equations solve is far more sensitive than the data themselves.
- Solving it. solve() still returns an answer, but it matches NIST's reference coefficients to only 7.6 to 8.5 significant digits, about half of what a double holds.
- The QR route. qr() factors X itself into an orthogonal Q and a triangular R, so the condition number is never squared. qr.coef() matches the reference values to 13.4 to 15.7 digits.
- What lm() does. lm() fits with QR too, its only fitting method, so it gets the same digits as qr.coef() without any extra work.
The run fits a five-point dose-response line three ways, then Longley’s data, which ships with R:
from session import run, version
version()
# R reads least_squares.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f least_squares.R
""")Output
# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f least_squares.R
> # Five doses and one response each.
> dose <- c(0, 1, 2, 4, 8)
> resp <- c(2.1, 3.9, 6.2, 9.8, 18.1)
> X <- cbind(intercept = 1, dose)
> b_normal <- solve(crossprod(X), crossprod(X, resp))
> b_qr <- qr.coef(qr(X), resp)
> b_lm <- coef(lm(resp ~ dose))
> cbind(normal = drop(b_normal), qr = b_qr, lm = b_lm)
normal qr lm
intercept 2.02 2.02 2.02
dose 2.00 2.00 2.00
> # Longley's data, which ships with R, is highly collinear.
> L <- cbind(1, as.matrix(longley[, 1:6]))
> y <- longley$Employed
> kappa(L, exact = TRUE)
[1] 23845862
> kappa(crossprod(L), exact = TRUE)
[1] 5.685668e+14
> # NIST's reference coefficients, in the units of R's copy.
> nist <- c(-3482.25863459582, 0.0150618722713733,
+ -0.0358191792925910, -0.0202022980381683,
+ -0.0103322686717359, -0.0511041056535807, 1.82915146461355)
> correct_digits <- function(b) {
+ round(-log10(abs(b - nist) / abs(nist)), 1)
+ }
> data.frame(
+ term = c("(Intercept)", names(longley)[1:6]),
+ normal = correct_digits(drop(solve(crossprod(L), crossprod(L, y)))),
+ qr = correct_digits(qr.coef(qr(L), y)),
+ lm = correct_digits(unname(coef(lm(Employed ~ ., data = longley)))),
+ row.names = NULL
+ )
term normal qr lm
1 (Intercept) 8.2 14.9 14.9
2 GNP.deflator 8.5 13.4 13.4
3 GNP 7.8 15.7 15.7
4 Unemployed 8.3 14.5 14.5
5 Armed.Forces 8.5 14.9 14.9
6 Population 7.6 13.6 13.6
7 Year 8.2 15.6 15.6
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 least_squares_session.py
On the dose-response data the three routes agree, an intercept of 2.02 and a slope of 2.00, because the problem is
easy. Longley’s data is not: its help page calls the regression highly collinear, and NIST’s Statistical Reference
Datasets publish reference coefficients for it, to many digits, as a test of least-squares software. The
condition number: A measure of how much a matrix can magnify small relative errors in its data or in rounding. kappa(X, exact = TRUE) computes it, and t(X) %*% X has the square of the condition number of X. of X, about 23.8 million, says how much relative errors can
grow, and crossprod(X) squares it to about 5.7 × 10¹⁴. Solving that system matched the reference values to only 7.6
to 8.5 significant digits. qr.coef(qr(X), y) works on X directly, never squares the condition number, and matched
them to 13.4 to 15.7 digits.
lm() gets the same digits as qr.coef() because its help page lists QR as the only fitting method it supports; An
Introduction to R recommends lm() over the low-level functions for modelling. So write lm() for models, qr()
and qr.coef() when you need the matrix algebra yourself, and never invert t(X) %*% X.
Key idea
Forming t(X) %*% X squares the condition number and throws away digits; QR works on X itself. On Longley’s data
that was about 8 correct digits against 13 to 16.
Checkpoint
Check: Least squares with QR
Practise: predict, then check
Six short lines on a 2 × 4 matrix follow. Settle on a value for every one of them, then compare with the run.
Before you look
For m <- matrix(1:8, nrow = 2), predict m[2, 3], dim(m[, 2]), dim(m[, 2, drop = FALSE]), colSums(m),
c(1, 1) %*% m and m * 2. Then type the lines into your own R.
from session import run, version
version()
# R reads practice.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f practice.R
""")Output
# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f practice.R
> m <- matrix(1:8, nrow = 2)
> m[2, 3]
[1] 6
> dim(m[, 2])
NULL
> dim(m[, 2, drop = FALSE])
[1] 2 1
> colSums(m)
[1] 3 7 11 15
> v <- c(1, 1)
> v %*% m
[,1] [,2] [,3] [,4]
[1,] 3 7 11 15
> m * 2
[,1] [,2] [,3] [,4]
[1,] 2 6 10 14
[2,] 4 8 12 16
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 practice_session.py
The values fill by column, so row 2 of column 3 holds 6; one column drops to a vector unless drop = FALSE keeps it
as a 2 × 1 matrix; colSums() adds each column; a row vector of two 1s times the matrix adds the two rows, the same
numbers as colSums(), as a 1 × 4 matrix; and * with a single number doubles every value.
Then the self-check: it starts with your weakest idea, explains every answer as you give it, and asks a missed question again before it ends:
Practice that adapts to you
Practise this lesson
16 questions about 6 ideas of this lesson. Each answer is explained at once, and what you get right is kept in this browser.
Interview questions
Warm-up: what is the difference between * and %*% for two matrices in R?
* is element-wise: it multiplies the values in the same position, so both matrices need the same dimensions, and
the result has them too. %*% is the matrix product of linear algebra: entry [i, j] of A %*% B is the sum of row
i of A times column j of B, so the number of columns of A must equal the number of rows of B, and a 2 × 3
matrix times a 3 × 4 matrix gives a 2 × 4 one. With A <- matrix(c(2, 1, 1, 3), 2) and
B <- matrix(c(1, 0, 4, 2), 2), (A * B)[1, 2] is 1 × 4 = 4, while (A %*% B)[1, 2] is 2 × 4 + 1 × 2 = 10. Use *
for arithmetic on values, and %*%, crossprod() and solve() for linear algebra.
Model answers open with any pass: ₹79 for 1 week ($2.99), one payment that never renews.
-
Fresher to mid
A helper turns each row of a count matrix into proportions with t(apply(m, 1, function(r) r / sum(r))). It passed its tests, then a caller passed m[keep, ] with one row kept and it stopped, and for a one-column matrix it returned the wrong shape. Explain both failures and write a version that works for every matrix.
-
Mid to senior
A script computes regression coefficients as solve(t(X) %% X) %% t(X) %*% y. On one data set it stops with "system is computationally singular", and on another its coefficients differ from lm() after the sixth digit. Explain both, and say how you would compute and check the coefficients instead.
Key takeaways
- A matrix is a vector with a
dimattribute, filled by column unlessbyrow = TRUE; it holds one type, and a recycling warning frommatrix()means the values did not fit. m[i, j]follows the vector rules in each dimension, and a single row or column drops to a vector unless you writedrop = FALSE, which every function taking a matrix should.rowSums(),colSums(),rowMeans()andcolMeans()are the fast summaries;apply()runs any function along a margin and returns one column per row when the function returns a vector.- Arrays add more indices, as in an 8 × 12 × 3 array of plate readings, and keep the matrix rules for slicing and summaries.
*multiplies values in place and%*%is the matrix product; solve systems withsolve(A, b), not an inverse, and compare results withall.equal().- For least squares, QR beats the normal equations: on Longley’s data about 13 to 16 correct digits against about 8,
and
lm()uses QR.
Terms in this lesson
- array
- A vector with any number of dimensions, set by its dim attribute. A matrix is the two-dimensional case, and arr[i, j, k] picks one value of a three-dimensional array.
- BLAS
- Basic Linear Algebra Subprograms: the library of routines for matrix products and related work that R calls. Builds of R can use different BLAS libraries, which can change the last digits of results.
- column-major order
- Storing a matrix column after column: all of column 1, then all of column 2. R stores matrices and arrays this way, so m[5] is the fifth value counted down the columns.
- condition number
- A measure of how much a matrix can magnify small relative errors in its data or in rounding. kappa(X, exact = TRUE) computes it, and t(X) %*% X has the square of the condition number of X.
- matrix
- A vector with two dimensions, set by its dim attribute, so that its values sit in rows and columns. Like every atomic vector it holds values of one type.
- matrix product
- The product of linear algebra, A %*% B in R: entry [i, j] is the sum of row i of A times column j of B, so the number of columns of A must equal the number of rows of B.
- normal equations
- The system t(X) %*% X %*% b = t(X) %*% y, whose solution b is the least-squares fit of y on X. Solving it directly works with the square of the condition number of X.
- QR decomposition
- Writing a matrix X as Q %*% R, with Q orthogonal and R triangular. Least squares solved this way, as lm() and qr.coef() do, works on X itself instead of t(X) %*% X.
References
- An Introduction to R: Arrays and matrices (chapter 5) (The R Core Team)
- An Introduction to R: Linear equations and inversion (section 5.7.2) (The R Core Team)
- An Introduction to R: Least squares fitting and the QR decomposition (section 5.7.5) (The R Core Team)
- Matrices (matrix) (The R Core Team)
- Multi-way Arrays (array) (The R Core Team)
- Extract or Replace Parts of an Object (Extract): the drop argument (The R Core Team)
- Apply Functions Over Array Margins (apply) (The R Core Team)
- Form Row and Column Sums and Means (colSums) (The R Core Team)
- Matrix Multiplication (matmult) (The R Core Team)
- Matrix Crossproduct (crossprod) (The R Core Team)
- Solve a System of Equations (solve) (The R Core Team)
- The QR Decomposition of a Matrix (qr) (The R Core Team)
- Compute or Estimate the Condition Number of a Matrix (kappa) (The R Core Team)
- The Cholesky Decomposition (chol) (The R Core Team)
- Spectral Decomposition of a Matrix (eigen) (The R Core Team)
- Singular Value Decomposition of a Matrix (svd) (The R Core Team)
- Principal Components Analysis (princomp) (The R Core Team)
- Fitting Linear Models (lm): the method argument (The R Core Team)
- Longley's Economic Regression Data (longley) (The R Core Team)
- StRD Linear Least Squares Regression: Longley (National Institute of Standards and Technology (NIST))
- R Installation and Administration: BLAS (section A.3.1) (The R Core Team)
Related tools
Report a problem with this lesson
Kept only in this browser. Your Learn progress