Matrices are the backbone of statistical computing. Whether you’re aware of it or not, every time you fit a linear regression model, perform matrix decomposition, or work with multivariate statistics, you’re relying on matrix operations happening behind the scenes.

I’ve found that understanding matrix operations deeply;not just using them;makes you a much better data scientist. It helps you understand why certain algorithms work, how to optimize performance, and how to debug numerical issues. In my work with large datasets and complex statistical models, solid matrix knowledge has saved me hours of debugging time.

Whether you’re working with spatial data, solving systems of equations, implementing machine learning algorithms, or optimizing statistical computations, mastering matrix operations in R is non-negotiable. R provides powerful built-in functions for matrix operations and integrates seamlessly with linear algebra libraries.

This comprehensive guide covers all essential matrix operations and linear algebra concepts with practical R implementations tested in R 4.0+.

Matrix Basics

Creating Matrices

# Create matrix from vector
m <- matrix(1:12, nrow = 3, ncol = 4)
print(m)
#      [,1] [,2] [,3] [,4]
# [1,]    1    4    7   10
# [2,]    2    5    8   11
# [3,]    3    6    9   12

# Matrix with specific values
m2 <- matrix(0, nrow = 2, ncol = 3)  # Matrix of zeros
m3 <- matrix(1, nrow = 3, ncol = 3)  # Matrix of ones

# Combine vectors as rows or columns
v1 <- c(1, 2, 3)
v2 <- c(4, 5, 6)
by_row <- rbind(v1, v2)
by_col <- cbind(v1, v2)

# Check matrix properties
nrow(m)   # Number of rows
ncol(m)   # Number of columns
dim(m)    # Dimensions

Identity Matrix

# Create identity matrix
I3 <- diag(3)  # 3x3 identity matrix
print(I3)
#      [,1] [,2] [,3]
# [1,]    1    0    0
# [2,]    0    1    0
# [3,]    0    0    1

# Create diagonal matrix from vector
D <- diag(c(2, 4, 6))  # Diagonal matrix with 2, 4, 6
print(D)
#      [,1] [,2] [,3]
# [1,]    2    0    0
# [2,]    0    4    0
# [3,]    0    0    6

# Extract diagonal from matrix
diag(m)  # Gets diagonal elements

Matrix Operations

Element-wise Operations

A <- matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2)
B <- matrix(c(5, 6, 7, 8), nrow = 2, ncol = 2)

# Element-wise multiplication
result <- A * B
print(result)
#      [,1] [,2]
# [1,]    5   21
# [2,]   12   32

# Element-wise addition
A + B
#      [,1] [,2]
# [1,]    6   10
# [2,]    8   12

# Element-wise subtraction
A - B
#      [,1] [,2]
# [1,]   -4   -4
# [2,]   -4   -4

# Element-wise division
A / B
#           [,1]       [,2]
# [1,] 0.2000000 0.4285714
# [2,] 0.3333333 0.5000000

# Element-wise power
A ^ 2
#      [,1] [,2]
# [1,]    1    9
# [2,]    4   16

Matrix Multiplication

# Matrix multiplication using %*%
A <- matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2)
B <- matrix(c(5, 6, 7, 8), nrow = 2, ncol = 2)

result <- A %*% B  # Matrix multiplication
print(result)
#      [,1] [,2]
# [1,]   19   22
# [2,]   43   50

# Verify: A[1,1] * B[1,1] + A[1,2] * B[2,1] = 1*5 + 3*6 = 23

# Matrix-vector multiplication
v <- c(2, 3)
A %*% v
#      [,1]
# [1,]   11
# [2,]   16

# Verify: 1*2 + 3*3 = 11, 2*2 + 4*3 = 16

# Transpose
t(A)  # Transpose matrix
#      [,1] [,2]
# [1,]    1    2
# [2,]    3    4

Outer and Cross Products

# Dot product (inner product)
u <- c(1, 2, 3)
v <- c(4, 5, 6)

dot_product <- u %*% v  # Same as sum(u * v)
print(dot_product)
# [1] 32

# Manual verification: 1*4 + 2*5 + 3*6 = 32

# Outer product
outer_product <- u %o% v
print(outer_product)
#      [,1] [,2] [,3]
# [1,]    4    5    6
# [2,]    8   10   12
# [3,]   12   15   18

# Alternative: outer(u, v, "*")

# Cross product (3D vectors only)
u3 <- c(1, 0, 0)
v3 <- c(0, 1, 0)

cross_product <- function(u, v) {
  c(u[2]*v[3] - u[3]*v[2],
    u[3]*v[1] - u[1]*v[3],
    u[1]*v[2] - u[2]*v[1])
}

cross <- cross_product(u3, v3)
print(cross)  # [1] 0 0 1

Matrix Properties and Decompositions

Matrix Inverse

# Create invertible matrix
A <- matrix(c(4, 7, 2, 6), nrow = 2, ncol = 2)

# Calculate inverse
A_inv <- solve(A)
print(A_inv)

# Verify: A %*% A_inv ≈ I
result <- A %*% A_inv
print(result)
#      [,1] [,2]
# [1,]    1    0
# [2,]    0    1

# Check determinant (must be non-zero for invertibility)
det(A)  # [1] 10

Determinant and Trace

A <- matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2)

# Determinant
det(A)  # [1] -2

# Trace (sum of diagonal elements)
trace <- function(x) sum(diag(x))
trace(A)  # [1] 5

# For identity matrix
trace(diag(3))  # [1] 3

Eigenvalues and Eigenvectors

# Create symmetric matrix
A <- matrix(c(4, -2, -2, 4), nrow = 2, ncol = 2)

# Calculate eigenvalues and eigenvectors
eig <- eigen(A)
print(eig)
# $values
# [1] 6 2
# $vectors
#            [,1]       [,2]
# [1,] -0.7071068  0.7071068
# [2,] -0.7071068 -0.7071068

# Access eigenvalues and eigenvectors
eigenvalues <- eig$values
eigenvectors <- eig$vectors

# Verify: A %*% v ≈ λ * v
v1 <- eigenvectors[, 1]
lambda1 <- eigenvalues[1]

A %*% v1  # Should equal lambda1 * v1
lambda1 * v1

Solving Linear Systems

Basic Linear Systems

# Solve Ax = b
A <- matrix(c(2, 1, 1, 3), nrow = 2, ncol = 2)
b <- c(5, 6)

# Solution
x <- solve(A, b)
print(x)  # [1] 2.6 -0.2

# Verify: A %*% x should equal b
A %*% x
# [1] 5 6 ✓

# Example: System of equations
# 2x + y = 5
# x + 3y = 6

# Solution: x = 2.6, y = -0.2

Least Squares Regression

# Data
X <- cbind(1, c(1, 2, 3, 4, 5))  # Design matrix with intercept
y <- c(2.1, 3.9, 6.2, 7.8, 10.1)  # Response

# Least squares solution: (X'X)^-1 X'y
beta <- solve(t(X) %*% X) %*% t(X) %*% y
print(beta)
# Intercept: 0.02
# Slope: 1.98

# Or use lm() function
model <- lm(y ~ X[, 2])
coef(model)

Matrix Decompositions

QR Decomposition

A <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 3, ncol = 2)

# QR decomposition
qr_decomp <- qr(A)
Q <- qr.Q(qr_decomp)  # Orthogonal matrix
R <- qr.R(qr_decomp)  # Upper triangular matrix

# Verify: A = Q %*% R
A
Q %*% R

# Check orthogonality: Q'Q = I
t(Q) %*% Q

Singular Value Decomposition (SVD)

A <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)

# SVD
svd_decomp <- svd(A)
U <- svd_decomp$u      # Left singular vectors
s <- svd_decomp$d      # Singular values
V <- svd_decomp$v      # Right singular vectors

# Reconstruct matrix
A_reconstructed <- U %*% diag(s) %*% t(V)
print(A_reconstructed)  # Should equal A

# Rank of matrix
rank <- length(s[s > 1e-10])
print(rank)

Cholesky Decomposition

# For positive definite matrix
A <- matrix(c(4, 2, 2, 3), nrow = 2, ncol = 2)

# Check if positive definite
all(eigen(A)$values > 0)  # [1] TRUE

# Cholesky decomposition: A = L L'
L <- chol(A)  # Returns upper triangular (transpose for lower)
print(L)

# Verify: A = t(L) %*% L
t(L) %*% L  # Should equal A

Norms and Distances

v <- c(3, 4)

# Euclidean norm (L2 norm)
euclidean_norm <- sqrt(sum(v^2))
print(euclidean_norm)  # [1] 5

# Using norm function from base R
norm(v, type = "2")  # Euclidean norm

# Manhattan norm (L1 norm)
manhattan_norm <- sum(abs(v))
print(manhattan_norm)  # [1] 7
norm(v, type = "1")

# Frobenius norm for matrices
A <- matrix(1:6, nrow = 2, ncol = 3)
frob_norm <- sqrt(sum(A^2))
print(frob_norm)
norm(A, type = "F")  # Frobenius norm

# Distance between vectors
u <- c(1, 2, 3)
v <- c(4, 5, 6)
distance <- sqrt(sum((u - v)^2))
print(distance)

Matrix Subsetting and Operations

A <- matrix(1:12, nrow = 3, ncol = 4)

# Extract rows and columns
A[1, ]        # First row
A[, 2]        # Second column
A[1:2, 2:3]   # Rows 1-2, columns 2-3
A[c(1, 3), ]  # Rows 1 and 3

# Replace elements
A[1, 1] <- 99
A[2, ] <- c(20, 30, 40, 50)

# Matrix binding
B <- matrix(13:16, nrow = 2, ncol = 2)
rbind(A, B)  # Row bind
cbind(A, B)  # Column bind

# Apply operations
rowSums(A)      # Sum each row
colSums(A)      # Sum each column
rowMeans(A)     # Mean of each row
colMeans(A)     # Mean of each column

Practical Applications

Solving Linear Regression

# Data
set.seed(123)
x <- 1:10
y <- 2*x + rnorm(10, 0, 2)

# Using matrix approach
X <- cbind(1, x)  # Design matrix
y <- as.matrix(y)

# Coefficients: (X'X)^-1 X'y
beta <- solve(t(X) %*% X) %*% t(X) %*% y
print(beta)

# Predictions
y_pred <- X %*% beta

# Compare with lm()
model <- lm(y ~ x)
coef(model)

Principal Component Analysis (PCA)

# Data
data <- matrix(rnorm(100), nrow = 20, ncol = 5)

# Center the data
data_centered <- scale(data, center = TRUE, scale = FALSE)

# Covariance matrix
cov_matrix <- (t(data_centered) %*% data_centered) / (nrow(data) - 1)

# Eigendecomposition
eig <- eigen(cov_matrix)
eigenvalues <- eig$values
eigenvectors <- eig$vectors

# Principal components
pc_scores <- data_centered %*% eigenvectors

# Variance explained
var_explained <- eigenvalues / sum(eigenvalues)
cumsum(var_explained)  # Cumulative variance

Transformations and Rotations

# 2D rotation matrix
theta <- pi / 4  # 45 degrees

rotation_matrix <- matrix(
  c(cos(theta), -sin(theta),
    sin(theta), cos(theta)),
  nrow = 2, ncol = 2
)

# Point to rotate
point <- c(1, 0)

# Rotated point
rotated <- rotation_matrix %*% point
print(rotated)

# Scaling matrix
scale_matrix <- diag(c(2, 0.5))  # Scale x by 2, y by 0.5
scaled <- scale_matrix %*% point

Best Practices

  1. Check matrix properties - Verify dimensions before operations
  2. Use appropriate operations - Element-wise () vs matrix multiplication (%%)
  3. Check invertibility - Verify determinant ≠ 0 before solving
  4. Handle singular matrices - Use SVD or pseudo-inverse for singular systems
  5. Numerical stability - Be aware of rounding errors in large operations
  6. Use efficient functions - Leverage built-in functions for speed
  7. Vectorize operations - Avoid loops with matrix operations
  8. Document decompositions - Record which method used for reproducibility

Troubleshooting Matrix Operations

Issue: “Non-conformable arguments” Error

Problem: You get an error like: Error in A %*% B : non-conformable arguments

Cause: Matrix dimensions don’t match for multiplication. For A %*% B, A’s columns must equal B’s rows.

Solution:

# Check dimensions first
A <- matrix(1:6, nrow = 2, ncol = 3)  # 2x3 matrix
B <- matrix(1:6, nrow = 2, ncol = 3)  # 2x3 matrix

# ❌ WRONG - Can't multiply 2x3 by 2x3
A %*% B  # Error!

# ✅ RIGHT - Transpose B to make it 3x2
A %*% t(B)  # Now 2x3 times 3x2 = 2x2

Issue: Unexpected results from inverse or solve()

Problem: solve(A) returns strange values or errors.

Causes & Solutions:

  1. Singular matrix - Matrix is not invertible (det = 0) 
    det(A)  # Check if this is 0 or very close to 0
  2. Near-singular matrix - Very small determinant causes numerical instability 
    # Use pseudo-inverse instead
A_pseudo <- MASS::ginv(A)
  3. Solving Ax=b incorrectly - Use solve(A, b) not solve(A) %*% b 
    # ✅ CORRECT - More accurate
x <- solve(A, b)

Issue: Eigenvalue/eigenvector confusion

Problem: Not sure which is which or getting unexpected values.

Explanation:

result <- eigen(A)
result$values   # Eigenvalues (scalars)
result$vectors  # Eigenvectors (columns are vectors)

# Verify: A %*% v = λ * v
lambda <- result$values[1]
v <- result$vectors[, 1]
all.equal(A %*% v, lambda * v)  # Should be TRUE

Issue: Performance issues with large matrices

Problem: Matrix operations running very slowly.

Solutions:

  1. Use efficient operations - crossprod(A) instead of t(A) %*% A
  2. Check data type - Use numeric, not character matrices
  3. Consider sparse matrices - Use Matrix package for sparse data
  4. Vectorize loops - Avoid for loops with matrix operations

Issue: NA or NaN values in results

Problem: Unexpected NA or NaN values appearing in matrix calculations.

Causes:

  1. Singular matrix - Some operations can’t be computed
  2. Numerical overflow - Very large numbers causing inf/NaN
  3. Missing data - NA values propagating through calculations

Solution:

# Check for problematic values
any(is.na(A))      # Check for NAs
any(is.infinite(A)) # Check for infinite values
det(A)              # If det ≈ 0, matrix is singular

Common Questions

Q: What’s the difference between * and %*%? A: * is element-wise multiplication, %*% is matrix multiplication (dot product of rows and columns)

Q: How do I check if a matrix is invertible? A: Check det(matrix) != 0 or verify all eigenvalues are non-zero

Q: What’s the best way to solve Ax = b? A: Use solve(A, b) - it’s optimized and handles numerical issues

Q: When should I use SVD vs inverse? A: SVD is more numerically stable and works for non-square matrices; inverse is simpler for square invertible matrices

Q: How do I compute matrix powers? A: Use %*% repeatedly: A %*% A for A², or use eigendecomposition for efficiency

Build on matrix operations for advanced analysis:

Download R Script

Get all code examples from this tutorial: matrices-linear-algebra-examples.R