Matrices and Arrays
|
This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production. The reference material consulted while preparing these pages is listed in the section bibliography. |
Checked against R 4.6.1 (2026-06-24) on 2026-10-20; R 4.6.2 is scheduled for 2026-10-29. Everything on this page is base R (the functions xtabs() and ftable() live in the stats package, which is attached in every session), except the last section, which uses the Matrix package, version 1.7-5, a recommended package that is installed with R. The lintr 3.4.0 default linters were used only to check the scripts.
R has no separate matrix type. A matrix is an ordinary vector that carries a dim attribute with two entries, and an array is the same thing with any number of entries. Everything else on this page follows from that sentence: the data are stored one column after another, so the first subscript moves fastest; the operators that work element by element on vectors work element by element on matrices; subsetting with one index ignores the shape; and the functions that treat the shape seriously (%*%, solve(), apply(), t()) are separate functions with their own names. This page goes through the construction of matrices and arrays, the class vector they have had since R 4.0.0, names and the drop argument, combining and reshaping, the matrix product and its relatives, the decompositions behind linear algebra, locating elements, the summaries over rows and columns, tables as arrays of counts, the BLAS and LAPACK libraries that do the numerical work, and the Matrix package for sparse matrices.
Every example below is a script that was run with Rscript on R 4.6.1 with warnings printed as they happen (options(warn = 1)), and the lines starting with #> are what it printed. The scripts pass the default lintr linters; the one idiom that a default linter dislikes, the use of the special variable .Generic inside a method in Operators on matrices, carries a # nolint: comment. Errors and warnings that a script shows on purpose are caught with tryCatch() and printed as their conditionMessage(), so the text is exactly what R says. The examples use the library scenario of the section: books, loans, shelf_today, loan_days and the book class come from the scenario.R file described on the section index, and the scripts that use them start with source("scenario.R"). Examples that do not need the shelf use small matrices with round numbers, so that every result can be checked by hand.
A matrix is a vector with dimensions
An Introduction to R defines an array as a multiply subscripted collection of data entries, a matrix as the special case with two subscripts, and says that a vector can be used as an array only if it has a dim attribute, a vector of non-negative integers (An Introduction to R: Arrays). The R Language Definition adds the rules of the representation: the content is stored in a vector in column-major order, R ensures that the length of the vector is the product of the extents, an extent may be zero, and a vector is not the same thing as a one-dimensional array, because the vector has no dim attribute and the array has one of length one (The R Language Definition: Dimensions). The order is the one of Fortran: for a dimension vector c(3, 4, 2) the data vector holds a[1, 1, 1], a[2, 1, 1], a[3, 1, 1], a[1, 2, 1] and so on, with the first subscript moving fastest and the last one slowest.
Creating a matrix
matrix(data, nrow, ncol, byrow, dimnames) builds a matrix from a vector (R documentation: matrix). The arguments and their rules are in the help page:
| Argument | Meaning |
|---|---|
|
the values, taken in order; a list is allowed and gives a list matrix; a classed object is coerced with |
|
the extents; when only one is given the other is worked out from the length of |
|
|
|
|
short |
too few values are recycled to fill the matrix; zero values give |
The first script shows the stored vector, the byrow fill, the inferred extent, the recycling and the warnings. matrix(1:3, nrow = 2, ncol = 3) is silent although the data do not line up with the columns: the warning is raised when the length of data does not divide the number of cells, here 2 * 3 = 6, and three does divide six. R 4.2.0 made matrix() warn in more cases where the length of data differs from nrow * ncol, and the environment variable _R_CHECK_MATRIX_DATA_, which R CMD check --as-cran sets, turns that warning into an error (R NEWS: Changes in R 4.2.0). The messages name the first extent that does not fit, or say that the data are longer than the matrix: data length [6] is not a sub-multiple or multiple of the number of rows [4], data length differs from size of matrix: [6 != 3 x 3].
# A matrix is a vector with a "dim" attribute; the fill is column-major
m <- matrix(1:6, nrow = 2, ncol = 3)
m
#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6
dim(m)
#> [1] 2 3
nrow(m)
#> [1] 2
ncol(m)
#> [1] 3
length(m)
#> [1] 6
typeof(m)
#> [1] "integer"
attributes(m)
#> $dim
#> [1] 2 3
#>
# byrow = TRUE fills the rows first; the stored vector does not change order
r <- matrix(1:6, nrow = 2, byrow = TRUE)
r
#> [,1] [,2] [,3]
#> [1,] 1 2 3
#> [2,] 4 5 6
as.vector(r)
#> [1] 1 4 2 5 3 6
# Only one dimension is needed: the other is inferred from the length
matrix(1:6, ncol = 2)
#> [,1] [,2]
#> [1,] 1 4
#> [2,] 2 5
#> [3,] 3 6
matrix(0, nrow = 2, ncol = 2)
#> [,1] [,2]
#> [1,] 0 0
#> [2,] 0 0
# Too few values are recycled; a non-multiple length warns
matrix(1:3, nrow = 2, ncol = 3)
#> [,1] [,2] [,3]
#> [1,] 1 3 2
#> [2,] 2 1 3
warning_text <- \(expr) tryCatch(expr, warning = \(w) conditionMessage(w))
warning_text(matrix(1:5, nrow = 2))
#> [1] "data length [5] is not a sub-multiple or multiple of the number of rows [2]"
warning_text(matrix(1:6, nrow = 4, ncol = 4))
#> [1] "data length [6] is not a sub-multiple or multiple of the number of rows [4]"
warning_text(matrix(1:6, ncol = 4))
#> [1] "data length [6] is not a sub-multiple or multiple of the number of columns [4]"
warning_text(matrix(1:6, nrow = 3, ncol = 3))
#> [1] "data length differs from size of matrix: [6 != 3 x 3]"
# Setting dim() on a vector turns it into a matrix without copying the data
v <- 1:6
dim(v) <- c(3, 2)
v
#> [,1] [,2]
#> [1,] 1 4
#> [2,] 2 5
#> [3,] 3 6
tryCatch(
dim(v) <- c(4, 2),
error = \(e) conditionMessage(e)
)
#> [1] "dims [product 8] do not match the length of object [6]"Source: R documentation: matrix
typeof(m) is "integer" and attributes(m) is a list with one element, dim: nothing else marks the object as a matrix. That is also why dim<- is the cheap way to turn a vector into a matrix. The help page for matrix() recommends it for a conversion of a vector, because it avoids duplicating the data and keeps the class of the vector, which matters for objects such as dates, and the same help page shows dimnames<- as the companion for the names. An assignment whose extents do not multiply to the length of the vector is an error: dims [product 8] do not match the length of object [6]. The manual points out the contrast with array(), which recycles a short vector and silently drops extra values, where dim<- insists on an exact match (An Introduction to R: The array() function).
Arrays of any rank
array(data, dim, dimnames) is the general constructor; matrix() is the rank-two case with a friendlier interface, and identical(matrix(1:6, 2), array(1:6, c(2, 3))) is TRUE (R documentation: array). The element a[i, j, k] of an array with extents c(2, 3, 4) is element i + 2 * (j - 1) + 6 * (k - 1) of the data vector, which the script checks. Indexing an array with one subscript per dimension follows the rules of Subsetting and Indexing, with the component of dimnames taking the place of names (The R Language Definition: Indexing matrices and arrays).
# array(data, dim, dimnames): the same column-major fill, with any rank
a <- array(1:24, dim = c(2, 3, 4))
dim(a)
#> [1] 2 3 4
a[2, 3, 4]
#> [1] 24
a[, , 1]
#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6
a[1, , ]
#> [,1] [,2] [,3] [,4]
#> [1,] 1 7 13 19
#> [2,] 3 9 15 21
#> [3,] 5 11 17 23
dim(a[1, , ])
#> [1] 3 4
# The element at (i, j, k) is at position i + 2 * (j - 1) + 6 * (k - 1)
a[2, 3, 4] == a[2 + 2 * (3 - 1) + 6 * (4 - 1)]
#> [1] TRUE
# Named dimensions
ar <- array(
1:8, dim = c(2, 2, 2),
dimnames = list(
genre = c("Fiction", "Poetry"),
decade = c("2000s", "2010s"),
format = c("print", "ebook")
)
)
ar
#> , , format = print
#>
#> decade
#> genre 2000s 2010s
#> Fiction 1 3
#> Poetry 2 4
#>
#> , , format = ebook
#>
#> decade
#> genre 2000s 2010s
#> Fiction 5 7
#> Poetry 6 8
#>
ar["Poetry", "2010s", "ebook"]
#> [1] 8
# A one-dimensional array is not a vector: it has dim and names via dimnames
one <- array(1:3, dim = 3, dimnames = list(c("a", "b", "c")))
one
#> a b c
#> 1 2 3
dim(one)
#> [1] 3
names(one)
#> [1] "a" "b" "c"
is.matrix(one)
#> [1] FALSE
# Short data is recycled; extra data is dropped
array(1:3, dim = c(2, 4))
#> [,1] [,2] [,3] [,4]
#> [1,] 1 3 2 1
#> [2,] 2 1 3 2
array(1:10, dim = c(2, 2))
#> [,1] [,2]
#> [1,] 1 3
#> [2,] 2 4
# A matrix is the rank-2 case and can be made with array()
identical(matrix(1:6, 2), array(1:6, c(2, 3)))
#> [1] TRUE
# dim<- reshapes without copying values; the order of elements stays
x <- 1:12
dim(x) <- c(3, 4)
x
#> [,1] [,2] [,3] [,4]
#> [1,] 1 4 7 10
#> [2,] 2 5 8 11
#> [3,] 3 6 9 12
dim(x) <- c(2, 2, 3)
x[, , 3]
#> [,1] [,2]
#> [1,] 9 11
#> [2,] 10 12
dim(x) <- NULL
x
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12
# An extent may be zero
dim(array(0, dim = c(0, 3)))
#> [1] 0 3Source: R documentation: array
Three details in the output matter later. The array printed with named dimnames shows each slice with a header such as , , format = print, and the names of the dimensions appear as row and column headings; those names are what apply() and marginSums() accept in place of numbers. A one-dimensional array has a dim of length one and its names() come from dimnames[[1]]; the manual says such arrays are usually treated like vectors, including when printing, but that the exceptions cause confusion, and is.matrix() is FALSE for them. Finally dim(x) <- NULL removes the shape and leaves the data vector, in the same order.
Extents: nrow, ncol, NROW and NCOL
nrow(x) and ncol(x) are dim(x)[1L] and dim(x)[2L], so they give NULL for a vector. NROW(x) and NCOL(x) treat a vector as a one-column matrix and give its length and 1, which is what code that accepts "a vector or a matrix" wants. R 4.4.0 changed NCOL(NULL) from 1 to 0, for consistency with cbind(), and NROW(NULL) was already 0 (R documentation: nrow). The script in Names and drop shows all four.
Matrices are arrays: the class vector
A matrix has had two classes since R 4.0.0. The NEWS entry for that release reads "matrix objects now also inherit from class "array"" and warns that this invalidates code which incorrectly assumes that the class of a matrix has length one (R NEWS: Changes in R 4.0.0). In the same release S3 methods for the class "array" began to be dispatched for matrix objects. The matrix help page repeats the fact in its See Also section: a matrix is the special case of a two-dimensional array and inherits(m, "array") is true for a matrix m.
The class is implicit: oldClass(m) is NULL because no class attribute is stored, and class(m) computes the vector from the dim attribute (R documentation: class). is.matrix(x) is true when x is a vector with a dim attribute of length two, which is why a data frame, which has a dim() method but no dim attribute, fails the test; is.array(x) is true for any dim attribute of positive length. The which = TRUE form of inherits() returns the position of each class in the class vector, 2 1 for c("array", "matrix"), because "matrix" comes first.
m <- matrix(1:6, nrow = 2)
a <- array(1:24, dim = c(2, 3, 4))
x <- 1:6
# Since R 4.0.0 a matrix has a class vector of length two
class(m)
#> [1] "matrix" "array"
class(a)
#> [1] "array"
class(x)
#> [1] "integer"
inherits(m, "array")
#> [1] TRUE
inherits(m, "matrix")
#> [1] TRUE
inherits(a, "matrix")
#> [1] FALSE
oldClass(m) # the implicit class is not stored as an attribute
#> NULL
inherits(m, c("array", "matrix"), which = TRUE)
#> [1] 2 1
# The type tests: is.matrix() needs exactly two dimensions
is.matrix(m)
#> [1] TRUE
is.array(m)
#> [1] TRUE
is.matrix(a)
#> [1] FALSE
is.array(a)
#> [1] TRUE
is.matrix(x)
#> [1] FALSE
is.array(array(1:3)) # a 1-d array is an array but not a matrix
#> [1] TRUE
is.matrix(data.frame(a = 1:2))
#> [1] FALSE
# The code that broke in 4.0.0: a test written for a length-one class
# (an error since R 4.2.0, a warning in 4.0.0 and 4.1.x)
old_test <- \(x) {
if (class(x) == "matrix") "a matrix" else "something else"
}
tryCatch(old_test(m), error = \(e) conditionMessage(e))
#> [1] "the condition has length > 1"
old_test(x)
#> [1] "something else"Source: R documentation: class
The last lines are the practical consequence. Code written before R 4.0.0 often tested class(x) == "matrix"; the comparison now returns two values for a matrix, and if with a condition of length greater than one is an error since R 4.2.0 (it was a warning before): the condition has length > 1. Use is.matrix(x) for "exactly two dimensions", inherits(x, "matrix") or is.array(x) for "has a shape", and never compare class(x) with a single string. The same reasoning applies to identical(class(x), "matrix"), which is now FALSE for every matrix. The attributes behind these classes are the subject of Attributes and Names.
The dispatch rule the NEWS entry describes is easy to see with a generic of your own. S3 dispatch on an object without a class attribute uses its implicit class, in the order matrix, array, the type and then the default; NextMethod() moves down that order. An "array" method therefore serves matrices, and a "matrix" method added later takes precedence over it for matrices while arrays of other ranks keep using the array method (see S3 Classes and Generics).
# S3 dispatch on a matrix tries "matrix", then "array", then the default
describe <- function(x) UseMethod("describe")
describe.array <- function(x) "an array method"
describe.default <- function(x) "the default method"
describe(matrix(1:4, 2))
#> [1] "an array method"
describe(array(1:8, c(2, 2, 2)))
#> [1] "an array method"
describe(1:4)
#> [1] "the default method"
describe.matrix <- function(x) paste("a matrix method, then", NextMethod())
describe(matrix(1:4, 2))
#> [1] "a matrix method, then an array method"
describe(array(1:8, c(2, 2, 2)))
#> [1] "an array method"Source: R NEWS: Changes in R 4.0.0
Names, dimnames and drop
The dimnames attribute is a list with one element per dimension, each either NULL or a character vector of the length of that extent. The list may itself be named, and the R Language Definition says those names are used as headings when an array is printed (The R Language Definition: Dimnames). rownames() and colnames() are the accessors for the first two elements (R documentation: colnames), and dimnames() with its replacement form sets all of them at once (R documentation: dimnames). A name that does not exist is an error for a matrix, unlike a missing name in a vector, which gives NA: subscript out of bounds. A replacement vector of the wrong length is an error too: length of 'dimnames' [1] not equal to array extent.
# Names for the rows and columns: dimnames is a list of two character vectors
visits <- matrix(
c(12, 7, 9, 4, 15, 6),
nrow = 2,
dimnames = list(c("Ana", "Ben"), c("Jan", "Feb", "Mar"))
)
visits
#> Jan Feb Mar
#> Ana 12 9 15
#> Ben 7 4 6
dimnames(visits)
#> [[1]]
#> [1] "Ana" "Ben"
#>
#> [[2]]
#> [1] "Jan" "Feb" "Mar"
#>
rownames(visits)
#> [1] "Ana" "Ben"
colnames(visits)
#> [1] "Jan" "Feb" "Mar"
# The list itself can be named: the names print as headers
named <- matrix(
1:4, nrow = 2,
dimnames = list(member = c("Ana", "Ben"), month = c("Jan", "Feb"))
)
named
#> month
#> member Jan Feb
#> Ana 1 3
#> Ben 2 4
names(dimnames(named))
#> [1] "member" "month"
# Assigning names to an existing matrix
m <- matrix(1:6, nrow = 2)
rownames(m) <- c("a", "b")
colnames(m) <- c("x", "y", "z")
m
#> x y z
#> a 1 3 5
#> b 2 4 6
dimnames(m) <- NULL
m
#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6
# Indexing by name, and the extents that dimnames must match
visits["Ana", "Feb"]
#> [1] 9
visits[, c("Jan", "Mar")]
#> Jan Mar
#> Ana 12 15
#> Ben 7 6
tryCatch(
rownames(m) <- c("a", "b", "c"),
error = \(e) conditionMessage(e)
)
#> [1] "length of 'dimnames' [1] not equal to array extent"
tryCatch(
visits["Eve", ],
error = \(e) conditionMessage(e)
)
#> [1] "subscript out of bounds"
# nrow() and ncol() need dim; NROW() and NCOL() treat a vector as one column
x <- c(5, 6, 7)
nrow(x)
#> NULL
NROW(x)
#> [1] 3
NCOL(x)
#> [1] 1
dim(x)
#> NULL
# NROW(NULL) is 0 and, since R 4.4.0, so is NCOL(NULL)
NROW(NULL)
#> [1] 0
NCOL(NULL)
#> [1] 0Source: R documentation: dimnames
drop = FALSE
Subsetting a matrix to a single row or column returns a plain vector, and an array with an extent of one loses that extent. R FAQ 7.5 says that this was determined, after much discussion, to be a feature, and gives the remedy: add drop = FALSE to the subscripting, and use it defensively in programs, because somerows <- mat[index, ] returns a vector instead of a matrix when index happens to have length one (R FAQ 7.5: Why do my matrices lose dimensions?). The function drop() is the explicit form of the same operation: it removes every extent of length one. See Subsetting and Indexing for the whole subscripting story; the script here shows the bug that drop = FALSE prevents, a helper that works for two rows and returns a vector for one.
m <- matrix(1:6, nrow = 2, dimnames = list(c("a", "b"), c("x", "y", "z")))
# A single row or column loses its dimension and becomes a plain vector
m[1, ]
#> x y z
#> 1 3 5
m[, 2]
#> a b
#> 3 4
dim(m[1, ])
#> NULL
# drop = FALSE keeps the matrix, with its names
m[1, , drop = FALSE]
#> x y z
#> a 1 3 5
m[, 2, drop = FALSE]
#> y
#> a 3
#> b 4
dim(m[1, , drop = FALSE])
#> [1] 1 3
# drop() removes extents of length one from an array on demand
one <- m[1, , drop = FALSE]
drop(one)
#> x y z
#> 1 3 5
# Code that works for three rows can break for one: a typical bug
first_two <- \(mat, rows) mat[rows, ]
dim(first_two(m, 1:2))
#> [1] 2 3
dim(first_two(m, 1))
#> NULL
dim(m[1, , drop = FALSE])
#> [1] 1 3
# The safe form takes the drop argument explicitly
first_two_safe <- \(mat, rows) mat[rows, , drop = FALSE]
dim(first_two_safe(m, 1))
#> [1] 1 3Combining: cbind() and rbind()
cbind() binds vectors and matrices side by side and rbind() stacks them (R documentation: cbind). The manual notes that the result always has matrix status, so cbind(x) and rbind(x) are the simplest ways to treat a vector explicitly as a column or a row (An Introduction to R: Forming partitioned matrices, cbind() and rbind()). The rules in the help page are worth knowing exactly:
-
If all arguments are vectors, the number of rows of the
cbind()result is the length of the longest one and shorter ones are recycled, with a warning only when the recycling is fractional (number of rows of result is not a multiple of vector length (arg 2)). -
If there are matrix arguments, they must agree on the number of rows (columns for
rbind()), otherwise the call is an error (number of rows of matrices must match (see arg 2)); vectors are recycled or cut to that length. -
The type of the result is the highest in the hierarchy raw, logical, integer, double, complex, character, list, so a character argument turns everything into text.
-
Classes of the inputs are discarded: a factor is replaced by its integer codes. Zero-length arguments, including
NULL, are ignored unless that would leave the result with no rows (columns). -
The column names of
cbind()come from thecolnamesof matrix arguments and otherwise from the names of the arguments; unnamed arguments are named by deparsing the expression, but only whendeparse.levelallows it:0never,1(the default) only when the expression is a symbol,2always.rbind()is the mirror image for row names. -
cbind()andrbind()are S3 generic, with methods for data frames; the dispatch is done internally and not withUseMethod(), so there is norbind.default. The data frame method ofcbind()is a wrapper fordata.frame(..., check.names = FALSE), andrbind()matches data frame columns by name.
# cbind() and rbind() glue vectors or matrices together
a <- 1:3
b <- c(10, 20, 30)
cbind(a, b)
#> a b
#> [1,] 1 10
#> [2,] 2 20
#> [3,] 3 30
rbind(a, b)
#> [,1] [,2] [,3]
#> a 1 2 3
#> b 10 20 30
cbind(a, b, deparse.level = 0)
#> [,1] [,2]
#> [1,] 1 10
#> [2,] 2 20
#> [3,] 3 30
# Names come from the argument names, or from the symbols given
cbind(first = a, second = b)
#> first second
#> [1,] 1 10
#> [2,] 2 20
#> [3,] 3 30
cbind(a, b + 1) # an expression gets no name at level 1
#> a
#> [1,] 1 11
#> [2,] 2 21
#> [3,] 3 31
cbind(a, b + 1, deparse.level = 2)
#> a b + 1
#> [1,] 1 11
#> [2,] 2 21
#> [3,] 3 31
# deparse.level, as in the help page example
dd <- 10
rbind(1:4, c = 2, "a++" = 10, dd, deparse.level = 0)
#> [,1] [,2] [,3] [,4]
#> 1 2 3 4
#> c 2 2 2 2
#> a++ 10 10 10 10
#> 10 10 10 10
rbind(1:4, c = 2, "a++" = 10, dd, deparse.level = 1)
#> [,1] [,2] [,3] [,4]
#> 1 2 3 4
#> c 2 2 2 2
#> a++ 10 10 10 10
#> dd 10 10 10 10
rbind(1:4, c = 2, "a++" = 10, dd, deparse.level = 2)
#> [,1] [,2] [,3] [,4]
#> 1:4 1 2 3 4
#> c 2 2 2 2
#> a++ 10 10 10 10
#> dd 10 10 10 10
# A shorter vector is recycled; fractional recycling warns
cbind(1:6, 0)
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 2 0
#> [3,] 3 0
#> [4,] 4 0
#> [5,] 5 0
#> [6,] 6 0
cbind(1:6, 1:2)
#> [,1] [,2]
#> [1,] 1 1
#> [2,] 2 2
#> [3,] 3 1
#> [4,] 4 2
#> [5,] 5 1
#> [6,] 6 2
tryCatch(
cbind(1:6, 1:4),
warning = \(w) conditionMessage(w)
)
#> [1] "number of rows of result is not a multiple of vector length (arg 2)"
# Matrices must agree on the other dimension
m <- matrix(1:6, nrow = 2)
cbind(m, 7:8)
#> [,1] [,2] [,3] [,4]
#> [1,] 1 3 5 7
#> [2,] 2 4 6 8
rbind(m, 7:9)
#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6
#> [3,] 7 8 9
tryCatch(
cbind(m, 1:3),
warning = \(w) conditionMessage(w)
)
#> [1] "number of rows of result is not a multiple of vector length (arg 2)"
tryCatch(
cbind(m, matrix(1:6, nrow = 3)),
error = \(e) conditionMessage(e)
)
#> [1] "number of rows of matrices must match (see arg 2)"
# The type of the result is the highest type of any argument
cbind(1:2, c("a", "b"))
#> [,1] [,2]
#> [1,] "1" "a"
#> [2,] "2" "b"
cbind(1:2, c(TRUE, FALSE))
#> [,1] [,2]
#> [1,] 1 1
#> [2,] 2 0
typeof(cbind(1:2, c(1.5, 2.5)))
#> [1] "double"
# Classes are discarded: a factor becomes its integer codes
f <- factor(c("x", "y"))
cbind(f, 1:2)
#> f
#> [1,] 1 1
#> [2,] 2 2
# Zero-length arguments (including NULL) are ignored
cbind(1:3, NULL)
#> [,1]
#> [1,] 1
#> [2,] 2
#> [3,] 3
cbind(NULL, NULL)
#> NULL
# Growing a matrix one row at a time works but copies it each time
acc <- NULL
for (i in 1:3) {
acc <- rbind(acc, c(i, i^2))
}
acc
#> [,1] [,2]
#> [1,] 1 1
#> [2,] 2 4
#> [3,] 3 9Source: R documentation: cbind
The loop at the end of the script is the pattern that older tutorials teach for building a matrix: start from NULL and rbind() a row per iteration. It works, and for a handful of rows it is fine, but every call allocates a new matrix and copies the rows so far, so the cost grows with the square of the number of rows. Preallocate with matrix(NA_real_, nrow = n, ncol = k) and assign m[i, ] <- row, or build a list of rows and bind them once with do.call(rbind, rows); see Performance, Profiling and Benchmarking and Memory and Copy Semantics.
Changing the shape: t(), aperm() and c()
t(x) transposes a matrix, and a plain vector is treated as a column, so t(1:3) is a one-row matrix and t(t(1:3)) is a column again (R documentation: t). For a data frame, t() first converts it with as.matrix(). aperm(a, perm) is the generalisation to any rank: the new dimension j is the old dimension perm[j], and the manual points out that aperm(A, c(2, 1)) is the transpose (An Introduction to R: Generalized transpose of an array). With no perm the dimensions are reversed, which makes aperm() equal to t() for a matrix. The dimension names, and the names of the dimensions, move with the dimensions, and perm may be given as those names. With resize = FALSE the data are rearranged but the dim attribute is kept as it was, which is rarely what you want (R documentation: aperm).
m <- matrix(1:6, nrow = 2, dimnames = list(c("a", "b"), c("x", "y", "z")))
t(m)
#> a b
#> x 1 2
#> y 3 4
#> z 5 6
dim(t(m))
#> [1] 3 2
# A vector is treated as a column, so t() gives a one-row matrix
t(1:3)
#> [,1] [,2] [,3]
#> [1,] 1 2 3
dim(t(1:3))
#> [1] 1 3
t(t(1:3))
#> [,1]
#> [1,] 1
#> [2,] 2
#> [3,] 3
# t() is aperm() with the permutation reversed
a <- array(1:24, dim = c(2, 3, 4))
dim(aperm(a))
#> [1] 4 3 2
dim(aperm(a, c(2, 1, 3)))
#> [1] 3 2 4
dim(aperm(a, c(3, 1, 2)))
#> [1] 4 2 3
# aperm() moves the elements: a[i, j, k] becomes b[k, i, j]
b <- aperm(a, c(3, 1, 2))
a[2, 3, 4] == b[4, 2, 3]
#> [1] TRUE
identical(t(m), aperm(m))
#> [1] TRUE
# Dimension names travel with the dimensions
ar <- array(
1:8, dim = c(2, 2, 2),
dimnames = list(
genre = c("Fiction", "Poetry"),
decade = c("2000s", "2010s"),
format = c("print", "ebook")
)
)
dimnames(aperm(ar, c("format", "genre", "decade")))
#> $format
#> [1] "print" "ebook"
#>
#> $genre
#> [1] "Fiction" "Poetry"
#>
#> $decade
#> [1] "2000s" "2010s"
#>
# The permutation must use every dimension exactly once
tryCatch(
aperm(a, c(1, 2)),
error = \(e) conditionMessage(e)
)
#> [1] "'perm' is of wrong length 2 (!= 3)"
tryCatch(
aperm(a, c(1, 1, 2)),
error = \(e) conditionMessage(e)
)
#> [1] "invalid 'perm' argument"
# resize = FALSE keeps the original dim and just rearranges the data
dim(aperm(a, c(2, 1, 3), resize = FALSE))
#> [1] 2 3 4Source: R documentation: aperm
The reverse operation of making an array is removing the shape. The manual says the official way to coerce an array back to a simple vector is as.vector(), and that c(X) with one argument gives a similar result as a side effect, because c() clears numeric objects of their dim and dimnames (An Introduction to R: The concatenation function, c(), with arrays). Both read the matrix column by column, which is why as.vector(matrix(1:6, nrow = 2, byrow = TRUE)) is 1 4 2 5 3 6; the script in Odds and ends shows it.
The diagonal: diag()
diag() has four distinct usages, which the help page lists, and the manual admits that making the meaning depend on the argument is confusing (R documentation: diag):
| Call | Result |
|---|---|
|
the diagonal of |
|
an identity matrix, not necessarily square |
|
the identity matrix of that size (the number is truncated with |
|
a matrix with |
The trap is the third row. The help page says that diag(x) can have unexpected effects if x is a vector that could be of length one, and recommends diag(x, nrow = length(x)) for consistent behaviour: diag(2.5) is a 2 by 2 identity matrix, not a 1 by 1 matrix holding 2.5. The replacement form diag<- sets the diagonal of an existing matrix from one value or from a vector of the right length, and it is an error to pass nrow or ncol together with a matrix: 'nrow' or 'ncol' cannot be specified when 'x' is a matrix. The identity matrix of the linear algebra books is diag(n); the function identity() is something else, the function that returns its argument, and it appears in apply() examples.
# Usage 1: a matrix gives its diagonal
m <- matrix(1:9, nrow = 3, dimnames = list(letters[1:3], letters[1:3]))
diag(m)
#> a b c
#> 1 5 9
diag(m, names = FALSE)
#> [1] 1 5 9
# Usage 2 and 3: a scalar gives an identity matrix of that size
diag(3)
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
diag(1, nrow = 2, ncol = 3)
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
diag(nrow = 2)
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1
# Usage 4: a vector gives a diagonal matrix
diag(c(2, 5, 7))
#> [,1] [,2] [,3]
#> [1,] 2 0 0
#> [2,] 0 5 0
#> [3,] 0 0 7
diag(2:1, 4)
#> [,1] [,2] [,3] [,4]
#> [1,] 2 0 0 0
#> [2,] 0 1 0 0
#> [3,] 0 0 2 0
#> [4,] 0 0 0 1
# The trap: a length-one vector is a size, not a diagonal
diag(2.5)
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1
diag(c(2.5), nrow = 1)
#> [,1]
#> [1,] 2.5
diag(5, nrow = 1)
#> [,1]
#> [1,] 5
# A safe wrapper names the size, so a length-one vector is a diagonal
scaling <- \(v) diag(v, nrow = length(v))
scaling(4)
#> [,1]
#> [1,] 4
scaling(c(4, 9))
#> [,1] [,2]
#> [1,] 4 0
#> [2,] 0 9
# Replacing the diagonal
z <- matrix(0, 3, 3)
diag(z) <- c(1, 2, 3)
z
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 2 0
#> [3,] 0 0 3
diag(z) <- 0.5
z
#> [,1] [,2] [,3]
#> [1,] 0.5 0.0 0.0
#> [2,] 0.0 0.5 0.0
#> [3,] 0.0 0.0 0.5
# Not allowed: nrow or ncol with a matrix argument
tryCatch(
diag(m, nrow = 2),
error = \(e) conditionMessage(e)
)
#> [1] "'nrow' or 'ncol' cannot be specified when 'x' is a matrix"Source: R documentation: diag
Operators on matrices
Element by element
The arithmetic, comparison and logical operators work element by element on arrays, exactly as on vectors, and keep the shape. A * B is therefore the matrix of element by element products, and A %*% B is the matrix product; the manual states the contrast in that form (An Introduction to R: Matrix multiplication), and the Language Definition warns that multiplying two matrices with * does not produce the usual matrix product (The R Language Definition: Operators). The rules for mixing arrays and vectors are in the section of An Introduction to R on the recycling rule, which admits that they are quirky and hard to find elsewhere: the expression is scanned from left to right, short vectors are extended by recycling to the size of the other operands, arrays that meet must have the same dim, a vector operand longer than an array operand is an error, and the result has the common dim of its array operands (An Introduction to R: Mixed vector and array arithmetic. The recycling rule). The Language Definition gives the compact version: for matrix plus matrix the dimensions must match, and for vector plus matrix R recycles first and then checks whether the dimensions fit (The R Language Definition: Dimensional attributes).
The script runs each rule once. m * c(1, 10) recycles the vector down the columns, because the data are column-major. Two matrices of different shapes are an error (non-conformable arrays); a vector whose length does not divide the number of cells is a warning (longer object length is not a multiple of shorter object length); and a vector longer than the matrix is an error once the warning has been silenced (dims [product 4] do not match the length of object [5]). Comparison keeps the shape too: m > 2 is a logical matrix, and which(m > 2) gives positions in the column-major vector.
The matrix product
x %*% y multiplies two matrices if they are conformable, that is if the number of columns of x equals the number of rows of y, and otherwise stops with non-conformable arguments (R documentation: matmult). A vector operand is promoted to a row or a column matrix, whichever makes the shapes conformable; the manual says that this is not always unambiguously possible. Two vectors of the same length give the inner product as a 1 by 1 matrix, which drop() turns into a number; the help page also points out that the names of a vector are not promoted to row or column names. Operands must be numeric or complex, so a character matrix gives requires numeric/complex matrix/vector arguments. The product is not commutative, and the script shows it with a matrix that swaps columns when it multiplies from the right and rows when it multiplies from the left.
m <- matrix(1:4, nrow = 2)
n <- matrix(c(10, 20, 30, 40), nrow = 2)
# * works element by element and needs the same dimensions
m * n
#> [,1] [,2]
#> [1,] 10 90
#> [2,] 40 160
m + n
#> [,1] [,2]
#> [1,] 11 33
#> [2,] 22 44
m * 2
#> [,1] [,2]
#> [1,] 2 6
#> [2,] 4 8
m * c(1, 10) # the vector is recycled down the columns
#> [,1] [,2]
#> [1,] 1 3
#> [2,] 20 40
# %*% is the matrix product, and it is not commutative
swap <- matrix(c(0, 1, 1, 0), nrow = 2)
m %*% swap # swaps the columns of m
#> [,1] [,2]
#> [1,] 3 1
#> [2,] 4 2
swap %*% m # swaps the rows of m
#> [,1] [,2]
#> [1,] 2 4
#> [2,] 1 3
m %*% n
#> [,1] [,2]
#> [1,] 70 150
#> [2,] 100 220
# Non-conformable shapes
tryCatch(
m * matrix(1:6, nrow = 3),
error = \(e) conditionMessage(e)
)
#> [1] "non-conformable arrays"
tryCatch(
m %*% matrix(1:6, nrow = 3),
error = \(e) conditionMessage(e)
)
#> [1] "non-conformable arguments"
tryCatch(
m + 1:3,
warning = \(w) conditionMessage(w)
)
#> [1] "longer object length is not a multiple of shorter object length"
# A vector longer than the matrix is an error once the warning is muted
tryCatch(
suppressWarnings(m + 1:5),
error = \(e) conditionMessage(e)
)
#> [1] "dims [product 4] do not match the length of object [5]"
tryCatch(
m %*% matrix(c("a", "b"), nrow = 2),
error = \(e) conditionMessage(e)
)
#> [1] "requires numeric/complex matrix/vector arguments"
# Vectors in %*%: promoted to a row or a column, whichever is conformable
x <- 1:3
x %*% x # inner product, a 1 x 1 matrix
#> [,1]
#> [1,] 14
drop(x %*% x)
#> [1] 14
x %*% t(x) # outer product as a 3 x 3 matrix
#> [,1] [,2] [,3]
#> [1,] 1 2 3
#> [2,] 2 4 6
#> [3,] 3 6 9
mm <- matrix(1:6, nrow = 2)
mm %*% x # x as a column
#> [,1]
#> [1,] 22
#> [2,] 28
c(1, 1) %*% mm # the vector as a row
#> [,1] [,2] [,3]
#> [1,] 3 7 11
# Comparison keeps the dimensions too
m > 2
#> [,1] [,2]
#> [1,] FALSE TRUE
#> [2,] FALSE TRUE
which(m > 2)
#> [1] 3 4Source: R documentation: matmult
The manual’s advice for a quadratic form is to avoid forming the inverse, as Linear equations shows. The help page for %*% adds that crossprod() and tcrossprod() are typically preferable for cross products, which the next section covers.
Methods for %*%: the matrixOps group generic
Since R 4.3.0 %*% is an S3 generic that belongs to a new group generic called matrixOps; the help page says methods are written for a function of two arguments named x and y and are also found for S4 classes (R NEWS: Changes in R 4.3.0). R 4.4.0 made crossprod() and tcrossprod() primitive and S3 generic as well, in the same group (R NEWS: Changes in R 4.4.0). A method for the whole group is named matrixOps.<class> and learns which operator was called from .Generic; NextMethod() continues with the default. This is the way a class such as a sparse matrix or a labelled matrix provides its own product without defining three separate methods. See S3 Group Generics and Operator Overloading for the other groups (Ops, Math, Summary).
# %*% belongs to the group generic matrixOps (R 4.3.0); crossprod() and
# tcrossprod() joined it in R 4.4.0. A method receives x and y.
counts <- function(x) structure(x, class = c("counts", class(x)))
matrixOps.counts <- function(x, y) {
cat("matrixOps method for", .Generic, "\n") # nolint: object_usage_linter.
NextMethod()
}
a <- counts(matrix(1:4, nrow = 2))
b <- matrix(1:4, nrow = 2)
unclass(a %*% b)
#> matrixOps method for %*%
#> [,1] [,2]
#> [1,] 7 15
#> [2,] 10 22
unclass(crossprod(a, b))
#> matrixOps method for crossprod
#> [,1] [,2]
#> [1,] 5 11
#> [2,] 11 25
unclass(tcrossprod(a))
#> matrixOps method for tcrossprod
#> [,1] [,2]
#> [1,] 10 14
#> [2,] 14 20
class(a %*% b)
#> matrixOps method for %*%
#> [1] "matrix" "array"Source: R documentation: groupGeneric
Cross products, outer and Kronecker products
crossprod(x, y) is t(x) %*% y, and tcrossprod(x, y) is x %*% t(y), but the transpose is not formed, so the call is faster; with one argument, crossprod(x) is t(x) %*% x, the form that appears in the normal equations of regression (R documentation: crossprod). The outer product of two arrays has the concatenated dimension vector of its operands, and its data are all possible products of an element of the first with one of the second; it is written a %o% b or outer(a, b), it is not commutative, and the multiplication can be replaced by any vectorised function of two arguments (An Introduction to R: The outer product of two arrays). The function must accept two vectors of the same length and return a vector of that length: the call outer(1:3, 1:2, \(a, b) 1) fails with dims [product 6] do not match the length of object [1] because the function returned one value, not six. kronecker(X, Y) (alias %x%) is the Kronecker product: each element of X multiplies a whole copy of Y, so the result has the product of the extents (R documentation: kronecker).
x <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 3)
y <- matrix(c(1, 0, 1, 1, 0, 1), nrow = 3)
# crossprod(x, y) is t(x) %*% y, computed without forming t(x)
crossprod(x, y)
#> [,1] [,2]
#> [1,] 4 4
#> [2,] 10 10
t(x) %*% y
#> [,1] [,2]
#> [1,] 4 4
#> [2,] 10 10
crossprod(x) # t(x) %*% x, symmetric
#> [,1] [,2]
#> [1,] 14 32
#> [2,] 32 77
tcrossprod(x) # x times its transpose
#> [,1] [,2] [,3]
#> [1,] 17 22 27
#> [2,] 22 29 36
#> [3,] 27 36 45
all.equal(crossprod(x), t(x) %*% x)
#> [1] TRUE
# For vectors, crossprod() gives a 1 x 1 matrix
crossprod(1:3)
#> [,1]
#> [1,] 14
# outer(): every pair of elements, and its %o% shorthand
1:3 %o% 1:4
#> [,1] [,2] [,3] [,4]
#> [1,] 1 2 3 4
#> [2,] 2 4 6 8
#> [3,] 3 6 9 12
outer(1:3, 1:4, "+")
#> [,1] [,2] [,3] [,4]
#> [1,] 2 3 4 5
#> [2,] 3 4 5 6
#> [3,] 4 5 6 7
outer(1:3, 1:4, \(a, b) a^b)
#> [,1] [,2] [,3] [,4]
#> [1,] 1 1 1 1
#> [2,] 2 4 8 16
#> [3,] 3 9 27 81
dim(outer(1:2, matrix(1:6, 2)))
#> [1] 2 2 3
tryCatch(
outer(1:3, 1:2, \(a, b) 1),
error = \(e) conditionMessage(e)
)
#> [1] "dims [product 6] do not match the length of object [1]"
# A table with names: outer() keeps them
loan_weeks <- c(Ana = 1, Ben = 2)
fee_cents <- c(a = 10, b = 20, c = 30)
loan_weeks %o% fee_cents
#> a b c
#> Ana 10 20 30
#> Ben 20 40 60
# kronecker(): each element of X scales a copy of Y
kronecker(diag(2), matrix(1:4, 2))
#> [,1] [,2] [,3] [,4]
#> [1,] 1 3 0 0
#> [2,] 2 4 0 0
#> [3,] 0 0 1 3
#> [4,] 0 0 2 4
diag(2) %x% matrix(1:4, 2)
#> [,1] [,2] [,3] [,4]
#> [1,] 1 3 0 0
#> [2,] 2 4 0 0
#> [3,] 0 0 1 3
#> [4,] 0 0 2 4
kronecker(1:2, matrix(1, 1, 2), FUN = "+")
#> [,1] [,2]
#> [1,] 2 2
#> [2,] 3 3Source: R documentation: outer
The options for the implementation of %*%, crossprod() and tcrossprod() are in BLAS and LAPACK. The operators table on Operators and Precedence lists %*%, %o% and %x% among the %name% operators, with their precedence.
Linear equations, determinants and conditioning
solve() and det()
solve(a, b) solves the linear system a x = b; with b missing, it returns the inverse of a. The default method is an interface to the LAPACK routines DGESV and ZGESV, the matrix must be square, and b can be a vector or a matrix (R documentation: solve). The manual’s advice is to prefer solve(a, b) to solve(a) %*% b: the inverse is rarely needed, and computing it is both inefficient and potentially unstable. It writes the quadratic form x %*% A^-1 %*% x as x %*% solve(A, x) for the same reason (An Introduction to R: Linear equations and inversion). The row and column names of the result come from the column names of a and of b, as the help page says; in the script solve(a, b) is named by the columns of a and solve(a) has the columns' names as rows and the rows' names as columns.
Two failures are distinct. A matrix that is exactly singular is detected by the factorisation itself, and the message names the pivot: Lapack routine dgesv: system is exactly singular: U[2,2] = 0. A matrix that is singular only to working precision passes the factorisation, and solve() then compares the estimated reciprocal condition number, rcond(), with the argument tol, whose default is .Machine$double.eps; if it is smaller, the error is system is computationally singular: reciprocal condition number = 2.48656e-17, where the number is the estimate for the matrix at hand and can differ in the last digits with the BLAS in use. The 12 by 12 Hilbert matrix is the standard badly conditioned example. The check is skipped when tol <= 0; the result is then dominated by rounding error. Shape errors name the dimensions: 'a' (2 x 3) must be square and 'b' (3 x 1) must be compatible with 'a' (2 x 2). The argument LINPACK of solve() is defunct. From R 4.3.0 solve() for complex inputs uses tol as well and checks for computational singularity, as it has long done for real inputs.
det(x) is a wrapper around determinant(x), which uses an LU decomposition and returns the modulus (by default on the logarithmic scale) and the sign separately; the help page adds that computing the determinant is often not what you should be doing to solve a given problem, and its logarithm = TRUE default exists because the determinant itself easily overflows or underflows (R documentation: det). A non-square argument is an error: 'x' must be a square matrix.
# A small symmetric positive definite system A x = b
a <- matrix(c(4, 2, 2, 3), nrow = 2)
b <- c(2, 1)
# solve(a, b) solves the system; solve(a) is the inverse
x <- solve(a, b)
x
#> [1] 0.5 0.0
a %*% x
#> [,1]
#> [1,] 2
#> [2,] 1
solve(a)
#> [,1] [,2]
#> [1,] 0.375 -0.25
#> [2,] -0.250 0.50
solve(a) %*% a
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1
all.equal(solve(a) %*% a, diag(2))
#> [1] TRUE
# Prefer solve(a, b) to solve(a) %*% b: it is faster and more accurate
all.equal(solve(a) %*% b, matrix(x))
#> [1] TRUE
# Names: the result takes its column names from a and b
dimnames(a) <- list(c("r1", "r2"), c("c1", "c2"))
solve(a, b)
#> c1 c2
#> 0.5 0.0
solve(a)
#> r1 r2
#> c1 0.375 -0.25
#> c2 -0.250 0.50
# The determinant, and why you rarely need it
det(a)
#> [1] 8
determinant(a)
#> $modulus
#> [1] 2.079442
#> attr(,"logarithm")
#> [1] TRUE
#>
#> $sign
#> [1] 1
#>
#> attr(,"class")
#> [1] "det"
determinant(a, logarithm = FALSE)$modulus
#> [1] 8
#> attr(,"logarithm")
#> [1] FALSE
# A singular matrix: the second row is twice the first
s <- matrix(c(1, 2, 2, 4), nrow = 2)
det(s)
#> [1] 0
tryCatch(
solve(s),
error = \(e) conditionMessage(e)
)
#> [1] "Lapack routine dgesv: system is exactly singular: U[2,2] = 0"
# The classic 3 x 3 with a dependent column
tryCatch(
solve(matrix(1:9, nrow = 3)),
error = \(e) conditionMessage(e)
)
#> [1] "Lapack routine dgesv: system is exactly singular: U[3,3] = 0"
# A matrix that is singular only to working precision: the 12 x 12 Hilbert
# matrix, the standard badly conditioned example
hilbert <- function(n) {
i <- seq_len(n)
1 / outer(i - 1, i, "+")
}
h12 <- hilbert(12)
rcond(h12)
#> [1] 2.486559e-17
tryCatch(
solve(h12),
error = \(e) conditionMessage(e)
)
#> [1] "system is computationally singular: reciprocal condition number = 2.48656e-17"
# tol = 0 skips the check; the result is then dominated by rounding error
dim(solve(h12, tol = 0))
#> [1] 12 12
# Shape errors
tryCatch(
solve(matrix(1:6, nrow = 2)),
error = \(e) conditionMessage(e)
)
#> [1] "'a' (2 x 3) must be square"
tryCatch(
solve(a, 1:3),
error = \(e) conditionMessage(e)
)
#> [1] "'b' (3 x 1) must be compatible with 'a' (2 x 2)"
tryCatch(
det(matrix(1:6, nrow = 2)),
error = \(e) conditionMessage(e)
)
#> [1] "'x' must be a square matrix"Source: R documentation: solve
norm(), rcond(), kappa() and the triangles
norm(x, type) computes a matrix norm: "O" (or "1") is the largest absolute column sum and the default, "I" the largest absolute row sum, "F" the Frobenius norm, "M" the largest absolute value of an element, and "2" the spectral norm, the largest singular value (R documentation: norm). R 4.4.0 made norm() work for complex matrices, and R 4.6.0 fixed the S4 method for type missing to use "O", as intended for years (R NEWS: Changes in R 4.6.0). The condition number of a matrix is the norm of the matrix times the norm of its inverse. rcond(x) estimates the reciprocal of the 1-norm condition number cheaply, which is what solve() uses, and kappa(x, exact = TRUE) computes the 2-norm condition number as the ratio of the largest to the smallest non-zero singular value; with the default exact = FALSE, kappa() is only a rough estimate (R documentation: kappa). The 4.4.0 release made kappa() and rcond() correct in more cases and made kappa(., norm = "2") warn that it then computes the 1-norm unless exact = TRUE.
lower.tri(x, diag = FALSE) and upper.tri(x, diag = FALSE) return logical matrices that mark the triangles, so they are indexes for reading or writing them (R documentation: lower.tri). A common use is to fill half of a symmetric matrix and mirror it with s + t(s); isSymmetric() checks the result (R documentation: isSymmetric).
m <- matrix(c(1, -2, 3, 4), nrow = 2)
norm(m) # the one norm, the largest column sum of abs values
#> [1] 7
norm(m, "I") # the largest row sum
#> [1] 6
norm(m, "F") # the Frobenius norm
#> [1] 5.477226
norm(m, "M") # the largest absolute element
#> [1] 4
norm(m, "2") # the spectral norm, the largest singular value
#> [1] 5.116673
all.equal(norm(m, "2"), svd(m)$d[1])
#> [1] TRUE
tryCatch(
norm(m, "Q"),
error = \(e) conditionMessage(e)
)
#> [1] "argument type[1]='Q' must be one of 'M','1','O','I','F' or 'E'"
# Condition numbers: rcond() estimates the reciprocal 1-norm condition number,
# kappa(exact = TRUE) is the ratio of the extreme singular values
a <- matrix(c(4, 2, 2, 3), nrow = 2)
rcond(a)
#> [1] 0.2222222
kappa(a, exact = TRUE)
#> [1] 3.866359
# lower.tri() and upper.tri() return logical masks
g <- matrix(1:9, nrow = 3)
lower.tri(g)
#> [,1] [,2] [,3]
#> [1,] FALSE FALSE FALSE
#> [2,] TRUE FALSE FALSE
#> [3,] TRUE TRUE FALSE
g[upper.tri(g)] <- 0
g
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 2 5 0
#> [3,] 3 6 9
upper.tri(g, diag = TRUE)
#> [,1] [,2] [,3]
#> [1,] TRUE TRUE TRUE
#> [2,] FALSE TRUE TRUE
#> [3,] FALSE FALSE TRUE
# Symmetrising a matrix that was filled below the diagonal
s <- matrix(0, 3, 3)
s[lower.tri(s)] <- 1:3
s <- s + t(s)
s
#> [,1] [,2] [,3]
#> [1,] 0 1 2
#> [2,] 1 0 3
#> [3,] 2 3 0
isSymmetric(s)
#> [1] TRUESource: R documentation: norm
Decompositions: qr(), chol(), eigen() and svd()
R’s matrix factorisations all return a plain list (or an object of class "qr"), and the manual introduces them in its section on matrix facilities (An Introduction to R: Matrix facilities).
The QR decomposition
qr(x) computes the QR decomposition and returns an object of class "qr" with the components qr, rank, qraux and pivot; qr.Q() and qr.R() extract the factors, and the qr.coef(), qr.fitted(), qr.resid() and qr.solve() family uses the decomposition without forming Q (R documentation: qr). The default code is the LINPACK routine with limited column pivoting and the tolerance tol = 1e-07 for detecting linear dependence of the columns; LAPACK = TRUE selects the LAPACK routine, which pivots differently and always reports full rank; qr.fitted() and qr.resid() support only the LINPACK interface. The decomposition does not assume full column rank: redundancies are found and reported in rank. qr.solve(a, b) solves a system, and for a non-square a it returns the least-squares solution, which is why the same coefficients come out of qr.solve(x, y), of the normal equations solve(crossprod(x), crossprod(x, y)) and of lm(). The QR route is the numerically safer one, because forming crossprod(x) squares the condition number; the manual says this QR-based interface is the older, low-level way to do least squares, which lm() and the facilities of Statistical Models have replaced for modelling (An Introduction to R: Least squares fitting and the QR decomposition). A rank-deficient matrix makes qr.solve() stop: singular matrix 'a' in solve. Signs in qr.R() and qr.Q() are fixed by the algorithm (the diagonal of R is negative in the script).
# A tall matrix: fit y = b0 + b1 * x by least squares with the QR decomposition
x <- cbind(1, c(1, 2, 3, 4))
y <- c(2.1, 3.9, 6.2, 7.8)
q <- qr(x)
q$rank
#> [1] 2
class(q)
#> [1] "qr"
names(q)
#> [1] "qr" "rank" "qraux" "pivot"
qr.R(q)
#> [,1] [,2]
#> [1,] -2 -5.000000
#> [2,] 0 -2.236068
round(qr.Q(q), 4)
#> [,1] [,2]
#> [1,] -0.5 0.6708
#> [2,] -0.5 0.2236
#> [3,] -0.5 -0.2236
#> [4,] -0.5 -0.6708
# The factors multiply back to x
all.equal(qr.Q(q) %*% qr.R(q), x)
#> [1] TRUE
# qr.solve() solves least squares for non-square systems
qr.solve(x, y)
#> [1] 0.15 1.94
# the same coefficients from the normal equations, which are less stable
solve(crossprod(x), crossprod(x, y))
#> [,1]
#> [1,] 0.15
#> [2,] 1.94
# and from lm(), which uses the QR decomposition internally
unname(coef(lm(y ~ x[, 2])))
#> [1] 0.15 1.94
# The pieces of the decomposition: fitted values and residuals
round(qr.fitted(q, y), 3)
#> [1] 2.09 4.03 5.97 7.91
round(qr.resid(q, y), 3)
#> [1] 0.01 -0.13 0.23 -0.11
# A rank-deficient matrix
d <- cbind(1, c(1, 2, 3), c(2, 4, 6))
qr(d)$rank
#> [1] 2
tryCatch(
qr.solve(d, c(1, 2, 3)),
error = \(e) conditionMessage(e)
)
#> [1] "singular matrix 'a' in solve"Source: R documentation: qr
Cholesky, eigenvalues and singular values
chol(x) returns the upper triangular factor R with R'R = x for a symmetric positive definite x; only the upper triangle of x is used, and an x that is not positive definite, or only positive semi-definite, is an error that names the leading minor: the leading minor of order 2 is not positive (R documentation: chol). chol2inv(r) inverts x from its Cholesky factor, which is cheaper than solve() when the factor is already at hand (R documentation: chol2inv). The argument pivot = TRUE handles positive semi-definite matrices; R 4.3.1 corrected the error or warning message for a matrix that is not positive definite, and R 4.6.0 made chol(x, pivot = TRUE) zero the trailing submatrix when the algorithm stops early, because the reconstruction recipe in the help page gave wrong results when the rank was two or more less than the dimension. A non-square input fails with 'a' must be a square matrix.
eigen(x) returns a list with values, sorted in decreasing order (by modulus for a non-symmetric matrix, whose values may be complex), and vectors, the unit-length eigenvectors as columns (R documentation: eigen). With symmetric unspecified, isSymmetric(x) decides; it is surer and typically faster to set it yourself, and with symmetric = TRUE only the lower triangle is used. only.values = TRUE skips the eigenvectors, which is the slow part for large matrices. A real matrix whose eigenvalues form a complex conjugate pair, such as a rotation, gives complex values. Missing, NaN or infinite entries are an error. The signs of the eigenvectors are not fixed by the mathematics, and the Installation and Administration manual says that the signs of eigenvectors and of singular vectors may differ between BLAS and LAPACK implementations without either being wrong (R Installation and Administration: Linear algebra); that is why the script prints eigenvalues and checks reconstructions instead of printing the vectors.
svd(x) returns d, the singular values in decreasing order, and the matrices u and v with orthonormal columns such that x = u %*% diag(d) %*% t(v) (R documentation: svd). The arguments nu and nv control how many columns of each are computed. Missing or infinite values are an error: infinite or missing values in 'x'. A non-square matrix is an error for eigen() (non-square matrix in 'eigen') but not for svd().
# Cholesky: a symmetric positive definite matrix is R'R with R upper triangular
a <- matrix(c(4, 2, 2, 3), nrow = 2)
r <- chol(a)
r
#> [,1] [,2]
#> [1,] 2 1.000000
#> [2,] 0 1.414214
all.equal(crossprod(r), a)
#> [1] TRUE
chol2inv(r)
#> [,1] [,2]
#> [1,] 0.375 -0.25
#> [2,] -0.250 0.50
all.equal(chol2inv(r), solve(a))
#> [1] TRUE
# Not positive definite: the error names the leading minor
tryCatch(
chol(matrix(c(1, 2, 2, 1), nrow = 2)),
error = \(e) conditionMessage(e)
)
#> [1] "the leading minor of order 2 is not positive"
tryCatch(
chol(matrix(1:6, nrow = 2)),
error = \(e) conditionMessage(e)
)
#> [1] "'a' must be a square matrix"
# Eigen decomposition of a symmetric matrix: real values, in decreasing order
e <- eigen(a)
round(e$values, 4)
#> [1] 5.5616 1.4384
all.equal(e$vectors %*% diag(e$values) %*% t(e$vectors), a)
#> [1] TRUE
all.equal(crossprod(e$vectors), diag(2))
#> [1] TRUE
# a symmetric matrix is detected; symmetric = FALSE forces the general code
isSymmetric(a)
#> [1] TRUE
round(eigen(a, symmetric = FALSE)$values, 4)
#> [1] 5.5616 1.4384
eigen(a, only.values = TRUE)$values
#> [1] 5.561553 1.438447
# A rotation has no real eigenvalues
rot <- matrix(c(0, 1, -1, 0), nrow = 2)
eigen(rot)$values
#> [1] 0+1i 0-1i
# Singular value decomposition
m <- matrix(c(3, 0, 0, 0, 4, 0), nrow = 3)
s <- svd(m)
s$d
#> [1] 4 3
dim(s$u)
#> [1] 3 2
dim(s$v)
#> [1] 2 2
all.equal(s$u %*% diag(s$d) %*% t(s$v), m)
#> [1] TRUE
# Missing values stop the decompositions
tryCatch(
svd(matrix(c(1, NA, 3, 4), nrow = 2)),
error = \(e) conditionMessage(e)
)
#> [1] "infinite or missing values in 'x'"
tryCatch(
eigen(matrix(1:6, nrow = 2)),
error = \(e) conditionMessage(e)
)
#> [1] "non-square matrix in 'eigen'"Source: R documentation: eigen
The manual’s three small exercises round the section off: the quadratic form without the inverse, the absolute determinant as the product of the singular values, the trace with diag(), and the distribution of the determinants of all 2 by 2 matrices with digit entries computed with outer() twice. The manual remarks that about 1 in 20 of those matrices is singular; the script gives the exact share as 0.057, and the loop-based alternative is, as the manual says, so inefficient as to be impractical.
# Three small exercises from the manual, section "Matrix facilities"
a <- matrix(c(4, 2, 2, 3), nrow = 2)
x <- c(1, 2)
# The quadratic form x' A^-1 x without forming the inverse
x %*% solve(a, x)
#> [,1]
#> [1,] 1.375
x %*% solve(a) %*% x
#> [,1]
#> [1,] 1.375
# The absolute value of the determinant from the singular values
absdet <- function(m) prod(svd(m)$d)
absdet(a)
#> [1] 8
det(a)
#> [1] 8
# The trace of a square matrix: no loop needed
tr <- function(m) sum(diag(m))
tr(a)
#> [1] 7
# Determinants of all 2 x 2 matrices with digit entries, with outer() twice
d <- outer(0:9, 0:9)
frequency <- table(outer(d, d, "-"))
length(frequency)
#> [1] 163
frequency[["0"]] / 10^4
#> [1] 0.057Locating elements: which(), row(), col() and friends
Subscripting with one index ignores the shape and uses the data vector, so m[1] usually differs a lot from m[1, ] and m[, 1], as the Language Definition points out. Two index forms respect the shape without naming each dimension. A numeric matrix with one column per dimension picks single cells, one per row: negative values are not allowed, zero rows are ignored and NA rows give NA (An Introduction to R: Index matrices). A logical matrix of the same shape picks the cells that are TRUE, as a vector. which(x, arr.ind = TRUE) goes the other way and converts the TRUE cells to row and column subscripts; arrayInd(ind, .dim) is the same conversion for any vector of positions and any dim (R documentation: which). row(m) and col(m) return matrices holding each cell’s row and column number, so row(m) > col(m) is the strict lower triangle written from first principles, and slice.index(a, d) generalises them to the d-th dimension of any array (R documentation: row, R documentation: slice.index). max.col(m) returns, for each row, the column of its largest element; its default ties.method = "random" breaks ties at random, treating entries within a relative tolerance of 1e-5 as tied, so pass "first" (the same as unname(apply(m, 1, which.max))) or "last" where the result must be reproducible (R documentation: max.col).
m <- matrix(
c(5, 8, 2, 9, 1, 7),
nrow = 3,
dimnames = list(c("a", "b", "c"), c("x", "y"))
)
m
#> x y
#> a 5 9
#> b 8 1
#> c 2 7
# which() on a matrix returns positions in the column-major vector
which(m > 6)
#> [1] 2 4 6
# arr.ind = TRUE returns row and column indices instead
which(m > 6, arr.ind = TRUE)
#> row col
#> b 2 1
#> a 1 2
#> c 3 2
which(m > 6, arr.ind = TRUE, useNames = FALSE)
#> [,1] [,2]
#> [1,] 2 1
#> [2,] 1 2
#> [3,] 3 2
# arrayInd() converts positions to subscripts for any dim
arrayInd(c(2, 4, 5), dim(m))
#> [,1] [,2]
#> [1,] 2 1
#> [2,] 1 2
#> [3,] 2 2
arrayInd(5, dim(m), dimnames(m), useNames = TRUE)
#> row col
#> b 2 2
# Matrix indexing: a two-column matrix of subscripts picks single cells
idx <- cbind(c(1, 3), c(2, 1))
idx
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 3 1
m[idx]
#> [1] 9 2
m[idx] <- 0
m
#> x y
#> a 5 0
#> b 8 1
#> c 0 7
# row() and col() give the row and column number of every cell
row(m)
#> [,1] [,2]
#> [1,] 1 1
#> [2,] 2 2
#> [3,] 3 3
col(m)
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 1 2
#> [3,] 1 2
m[row(m) == col(m)] # the diagonal, from first principles
#> [1] 5 1
m[row(m) > col(m)] # the strict lower triangle
#> [1] 8 0 7
# slice.index() generalises row() and col() to any array
a <- array(1:24, dim = c(2, 3, 4))
slice.index(a, 3)[, , 2]
#> [,1] [,2] [,3]
#> [1,] 2 2 2
#> [2,] 2 2 2
sum(a[slice.index(a, 3) == 2])
#> [1] 57
# max.col() gives the column of the largest value in each row
v <- matrix(c(3, 9, 2, 4, 1, 8), nrow = 2, byrow = TRUE)
v
#> [,1] [,2] [,3]
#> [1,] 3 9 2
#> [2,] 4 1 8
max.col(v, ties.method = "first")
#> [1] 2 3
apply(v, 1, which.max)
#> [1] 2 3Source: R documentation: which
Summaries over rows and columns: apply(), sweep() and scale()
apply()
apply(X, MARGIN, FUN, ...) calls FUN on each slice of X: MARGIN = 1 is the rows, 2 the columns, c(1, 2) every cell, and for a higher-rank array c(1, 3) keeps the first and third dimensions and applies FUN to what is left, a matrix for each pair. MARGIN can be a character vector when the dimensions are named (R documentation: apply). The shape of the result follows the help page. If every call returns a vector of length n, the result is an array with dim c(n, dim(X)[MARGIN]) when n > 1, so the answers sit in columns, one per slice, which is why a function that returns a vector per row produces the transpose of what a reader expects; when n is 1 the result is a vector for a single margin. If the calls return different lengths, the result is a list. Since R 4.1.0 simplify = FALSE always returns a list (R NEWS: Changes in R 4.1.0). Arguments after FUN are passed on to FUN, which also means that an extra argument must not partially match MARGIN or FUN. A data frame is converted with as.matrix() first, which turns it into text when one column is not numeric (see Data frames and matrices). apply() is the matrix member of the family described on The apply Family and Functionals; it is a loop in R code, not a vectorised kernel.
The first script computes the loans per member per month of 2026 with table() and applies functions over its margins. The table is an integer array, and apply(per_month, 1, which.max) returns for each member the number of the month with most loans.
source("scenario.R")
# Loans per member per month, January to May 2026
loans$member <- tolower(trimws(loans$member))
recent <- loans[
!is.na(loans$loan_date) &
loans$loan_date >= as.Date("2026-01-01") &
nzchar(loans$member),
]
month <- format(recent$loan_date, "%m")
per_month <- table(member = recent$member, month = month)
per_month
#> month
#> member 01 02 03 04 05
#> ana ruiz 1 0 1 0 2
#> chloe martin 1 0 1 1 0
#> dev patel 0 0 1 2 0
#> eva lindqvist 0 0 0 1 0
#> farid karimi 1 0 0 0 0
#> grace liu 1 2 0 0 0
#> hugo silva 0 1 0 0 1
# apply(X, MARGIN, FUN): MARGIN = 1 for rows, 2 for columns
apply(per_month, 1, sum)
#> ana ruiz chloe martin dev patel eva lindqvist farid karimi
#> 4 3 3 1 1
#> grace liu hugo silva
#> 3 2
apply(per_month, 2, sum)
#> 01 02 03 04 05
#> 4 3 3 4 3
apply(per_month, 1, max)
#> ana ruiz chloe martin dev patel eva lindqvist farid karimi
#> 2 1 2 1 1
#> grace liu hugo silva
#> 2 1
apply(per_month, 1, which.max) # the month number with most loans
#> ana ruiz chloe martin dev patel eva lindqvist farid karimi
#> 5 1 4 4 1
#> grace liu hugo silva
#> 2 2
apply(per_month, 2, range) # a matrix: FUN returns two values
#> month
#> 01 02 03 04 05
#> [1,] 0 0 0 0 0
#> [2,] 1 2 1 2 2
# The faster, dedicated functions give the same answers
rowSums(per_month)
#> ana ruiz chloe martin dev patel eva lindqvist farid karimi
#> 4 3 3 1 1
#> grace liu hugo silva
#> 3 2
colMeans(per_month)
#> 01 02 03 04 05
#> 0.5714286 0.4285714 0.4285714 0.5714286 0.4285714
all.equal(apply(per_month, 1, sum), rowSums(per_month))
#> [1] TRUE
# but the types differ: sum() of integers is an integer, rowSums() a double
typeof(apply(per_month, 1, sum))
#> [1] "integer"
typeof(rowSums(per_month))
#> [1] "double"
# A FUN that returns a vector gives a matrix with one column per row of X,
# so the result is the transpose of what you may expect
shares <- apply(per_month, 1, \(r) r / sum(r))
dim(shares)
#> [1] 5 7
round(t(shares), 2)
#> month
#> member 01 02 03 04 05
#> ana ruiz 0.25 0.00 0.25 0.00 0.5
#> chloe martin 0.33 0.00 0.33 0.33 0.0
#> dev patel 0.00 0.00 0.33 0.67 0.0
#> eva lindqvist 0.00 0.00 0.00 1.00 0.0
#> farid karimi 1.00 0.00 0.00 0.00 0.0
#> grace liu 0.33 0.67 0.00 0.00 0.0
#> hugo silva 0.00 0.50 0.00 0.00 0.5
# Proportions: of each row, of each column, of the whole table
round(proportions(per_month, 1), 2)
#> month
#> member 01 02 03 04 05
#> ana ruiz 0.25 0.00 0.25 0.00 0.50
#> chloe martin 0.33 0.00 0.33 0.33 0.00
#> dev patel 0.00 0.00 0.33 0.67 0.00
#> eva lindqvist 0.00 0.00 0.00 1.00 0.00
#> farid karimi 1.00 0.00 0.00 0.00 0.00
#> grace liu 0.33 0.67 0.00 0.00 0.00
#> hugo silva 0.00 0.50 0.00 0.00 0.50
round(proportions(per_month, 2), 2)
#> month
#> member 01 02 03 04 05
#> ana ruiz 0.25 0.00 0.33 0.00 0.67
#> chloe martin 0.25 0.00 0.33 0.25 0.00
#> dev patel 0.00 0.00 0.33 0.50 0.00
#> eva lindqvist 0.00 0.00 0.00 0.25 0.00
#> farid karimi 0.25 0.00 0.00 0.00 0.00
#> grace liu 0.25 0.67 0.00 0.00 0.00
#> hugo silva 0.00 0.33 0.00 0.00 0.33
round(proportions(per_month), 3)
#> month
#> member 01 02 03 04 05
#> ana ruiz 0.059 0.000 0.059 0.000 0.118
#> chloe martin 0.059 0.000 0.059 0.059 0.000
#> dev patel 0.000 0.000 0.059 0.118 0.000
#> eva lindqvist 0.000 0.000 0.000 0.059 0.000
#> farid karimi 0.059 0.000 0.000 0.000 0.000
#> grace liu 0.059 0.118 0.000 0.000 0.000
#> hugo silva 0.000 0.059 0.000 0.000 0.059
round(prop.table(per_month, 1), 2)
#> month
#> member 01 02 03 04 05
#> ana ruiz 0.25 0.00 0.25 0.00 0.50
#> chloe martin 0.33 0.00 0.33 0.33 0.00
#> dev patel 0.00 0.00 0.33 0.67 0.00
#> eva lindqvist 0.00 0.00 0.00 1.00 0.00
#> farid karimi 1.00 0.00 0.00 0.00 0.00
#> grace liu 0.33 0.67 0.00 0.00 0.00
#> hugo silva 0.00 0.50 0.00 0.00 0.50
marginSums(per_month, 1)
#> member
#> ana ruiz chloe martin dev patel eva lindqvist farid karimi
#> 4 3 3 1 1
#> grace liu hugo silva
#> 3 2Source: R documentation: apply
rowSums(), colSums(), rowMeans() and colMeans() are equivalent to apply() with sum or mean over the right margin, and the help page says they are much faster, and their dims argument chooses how many leading dimensions count as rows (R documentation: colSums). They accept numeric data frames and arrays of two or more dimensions, and a vector is an error: 'x' must be an array of at least two dimensions. The results are always double, which is why the script finds the values equal with all.equal() but not identical(): apply() with sum over an integer table returns integers. Missing values are omitted per column or per row when na.rm = TRUE, so the column means may not be over the same rows; select complete rows first when that matters. asplit(x, MARGIN) (R 3.6.0) splits an array into a list-like array of slices, a typed alternative with vapply() when FUN has a known return type. rowsum(x, group) is the grouped relative: it sums the rows within each group, as the last script of this page shows.
m <- matrix(1:6, nrow = 2, dimnames = list(c("a", "b"), c("x", "y", "z")))
# apply() simplifies when every result has the same length
apply(m, 1, range)
#> a b
#> [1,] 1 2
#> [2,] 5 6
# and returns a list when the lengths differ
ragged <- apply(m, 1, \(r) r[r > 3])
ragged
#> $a
#> z
#> 5
#>
#> $b
#> y z
#> 4 6
#>
class(ragged)
#> [1] "list"
# simplify = FALSE (R 4.1.0) always returns a list
apply(m, 1, range, simplify = FALSE)
#> $a
#> [1] 1 5
#>
#> $b
#> [1] 2 6
#>
# Both margins: FUN gets one element at a time, the result is a matrix
apply(m, c(1, 2), \(v) v * 10)
#> x y z
#> a 10 30 50
#> b 20 40 60
# Extra arguments are passed on after FUN
apply(m, 2, quantile, probs = 0.5)
#> x y z
#> 1.5 3.5 5.5
apply(m, 1, paste, collapse = "-")
#> a b
#> "1-3-5" "2-4-6"
# A matrix with no rows still has two columns: FUN is called on empty vectors
apply(matrix(numeric(0), nrow = 0, ncol = 2), 2, sum)
#> [1] 0 0
# asplit() splits into a list-like array of slices, which vapply() can type
rows <- asplit(m, 1)
class(rows)
#> [1] "array"
vapply(rows, sum, numeric(1))
#> a b
#> 9 12
# For sums and means the dedicated functions are the idiomatic choice
rowSums(m)
#> a b
#> 9 12
colMeans(m)
#> x y z
#> 1.5 3.5 5.5
# apply() needs a dimension, and rowSums() needs at least two
tryCatch(
apply(1:3, 1, sum),
error = \(e) conditionMessage(e)
)
#> [1] "dim(X) must have a positive length"
tryCatch(
rowSums(1:3),
error = \(e) conditionMessage(e)
)
#> [1] "'x' must be an array of at least two dimensions"Source: R documentation: apply
sweep() and scale()
sweep(x, MARGIN, STATS, FUN = "-") applies FUN between the array and a summary statistic that is recycled along the margin, which is how a column of means is subtracted from each column or a row of totals divides each row (R documentation: sweep). With the default check.margin = TRUE it warns when STATS does not recycle exactly: STATS does not recycle exactly across MARGIN. MARGIN can be a dimension name, as for apply(), since R 4.0.0, which also gave slice.index() the same ability. scale(x, center = TRUE, scale = TRUE) centres and scales the columns of a numeric matrix: with center = TRUE it subtracts the column means and with scale = TRUE divides by the column standard deviations, and when it does not centre, the divisor is the root mean square of the column, computed with n - 1 (R documentation: scale). The result is a matrix with extra attributes scaled:center and scaled:scale that record what was used, so the transformation can be undone or applied to new data with sweep(); they print after the matrix, which surprises people the first time. Both arguments can also be numeric vectors of per-column values. scale() works on columns; to standardise rows, transpose twice.
m <- matrix(
c(10, 20, 30, 12, 18, 33, 8, 25, 27),
nrow = 3,
dimnames = list(member = c("Ana", "Ben", "Eva"), week = c("w1", "w2", "w3"))
)
m
#> week
#> member w1 w2 w3
#> Ana 10 12 8
#> Ben 20 18 25
#> Eva 30 33 27
# sweep(x, MARGIN, STATS, FUN): subtract the column means
means <- colMeans(m)
sweep(m, 2, means)
#> week
#> member w1 w2 w3
#> Ana -10 -9 -12
#> Ben 0 -3 5
#> Eva 10 12 7
# STATS along the rows: divide every row by its total
sweep(m, 1, rowSums(m), "/")
#> week
#> member w1 w2 w3
#> Ana 0.3333333 0.4000000 0.2666667
#> Ben 0.3174603 0.2857143 0.3968254
#> Eva 0.3333333 0.3666667 0.3000000
# The margin can be given by the name of the dimension (R 4.0.0)
sweep(m, "week", means)
#> week
#> member w1 w2 w3
#> Ana -10 -9 -12
#> Ben 0 -3 5
#> Eva 10 12 7
# A STATS that does not fit warns
tryCatch(
sweep(m, 2, 1:2),
warning = \(w) conditionMessage(w)
)
#> [1] "STATS does not recycle exactly across MARGIN"
# scale() centres and scales the columns and records what it used
z <- scale(m)
round(z, 3)
#> week
#> member w1 w2 w3
#> Ana -1 -0.832 -1.149
#> Ben 0 -0.277 0.479
#> Eva 1 1.109 0.670
#> attr(,"scaled:center")
#> w1 w2 w3
#> 20 21 20
#> attr(,"scaled:scale")
#> w1 w2 w3
#> 10.00000 10.81665 10.44031
class(z)
#> [1] "matrix" "array"
attributes(z)[c("scaled:center", "scaled:scale")]
#> $`scaled:center`
#> w1 w2 w3
#> 20 21 20
#>
#> $`scaled:scale`
#> w1 w2 w3
#> 10.00000 10.81665 10.44031
#>
round(colMeans(z), 10)
#> w1 w2 w3
#> 0 0 0
apply(z, 2, sd)
#> w1 w2 w3
#> 1 1 1
# center and scale can be turned off or supplied
scale(m, center = TRUE, scale = FALSE)
#> week
#> member w1 w2 w3
#> Ana -10 -9 -12
#> Ben 0 -3 5
#> Eva 10 12 7
#> attr(,"scaled:center")
#> w1 w2 w3
#> 20 21 20
scale(m, center = c(20, 20, 20), scale = FALSE)
#> week
#> member w1 w2 w3
#> Ana -10 -8 -12
#> Ben 0 -2 5
#> Eva 10 13 7
#> attr(,"scaled:center")
#> [1] 20 20 20
# Without centring, scale = TRUE divides by the root mean square
attr(scale(m, center = FALSE), "scaled:scale")
#> w1 w2 w3
#> 26.45751 27.90161 26.62705
sqrt(colSums(m^2) / (nrow(m) - 1))
#> w1 w2 w3
#> 26.45751 27.90161 26.62705
# scale() works on columns: to scale the rows, transpose twice
round(t(scale(t(m))), 3)
#> week
#> member w1 w2 w3
#> Ana 0.000 1.000 -1.000
#> Ben -0.277 -0.832 1.109
#> Eva 0.000 1.000 -1.000
#> attr(,"scaled:center")
#> Ana Ben Eva
#> 10 21 30
#> attr(,"scaled:scale")
#> Ana Ben Eva
#> 2.000000 3.605551 3.000000
# The attributes survive and can undo the scaling
back <- sweep(z, 2, attr(z, "scaled:scale"), "*")
back <- sweep(back, 2, attr(z, "scaled:center"), "+")
all.equal(c(back), c(m))
#> [1] TRUESource: R documentation: scale
Proportions and margins
proportions(x, margin) divides by the sum over the margin (margin = 1 makes rows sum to one, 2 columns, none the whole table), and marginSums(x, margin) gives the sums themselves; both arrived in R 4.0.0 to replace the unfortunately named prop.table() and margin.table(), which stay as aliases for compatibility, and they accept the names of the margins (R documentation: proportions). Appearing in the script above, proportions(per_month, 1) and prop.table(per_month, 1) give the same numbers.
Tables as arrays
The last section of the matrix chapter of An Introduction to R is about frequency tables: table() takes equal-length factors and returns a k-way array of frequencies, ordered and labelled by the levels of the factors, and for one factor it gives the same counts as tapply(f, f, length) more conveniently (An Introduction to R: Frequency tables from factors). The manual builds its income-by-state example with cut(); the script uses the same pattern with the shelf, a genre factor and page-count bands. table() documents the arguments that matter in practice: useNA = "ifany" (or "always") adds an NA row or column, because NA is not a level and is otherwise left out; dnn names the dimensions; and a factor with an unused level still gets a row of zeros (R documentation: table).
source("scenario.R")
# table() on one factor: the counts follow the order of the levels
genre <- factor(tolower(trimws(books$genre)))
levels(genre)
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
table(genre)
#> genre
#> classic fiction mystery non-fiction poetry
#> 2 4 2 3 2
#> science fiction
#> 3
# the same counts with tapply(), as the manual shows
tapply(genre, genre, length)
#> classic fiction mystery non-fiction poetry
#> 2 4 2 3 2
#> science fiction
#> 3
# A second classifying factor from cut(): page-count bands
band <- cut(books$pages, breaks = c(0, 200, 300, Inf))
levels(band)
#> [1] "(0,200]" "(200,300]" "(300,Inf]"
table(genre, band)
#> band
#> genre (0,200] (200,300] (300,Inf]
#> classic 0 1 1
#> fiction 1 2 0
#> mystery 0 1 1
#> non-fiction 1 0 1
#> poetry 2 0 0
#> science fiction 0 0 3
# NA is not a level: ask for it explicitly
table(genre, band, useNA = "ifany")
#> band
#> genre (0,200] (200,300] (300,Inf] <NA>
#> classic 0 1 1 0
#> fiction 1 2 0 1
#> mystery 0 1 1 0
#> non-fiction 1 0 1 1
#> poetry 2 0 0 0
#> science fiction 0 0 3 0
# A level with no observation still gets a row or a column
genre2 <- factor(genre, levels = c(levels(genre), "drama"))
table(genre2)["drama"]
#> drama
#> 0
table(droplevels(genre2))
#>
#> classic fiction mystery non-fiction poetry
#> 2 4 2 3 2
#> science fiction
#> 3
# dnn names the dimensions without naming the arguments
table(genre, band, dnn = c("Genre", "Length"))
#> Length
#> Genre (0,200] (200,300] (300,Inf]
#> classic 0 1 1
#> fiction 1 2 0
#> mystery 0 1 1
#> non-fiction 1 0 1
#> poetry 2 0 0
#> science fiction 0 0 3
# The result is an array of counts: sum and margins
tab <- table(genre, band)
sum(tab)
#> [1] 14
dim(tab)
#> [1] 6 3
marginSums(tab, 1)
#> genre
#> classic fiction mystery non-fiction poetry
#> 2 3 2 2 2
#> science fiction
#> 3A table is an integer array with the class "table" and the dimnames (named after the arguments, or after dnn). Because the stored class is "table" and not a class vector with "array", inherits(tab, "array") is FALSE although is.matrix(tab) is TRUE for a two-way table; rowSums(), apply(), t(), indexing and %*% all treat it as the array it is. The next script crosses the genre of each book, after fixing the case and blanks of the labels, with the decade of publication. The genre column of the shelf scenario has nine spellings of six genres, so the cleaning step comes first; without it each spelling would be its own row. The books with an unknown year are left out of a two-way table unless useNA is used, and table(books$genre, is.na(books$pages)), the idiom for counting missing values by group, is the call that gives the manual’s two-way table of classifying factors.
source("scenario.R")
# Clean the genre labels, then cross-tabulate genre by decade of publication
books$genre <- tolower(trimws(books$genre))
decade <- 10 * (books$year %/% 10)
by_decade <- table(books$genre, decade)
by_decade
#> decade
#> 1810 1970 1990 2000 2010 2020
#> classic 1 0 1 0 0 0
#> fiction 0 0 0 0 1 3
#> mystery 0 0 0 0 2 0
#> non-fiction 0 0 0 1 1 1
#> poetry 0 0 0 1 0 0
#> science fiction 0 1 0 0 2 0
# A table is an integer array with a class and named dimnames
class(by_decade)
#> [1] "table"
typeof(by_decade)
#> [1] "integer"
dim(by_decade)
#> [1] 6 6
dimnames(by_decade)
#> [[1]]
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
#>
#> $decade
#> [1] "1810" "1970" "1990" "2000" "2010" "2020"
#>
is.matrix(by_decade)
#> [1] TRUE
inherits(by_decade, "array")
#> [1] FALSE
# It behaves like a matrix: indexing, margins and arithmetic work
by_decade["fiction", "2020"]
#> [1] 3
by_decade[, c("2010", "2020")]
#> decade
#> 2010 2020
#> classic 0 0
#> fiction 1 3
#> mystery 2 0
#> non-fiction 1 1
#> poetry 0 0
#> science fiction 2 0
rowSums(by_decade)
#> classic fiction mystery non-fiction poetry
#> 2 4 2 3 1
#> science fiction
#> 3
colSums(by_decade)
#> 1810 1970 1990 2000 2010 2020
#> 1 1 1 2 6 4
sum(by_decade)
#> [1] 15
# The row and column names come from the argument names: label them
table(genre = books$genre, decade = decade)
#> decade
#> genre 1810 1970 1990 2000 2010 2020
#> classic 1 0 1 0 0 0
#> fiction 0 0 0 0 1 3
#> mystery 0 0 0 0 2 0
#> non-fiction 0 0 0 1 1 1
#> poetry 0 0 0 1 0 0
#> science fiction 0 1 0 0 2 0
# NA years are left out unless asked for
sum(table(decade))
#> [1] 15
sum(table(decade, useNA = "ifany"))
#> [1] 16
table(decade, useNA = "ifany")
#> decade
#> 1810 1970 1990 2000 2010 2020 <NA>
#> 1 1 1 2 6 4 1
# Section 5.10 of the manual: a table of two classifying factors
table(books$genre, is.na(books$pages))
#>
#> FALSE TRUE
#> classic 2 0
#> fiction 3 1
#> mystery 2 0
#> non-fiction 2 1
#> poetry 2 0
#> science fiction 3 0
table(genre = books$genre, missing_pages = is.na(books$pages))
#> missing_pages
#> genre FALSE TRUE
#> classic 2 0
#> fiction 3 1
#> mystery 2 0
#> non-fiction 2 1
#> poetry 2 0
#> science fiction 3 0Source: R documentation: table
The formula interface is xtabs(). With ~ a + b it counts, like table(); with a variable on the left it sums that variable over the cells, and the result has the classes "xtabs" and "table" (R documentation: xtabs). A cell whose sum contains an NA is NA, which the print method for "xtabs" shows as a blank rather than as NA; na.rm = TRUE sums the rest. ftable() flattens a table with three or more factors into a two-dimensional layout, nesting the row factors, which is easier to read than a stack of slices (R documentation: ftable). addmargins() appends the totals, with the name Sum (R documentation: addmargins), and as.data.frame() turns a table into its long form, one row per cell, with a Freq column.
source("scenario.R")
# xtabs(): a formula interface to the same cross-tabulation
books$genre <- tolower(trimws(books$genre))
books$decade <- 10 * (books$year %/% 10)
xt <- xtabs(~ genre + decade, data = books)
xt
#> decade
#> genre 1810 1970 1990 2000 2010 2020
#> classic 1 0 1 0 0 0
#> fiction 0 0 0 0 1 3
#> mystery 0 0 0 0 2 0
#> non-fiction 0 0 0 1 1 1
#> poetry 0 0 0 1 0 0
#> science fiction 0 1 0 0 2 0
class(xt)
#> [1] "xtabs" "table"
# With a left-hand side it sums a variable instead of counting rows
xtabs(pages ~ genre, data = books) # NA prints as a blank
#> genre
#> classic fiction mystery non-fiction poetry
#> 712 545 184
#> science fiction
#> 1119
xtabs(pages ~ genre, data = books, na.rm = TRUE)
#> genre
#> classic fiction mystery non-fiction poetry
#> 712 668 545 527 184
#> science fiction
#> 1119
# ftable() flattens a three-way table to a two-dimensional layout
books$long <- books$pages > 300
ft <- ftable(xtabs(~ genre + long + decade, data = books))
ft
#> decade 1810 1970 1990 2000 2010 2020
#> genre long
#> classic FALSE 0 0 1 0 0 0
#> TRUE 1 0 0 0 0 0
#> fiction FALSE 0 0 0 0 1 2
#> TRUE 0 0 0 0 0 0
#> mystery FALSE 0 0 0 0 1 0
#> TRUE 0 0 0 0 1 0
#> non-fiction FALSE 0 0 0 0 0 1
#> TRUE 0 0 0 0 1 0
#> poetry FALSE 0 0 0 1 0 0
#> TRUE 0 0 0 0 0 0
#> science fiction FALSE 0 0 0 0 0 0
#> TRUE 0 1 0 0 2 0
# addmargins() appends the totals
addmargins(xt)
#> decade
#> genre 1810 1970 1990 2000 2010 2020 Sum
#> classic 1 0 1 0 0 0 2
#> fiction 0 0 0 0 1 3 4
#> mystery 0 0 0 0 2 0 2
#> non-fiction 0 0 0 1 1 1 3
#> poetry 0 0 0 1 0 0 1
#> science fiction 0 1 0 0 2 0 3
#> Sum 1 1 1 2 6 4 15
# as.data.frame() is the long form of a table
head(as.data.frame(xt), 4)
#> genre decade Freq
#> 1 classic 1810 1
#> 2 fiction 1810 0
#> 3 mystery 1810 0
#> 4 non-fiction 1810 0
# and matrix methods still work on it
t(xt)[, "classic"]
#> 1810 1970 1990 2000 2010 2020
#> 1 0 1 0 0 0Source: R documentation: xtabs
A three-way table is where apply() with MARGIN = c(1, 3) earns its keep. The script tabulates member by month by year and keeps two of the three dimensions; marginSums() does the same and keeps the names. Subsetting the third dimension with a name gives the member-by-month table of one year, and aperm() reorders the dimensions of a table, names included.
source("scenario.R")
# A three-way table: member by month by year, from the valid loans
loans$member <- tolower(trimws(loans$member))
ok <- !is.na(loans$loan_date) & nzchar(loans$member)
d <- loans[ok, ]
cube <- table(
member = d$member,
month = format(d$loan_date, "%m"),
year = format(d$loan_date, "%Y")
)
dim(cube)
#> [1] 8 12 2
dimnames(cube)$year
#> [1] "2025" "2026"
# Which margins does apply() keep? MARGIN = c(1, 3): member and year
by_year <- apply(cube, c(1, 3), sum)
by_year
#> year
#> member 2025 2026
#> ana ruiz 9 4
#> ben okafor 5 0
#> chloe martin 6 3
#> dev patel 3 3
#> eva lindqvist 5 1
#> farid karimi 1 1
#> grace liu 4 3
#> hugo silva 8 2
# MARGIN = 3 sums over members and months
apply(cube, 3, sum)
#> 2025 2026
#> 41 17
# MARGIN = c(2, 3) keeps the month and year: a 12 x 2 matrix
dim(apply(cube, c(2, 3), sum))
#> [1] 12 2
# marginSums() does the same for tables, and names the margins
marginSums(cube, c("member", "year"))
#> year
#> member 2025 2026
#> ana ruiz 9 4
#> ben okafor 5 0
#> chloe martin 6 3
#> dev patel 3 3
#> eva lindqvist 5 1
#> farid karimi 1 1
#> grace liu 4 3
#> hugo silva 8 2
all(marginSums(cube, c(1, 3)) == by_year)
#> [1] TRUE
# One slice of the cube is the member-by-month table of that year
cube[, , "2026"][, 1:5]
#> month
#> member 01 02 03 04 05
#> ana ruiz 1 0 1 0 2
#> ben okafor 0 0 0 0 0
#> chloe martin 1 0 1 1 0
#> dev patel 0 0 1 2 0
#> eva lindqvist 0 0 0 1 0
#> farid karimi 1 0 0 0 0
#> grace liu 1 2 0 0 0
#> hugo silva 0 1 0 0 1
# aperm() reorders the dimensions of a table; the names move along
dimnames(aperm(cube, c("year", "member", "month")))$year
#> [1] "2025" "2026"Source: R documentation: marginSums
Data frames and matrices
A data frame is not a matrix, although it has dim(), nrow() and t(). as.matrix(df) is the conversion, and its type follows the columns: all numeric columns give a numeric matrix of the highest type, all logical a logical matrix, and as soon as one column is not numeric, logical or complex the result is a character matrix, with factors converted by as.vector() and other columns by format(), so numbers become padded strings. The row names become the row names of the matrix only when they are not the automatic ones (R documentation: matrix). data.matrix(df) always returns numbers: factors become their integer codes, and since R 4.0.0 character columns are first converted to factors and then to integers (R documentation: data.matrix). That is the safe route to a numeric matrix for a mixed data frame, but the codes of a character column depend on the sorted values and mean nothing numerically. apply() calls as.matrix() on a data frame, so apply(df, 2, class) reports "character" for every column of a mixed data frame; vapply(df, class, character(1)) looks at the columns as they are. See Data Frames for the other direction.
source("scenario.R")
# as.matrix() on a data frame: all-numeric columns give a numeric matrix
nums <- books[1:4, c("year", "pages")]
as.matrix(nums)
#> year pages
#> 1 1974 387
#> 2 1813 432
#> 3 2019 214
#> 4 2021 198
typeof(as.matrix(nums))
#> [1] "integer"
# One character column and every cell becomes text; numbers are formatted
as.matrix(books[1:3, c("title", "year", "pages")])
#> title year pages
#> 1 "The Dispossessed" "1974" "387"
#> 2 "Pride and Prejudice" "1813" "432"
#> 3 " the quiet streets " "2019" "214"
# Row names of the data frame become row names of the matrix when not automatic
rownames(nums) <- c("dispossessed", "pride", "quiet", "loud")
as.matrix(nums)
#> year pages
#> dispossessed 1974 387
#> pride 1813 432
#> quiet 2019 214
#> loud 2021 198
# data.matrix() converts everything to numbers, factors to their codes
data.matrix(books[1:3, c("title", "year", "pages")])
#> title year pages
#> 1 3 1974 387
#> 2 2 1813 432
#> 3 1 2019 214
# A data frame is not a matrix, but it has dim()
is.matrix(nums)
#> [1] FALSE
dim(nums)
#> [1] 4 2
t(nums)
#> dispossessed pride quiet loud
#> year 1974 1813 2019 2021
#> pages 387 432 214 198
# apply() calls as.matrix() on a data frame first
apply(books[1:3, c("title", "year")], 2, class)
#> title year
#> "character" "character"
vapply(books[1:3, c("title", "year")], class, character(1))
#> title year
#> "character" "integer"
# The reverse: a matrix to a data frame
as.data.frame(as.matrix(nums))
#> year pages
#> dispossessed 1974 387
#> pride 1813 432
#> quiet 2019 214
#> loud 2021 198Source: R documentation: data.matrix
Odds and ends
m[] <- value assigns to every cell and keeps the dimensions, while m <- value replaces the object; m[m == 0] <- 9 assigns by a logical matrix. Matrices can hold any atomic type, and a list gives a list matrix, whose cells hold values of different types and print by type and length. The identity matrix is diag(n), and the script shows that identity() is the function that returns its argument. Two matrices with the same values and different dim are not identical(). The script ends with rowsum(), which sums rows within groups.
m <- matrix(1:6, nrow = 2)
# m[] <- keeps the dimensions; m <- replaces the object
m[] <- 0L
m
#> [,1] [,2] [,3]
#> [1,] 0 0 0
#> [2,] 0 0 0
m[m == 0L] <- 9L # logical indexing returns a vector, assigns in place
m
#> [,1] [,2] [,3]
#> [1,] 9 9 9
#> [2,] 9 9 9
m[m > 5L]
#> [1] 9 9 9 9 9 9
# c() drops the dimensions and reads column by column
c(matrix(1:6, nrow = 2))
#> [1] 1 2 3 4 5 6
as.vector(matrix(1:6, nrow = 2, byrow = TRUE))
#> [1] 1 4 2 5 3 6
# Matrices of other types
matrix(c("a", "b", "c", "d"), nrow = 2)
#> [,1] [,2]
#> [1,] "a" "c"
#> [2,] "b" "d"
matrix(c(TRUE, FALSE, NA, TRUE), nrow = 2)
#> [,1] [,2]
#> [1,] TRUE NA
#> [2,] FALSE TRUE
matrix(list(1, "a", TRUE, NULL), nrow = 2)
#> [,1] [,2]
#> [1,] 1 TRUE
#> [2,] "a" NULL
# identity() is the function x -> x; diag(n) is the identity matrix
identity(5)
#> [1] 5
diag(2)
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1
# Same values, different dimensions: not identical
identical(matrix(1:4, 2), matrix(1:4, 1))
#> [1] FALSE
# rowsum() groups the rows by a key; compare with rowSums()
rowsum(c(1, 2, 3, 4), group = c("a", "b", "a", "b"))
#> [,1]
#> a 4
#> b 6Source: R documentation: rowsum
BLAS and LAPACK
R does not implement matrix algebra itself. The Installation and Administration manual states in its appendix on linear algebra that the routines make use of BLAS (Basic Linear Algebra Subprograms) and, for most of them, of LAPACK (Linear Algebra PACKage), that the R sources contain reference Fortran implementations of both, and that they can be replaced by external libraries, usually tuned for specific CPUs (R Installation and Administration: Linear algebra). %*%, crossprod(), solve(), qr(LAPACK = TRUE), chol(), eigen() and svd() end up there. The practical consequences are in the same appendix:
-
An external BLAS gives results that are not identical to the reference implementation. Some differences are benign, such as the signs of singular vectors and eigenvectors, but the optimised routines can be less accurate and, for LAPACK, older, with fewer corrections.
-
R relies on ISO/IEC 60559 (IEEE 754) arithmetic, and a library that assumes a term with a zero factor is zero can break that, since
x * 0can beNaN. -
Multi-threaded BLAS libraries can use more CPU time, and even more elapsed time, than a single-threaded one, and BLAS and LAPACK calls inside threaded code (OpenMP sections in packages such as
mgcv) may be unsafe with some external libraries. See Parallel and Asynchronous Computing. -
BLAS has been unchanged for decades, until LAPACK 3.12.1 (January 2025) added routines such as
dgemmtr. R 4.5.0 updated its bundled sources to that release, stated that it is expected to be the last update of the bundled BLAS and LAPACK, and that building against external libraries will be required in future; BLAS libraries are therefore no longer freely interchangeable (R NEWS: Changes in R 4.5.0). -
The choice is made when R is built (
--with-blas,--with-lapack) or, on Debian-style systems with the alternatives mechanism and with FlexiBLAS, switched afterwards; see the BLAS and LAPACK subsections of the appendix (R Installation and Administration: BLAS, R Installation and Administration: LAPACK).
From inside R, La_version() reports the version of LAPACK in use at run time (R documentation: La_version), La_library() the shared object file that contains it (R documentation: La_library), and extSoftVersion() a named character vector with the versions of third-party libraries, among them the BLAS entry with the name of the library file (R documentation: extSoftVersion). The file names are machine specific, so the script prints the LAPACK version, which is the number the manual’s text is about, and shows that La_library() returns a string without printing it. The file is typically libRlapack.so (libRlapack.dylib on macOS) for the internal LAPACK and a vendor name for an external one, and the detection does not work on Windows or for the Accelerate framework on macOS. sessionInfo() includes La_version() since R 4.3.0, and since R 4.5.0 it still reports it when La_library() is empty (R documentation: sessionInfo). The version printed below, 3.12.0, is that of the container this page was checked in, and on your machine it can differ.
The option matprod selects how %*%, crossprod() and tcrossprod() are computed for double and complex vectors: "default" uses BLAS but falls back to a plain three-loop algorithm when the input may contain NaN or Inf, so that they propagate correctly; "internal" always uses the three-loop code, with a long double accumulator where available; "blas" always uses BLAS without checks, which the help page says needs extreme care; and "default.simd" is experimental (R documentation: options).
# Which LAPACK is behind solve(), qr(LAPACK = TRUE), eigen() and svd()?
La_version()
#> [1] "3.12.0"
is.character(La_library()) # the value is a path: not printed here
#> [1] TRUE
# extSoftVersion() reports the third-party libraries; BLAS is among them
names(extSoftVersion())
#> [1] "zlib" "bzlib" "xz" "libdeflate" "zstd"
#> [6] "PCRE" "ICU" "TRE" "iconv" "readline"
#> [11] "BLAS"
# How %*%, crossprod() and tcrossprod() treat NaN and Inf
getOption("matprod")
#> [1] "default"
# The results agree up to rounding whichever BLAS is used
a <- matrix(c(4, 2, 2, 3), nrow = 2)
all.equal(solve(a) %*% a, diag(2))
#> [1] TRUESource: R documentation: La_version
Sparse matrices: the Matrix package
A shelf of 16 books and 8 members is small, but the table of members by books of a real library is mostly zeros, and a dense matrix stores every zero. The recommended package Matrix, version 1.7-5 on this page, provides S4 classes for sparse and dense matrices with the same operators: %*%, t(), crossprod(), solve(), chol() and qr() have methods, and Matrix::sparseMatrix(i, j, x, dims, dimnames) builds a column-compressed sparse matrix from the row and column positions of the non-zero entries, adding up entries that repeat the same position (Matrix package documentation: sparseMatrix). Matrix::Matrix(data, nrow, ncol, sparse) is the friendlier constructor, which chooses a class from the content: a diagonal matrix becomes a "ddiMatrix", and sparse = TRUE asks for a sparse class (Matrix package documentation: Matrix). Class names such as "dgCMatrix" encode the content: d for double, g for general, and C for column-compressed sparse (Matrix package documentation: Matrix-class). A sparse matrix prints a dot for each structural zero. The script counts the loans by member and ISBN, takes the products tcrossprod(borrowed) that stay sparse and symmetric ("dsCMatrix"), and converts a corner to a base matrix with as.matrix(). The package is the right tool when the data are large and mostly zero; for the small dense matrices of this page, base R is enough. The package page is Matrix on CRAN.
source("scenario.R")
# Who borrowed which book? Count loans by member and ISBN as a sparse matrix
loans$member <- tolower(trimws(loans$member))
known <- loans$isbn %in% normalise_isbn(books$isbn)
d <- loans[nzchar(loans$member) & known, ]
who <- factor(d$member)
what <- factor(d$isbn)
borrowed <- Matrix::sparseMatrix(
i = as.integer(who),
j = as.integer(what),
x = 1,
dimnames = list(levels(who), levels(what))
)
class(borrowed)[[1]]
#> [1] "dgCMatrix"
methods::is(borrowed, "sparseMatrix")
#> [1] TRUE
dim(borrowed)
#> [1] 8 15
Matrix::nnzero(borrowed) # stored non-zero entries (duplicates summed)
#> [1] 46
sum(borrowed) # the number of loans used
#> [1] 58
# A few cells, shown with dots for the zeros
borrowed[1:4, 1:6]
#> 4 x 6 sparse Matrix of class "dgCMatrix"
#> 9780001000100 9780002000291 9780003000313 9780004000404
#> ana ruiz 2 . . 1
#> ben okafor . . . .
#> chloe martin 2 1 . 1
#> dev patel 1 1 . .
#> 9780005000571 9780006000624
#> ana ruiz . 2
#> ben okafor . .
#> chloe martin 1 .
#> dev patel . 1
# Products stay sparse: cell (i, j) adds up count_i * count_j over the books
# that members i and j both borrowed
shared <- Matrix::tcrossprod(borrowed)
class(shared)[[1]]
#> [1] "dsCMatrix"
dim(shared)
#> [1] 8 8
as.matrix(shared)[1:3, 1:3]
#> ana ruiz ben okafor chloe martin
#> ana ruiz 28 4 13
#> ben okafor 4 6 1
#> chloe martin 13 1 13
# Matrix() picks a class from the contents, sparse = TRUE asks for sparse
small <- Matrix::Matrix(c(0, 0, 3, 0, 5, 0), nrow = 3, sparse = TRUE)
small
#> 3 x 2 sparse Matrix of class "dgCMatrix"
#>
#> [1,] . .
#> [2,] . 5
#> [3,] 3 .
class(small)[[1]]
#> [1] "dgCMatrix"
class(Matrix::Matrix(diag(3)))[[1]] # a diagonal matrix gets its own class
#> [1] "ddiMatrix"Common errors and warnings
The messages below were produced by the examples of this page, in R 4.6.1.
| Message | Typical cause | Fix |
|---|---|---|
|
|
make the data fill the matrix, or recycle on purpose with a length that divides it |
|
data longer than the matrix, or a size that is not a multiple |
check |
|
|
|
|
|
supply one name per row or column |
|
a row or column number beyond the extent, or a name that is not in |
check |
|
|
compare |
|
|
make the dimensions equal, or use a vector |
|
a vector that does not divide the number of cells |
recycle with a divisor of the length |
|
|
use a vector of the right length |
|
|
make the matrices conformable |
|
|
pass each dimension number once |
|
|
|
|
|
convert with |
|
|
check the rank with |
|
|
rescale, regularise or use a decomposition suited to the problem; |
|
|
check the shapes |
|
|
|
|
|
check symmetry and definiteness, or try |
|
|
use |
|
|
drop dependent columns or use |
|
|
use |
|
|
pass one value per row or column |
|
|
|
|
|
use |
|
|
|
Version notes
| Version | Change |
|---|---|
R 3.6.0 |
|
R 4.0.0 |
Matrices inherit from class |
R 4.1.0 |
|
R 4.2.0 |
|
R 4.3.0 |
|
R 4.3.1 |
|
R 4.4.0 |
|
R 4.5.0 |
The bundled BLAS and LAPACK are updated to LAPACK 3.12.1 and are announced as the last update; BLAS libraries are no longer interchangeable because of the added |
R 4.6.0 |
|
Where older tutorials differ
| Older tutorial | Current practice | Why | Version |
|---|---|---|---|
|
|
the class vector has length two, and a condition of length two is an error |
R 4.0.0, R 4.2.0 |
|
|
better names and named margins; the old names remain as aliases |
R 4.0.0 |
|
|
the inverse is rarely needed and computing it is less accurate and slower |
always |
|
|
no transpose is formed |
always |
|
|
the help page says the dedicated functions are a lot faster |
always |
|
|
a single row or column loses its dimension (R FAQ 7.5) |
always |
growing a matrix with |
preallocate and assign, or |
every |
always |
|
|
a length-one |
always |
|
|
the manual says |
always |
|
call without the argument |
the |
R 3.1.0, R 4.1.0 |
"install ATLAS or MKL and everything is faster" |
check |
alternative BLAS and LAPACK libraries do not give identical results and may be unsafe in threaded code |
always |
For the full list of changes by release see What’s New and Migrating Older Code.
Where to go next
-
Atomic Vectors and Types for the vectors that a matrix is made of.
-
Subsetting and Indexing for the full rules of
[,[[and index matrices. -
Attributes and Names for
dim,dimnamesand the other attributes. -
Factors for the classifying variables that
table()cross-tabulates. -
Data Frames for the table-shaped type that keeps columns of different types.
-
The apply Family and Functionals for
lapply(),vapply()andMap()besideapply(). -
Numbers and Arithmetic for floating-point behaviour that decides what
solve()can do. -
Statistical Models for
lm(), which uses the QR decomposition of the model matrix. -
Performance, Profiling and Benchmarking for measuring matrix code before changing it.
-
Interfacing C, C++ and Fortran for calling BLAS and LAPACK routines from compiled code.
References
-
An Introduction to R: Arrays and matrices, in particular Arrays, Index matrices, The array() function, The outer product of two arrays, Generalized transpose of an array, Matrix facilities, Forming partitioned matrices and Frequency tables from factors
-
The R Language Definition: Dimensions, Dimnames, Dimensional attributes and Indexing matrices and arrays
-
R Installation and Administration: Linear algebra (BLAS and LAPACK)
-
R FAQ 7.5: Why do my matrices lose dimensions?
-
R NEWS: Changes in R 3.6.0, 4.0.0, 4.1.0, 4.2.0, 4.3.0, 4.3.1, 4.4.0, 4.5.0 and 4.6.0
-
matrix (including as.matrix and is.matrix)