Chapter 11 Matrix Operations

What You’ll Learn:

  • Matrix creation and structure
  • Matrix algebra operations
  • Dimension requirements
  • Transpose and multiplication
  • Common matrix errors

Key Errors Covered: 12+ matrix operation errors

Difficulty: ⭐⭐ Intermediate

11.1 Introduction

Matrices are fundamental to many R operations, especially statistics and linear algebra:

A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)

# Try to add them
A + B
#> Error in A + B: non-conformable arrays

🔴 ERROR

Error in A + B : non-conformable arrays

Let’s master matrix operations and avoid dimension mismatches.

11.2 Matrix Basics

💡 Key Insight: Matrices vs Data Frames

# Matrix: all same type
mat <- matrix(1:6, nrow = 2, ncol = 3)
typeof(mat)      # "integer"
#> [1] "integer"
is.matrix(mat)
#> [1] TRUE
is.data.frame(mat)
#> [1] FALSE

# Data frame: can mix types
df <- data.frame(
  x = 1:2,
  y = c("a", "b")
)
is.matrix(df)
#> [1] FALSE
is.data.frame(df)
#> [1] TRUE

# Can convert
as.matrix(df)    # Coerces to character!
#>      x   y  
#> [1,] "1" "a"
#> [2,] "2" "b"
as.data.frame(mat)
#>   V1 V2 V3
#> 1  1  3  5
#> 2  2  4  6

# Matrix properties
dim(mat)         # 2 3 (rows, cols)
#> [1] 2 3
nrow(mat)
#> [1] 2
ncol(mat)
#> [1] 3
length(mat)      # 6 (total elements)
#> [1] 6

Key differences: - Matrices: All same type, 2D array - Data frames: Can mix types, list of vectors

11.3 Error #1: non-conformable arrays

⭐⭐ INTERMEDIATE 📏 DIMENSION

11.3.1 The Error

A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:10, nrow = 2, ncol = 5)

A + B  # Different dimensions!
#> Error in A + B: non-conformable arrays

🔴 ERROR

Error in A + B : non-conformable arrays

11.3.2 What It Means

For element-wise operations (+, -, *, /), matrices must have the same dimensions.

11.3.3 Conformability Rules

# Same dimensions - OK
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(7:12, nrow = 2, ncol = 3)
A + B
#>      [,1] [,2] [,3]
#> [1,]    8   12   16
#> [2,]   10   14   18

# Scalar - OK (recycled)
A + 10
#>      [,1] [,2] [,3]
#> [1,]   11   13   15
#> [2,]   12   14   16

# Vector recycling
A + c(1, 2)  # Recycles down columns
#>      [,1] [,2] [,3]
#> [1,]    2    4    6
#> [2,]    4    6    8

# But these fail:
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)
A + B  # Different dimensions
#> Error in A + B: non-conformable arrays

11.3.4 Common Causes

11.3.4.1 Cause 1: Transposed Matrix

A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)

# B is transpose of A shape
A + B  # Error
#> Error in A + B: non-conformable arrays
# Fix: transpose one
A + t(B)
#>      [,1] [,2] [,3]
#> [1,]    2    5    8
#> [2,]    6    9   12

11.3.4.2 Cause 2: Wrong Construction

# Meant to be same size
A <- matrix(1:6, nrow = 2)  # 2x3
B <- matrix(1:8, nrow = 2)  # 2x4

A + B  # Error
#> Error in A + B: non-conformable arrays

11.3.4.3 Cause 3: After Subsetting

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

# Subset changes dimensions
A_sub <- A[, 1:2]  # Now 3x2
B_sub <- B[, 1:3]  # Now 3x3

A_sub + B_sub  # Error!
#> Error in A_sub + B_sub: non-conformable arrays

11.3.5 Solutions

SOLUTION 1: Check Dimensions First

safe_matrix_add <- function(A, B) {
  if (!identical(dim(A), dim(B))) {
    stop("Matrices have different dimensions: ",
         paste(dim(A), collapse = "x"), " vs ",
         paste(dim(B), collapse = "x"))
  }
  
  return(A + B)
}

# Test
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(7:12, nrow = 2, ncol = 3)
safe_matrix_add(A, B)
#>      [,1] [,2] [,3]
#> [1,]    8   12   16
#> [2,]   10   14   18

SOLUTION 2: Reshape to Match

A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)

# Transpose B to match
B_matched <- t(B)
A + B_matched
#>      [,1] [,2] [,3]
#> [1,]    2    5    8
#> [2,]    6    9   12

# Or reshape
B_reshaped <- matrix(B, nrow = 2, ncol = 3)
A + B_reshaped
#>      [,1] [,2] [,3]
#> [1,]    2    6   10
#> [2,]    4    8   12

SOLUTION 3: Extract Common Dimensions

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

# Find common dimensions
common_rows <- min(nrow(A), nrow(B))
common_cols <- min(ncol(A), ncol(B))

# Extract submatrices
A_sub <- A[1:common_rows, 1:common_cols]
B_sub <- B[1:common_rows, 1:common_cols]

A_sub + B_sub
#>      [,1] [,2] [,3] [,4]
#> [1,]    2    8   14   20
#> [2,]    4   10   16   22
#> [3,]    6   12   18   24

11.4 Error #2: non-conformable arguments

⭐⭐ INTERMEDIATE 📏 DIMENSION

11.4.1 The Error

A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 2, ncol = 3)

# Try matrix multiplication
A %*% B
#> Error in A %*% B: non-conformable arguments

🔴 ERROR

Error in A %*% B : non-conformable arguments

11.4.2 What It Means

For matrix multiplication (%*%), the number of columns in A must equal the number of rows in B.

11.4.3 Matrix Multiplication Rules

💡 Matrix Multiplication Requirements

For A %*% B: - A must be m × n - B must be n × p - Result will be m × p

# A is 2×3, B is 3×2 - OK
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)
dim(A)  # 2 3
#> [1] 2 3
dim(B)  # 3 2
#> [1] 3 2

result <- A %*% B
dim(result)  # 2 2 (outer dimensions)
#> [1] 2 2

Rule: Inner dimensions must match, outer dimensions form result.

(2 × 3) %*% (3 × 2) = (2 × 2)
     ↑       ↑
     └───────┘ must match

11.4.4 Common Causes

11.4.4.1 Cause 1: Wrong Order

A <- matrix(1:6, nrow = 2, ncol = 3)   # 2×3
B <- matrix(1:10, nrow = 2, ncol = 5)  # 2×5

# A has 3 cols, B has 2 rows - mismatch!
A %*% B
#> Error in A %*% B: non-conformable arguments
# But reverse works!
B %*% A  # 2×5 times 5×3... wait, A is 2×3
#> Error in B %*% A: non-conformable arguments

# Need to transpose A
B %*% t(A)  # 2×5 times 5×2 = 2×2
#> Error in B %*% t(A): non-conformable arguments

11.4.4.2 Cause 2: Using Element-wise Instead

A <- matrix(1:4, nrow = 2)
B <- matrix(5:8, nrow = 2)

# Element-wise multiplication (different!)
A * B  # Hadamard product
#>      [,1] [,2]
#> [1,]    5   21
#> [2,]   12   32

# Matrix multiplication
A %*% t(B)  # Need to transpose
#>      [,1] [,2]
#> [1,]   26   30
#> [2,]   38   44

11.4.5 Solutions

SOLUTION 1: Check Conformability

can_multiply <- function(A, B) {
  ncol(A) == nrow(B)
}

safe_matrix_mult <- function(A, B) {
  if (!can_multiply(A, B)) {
    stop("Cannot multiply: A is ", nrow(A), "×", ncol(A),
         ", B is ", nrow(B), "×", ncol(B),
         "\nNeed ncol(A) = nrow(B)")
  }
  
  return(A %*% B)
}

# Test
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)
safe_matrix_mult(A, B)
#>      [,1] [,2]
#> [1,]   22   49
#> [2,]   28   64

SOLUTION 2: Auto-transpose if Needed

smart_mult <- function(A, B) {
  # Try as-is
  if (ncol(A) == nrow(B)) {
    return(A %*% B)
  }
  
  # Try transposing B
  if (ncol(A) == ncol(B)) {
    message("Transposing B")
    return(A %*% t(B))
  }
  
  # Try transposing A
  if (nrow(A) == nrow(B)) {
    message("Transposing A")
    return(t(A) %*% B)
  }
  
  stop("Matrices not conformable in any configuration")
}

# Test
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 2, ncol = 3)
smart_mult(A, B)  # Transposes B
#> Transposing B
#>      [,1] [,2]
#> [1,]   35   44
#> [2,]   44   56

11.5 Error #3: system is computationally singular

⭐⭐⭐ ADVANCED 🧮 MATH

11.5.1 The Error

# Singular matrix (not invertible)
A <- matrix(c(1, 2, 2, 4), nrow = 2)
A
#>      [,1] [,2]
#> [1,]    1    2
#> [2,]    2    4

solve(A)  # Try to invert
#> Error in solve.default(A): Lapack routine dgesv: system is exactly singular: U[2,2] = 0

🔴 ERROR

Error in solve.default(A) : 
  system is computationally singular: reciprocal condition number = 0

11.5.2 What It Means

The matrix is singular (non-invertible). Its determinant is 0 (or very close to 0).

11.5.3 Why Matrices Become Singular

# Linearly dependent rows
A <- matrix(c(1, 2, 2, 4), nrow = 2)
A
#>      [,1] [,2]
#> [1,]    1    2
#> [2,]    2    4
# Row 2 = 2 * Row 1

det(A)  # 0 (singular)
#> [1] 0

# Compare to invertible matrix
B <- matrix(c(1, 2, 3, 4), nrow = 2)
det(B)  # -2 (non-zero, invertible)
#> [1] -2
solve(B)  # Works
#>      [,1] [,2]
#> [1,]   -2  1.5
#> [2,]    1 -0.5

11.5.4 Common Causes

11.5.4.1 Cause 1: Perfect Collinearity

# Data with perfect correlation
x1 <- 1:5
x2 <- 2 * x1  # Perfectly correlated

X <- cbind(1, x1, x2)  # Design matrix
A <- t(X) %*% X  # X'X matrix

solve(A)  # Singular!
#> Error in solve.default(A): Lapack routine dgesv: system is exactly singular: U[3,3] = 0

11.5.4.2 Cause 2: More Variables Than Observations

# 3 observations, 5 variables
X <- matrix(rnorm(15), nrow = 3, ncol = 5)
A <- t(X) %*% X  # 5×5 matrix

solve(A)  # Singular!
#> Error in solve.default(A): Lapack routine dgesv: system is exactly singular: U[5,5] = 0

11.5.4.3 Cause 3: Numerical Issues

# Very small numbers can cause numerical singularity
A <- matrix(c(1, 1e-10, 1e-10, 1), nrow = 2)
solve(A)  # May fail due to numerical precision
#>        [,1]   [,2]
#> [1,]  1e+00 -1e-10
#> [2,] -1e-10  1e+00

11.5.5 Solutions

SOLUTION 1: Check Before Inverting

safe_solve <- function(A, tol = 1e-10) {
  # Check if square
  if (nrow(A) != ncol(A)) {
    stop("Matrix must be square")
  }
  
  # Check determinant
  d <- det(A)
  
  if (abs(d) < tol) {
    stop("Matrix is singular (det = ", d, ")")
  }
  
  return(solve(A))
}

# Test
B <- matrix(c(1, 2, 3, 4), nrow = 2)
safe_solve(B)  # Works
#>      [,1] [,2]
#> [1,]   -2  1.5
#> [2,]    1 -0.5

A <- matrix(c(1, 2, 2, 4), nrow = 2)
safe_solve(A)  # Clear error message
#> Error in safe_solve(A): Matrix is singular (det = 0)

SOLUTION 2: Use Generalized Inverse

library(MASS)

# Singular matrix
A <- matrix(c(1, 2, 2, 4), nrow = 2)

# Moore-Penrose generalized inverse
A_inv <- ginv(A)
A_inv
#>      [,1] [,2]
#> [1,] 0.04 0.08
#> [2,] 0.08 0.16

# Check: A %*% ginv(A) %*% A = A
all.equal(A, A %*% A_inv %*% A)
#> [1] TRUE

SOLUTION 3: Remove Collinear Variables

# Detect and remove collinear columns
remove_collinear <- function(X, threshold = 0.99) {
  cor_matrix <- cor(X)
  
  # Find highly correlated pairs
  high_cor <- which(abs(cor_matrix) > threshold & 
                    upper.tri(cor_matrix, diag = FALSE),
                    arr.ind = TRUE)
  
  if (nrow(high_cor) > 0) {
    # Remove second column of correlated pairs
    remove_cols <- unique(high_cor[, 2])
    message("Removing collinear columns: ", 
            paste(remove_cols, collapse = ", "))
    X <- X[, -remove_cols]
  }
  
  return(X)
}

# Test
x1 <- 1:5
x2 <- 2 * x1
x3 <- rnorm(5)
X <- cbind(x1, x2, x3)

X_clean <- remove_collinear(X)
#> Removing collinear columns: 2
ncol(X_clean)  # One less column
#> [1] 2

11.6 Matrix Creation Errors

⚠️ Common Pitfall: Matrix Filling

# Matrix fills by COLUMN (default)
matrix(1:6, nrow = 2, ncol = 3)
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    2    4    6

# To fill by row:
matrix(1:6, nrow = 2, ncol = 3, byrow = TRUE)
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    4    5    6

# This catches many people!
matrix(c(1, 2, 3,
         4, 5, 6), nrow = 2, ncol = 3)
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    2    4    6
# NOT what you might expect!

# Want row-wise? Use byrow:
matrix(c(1, 2, 3,
         4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE)
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    4    5    6

11.7 Matrix Operations Reference

🎯 Best Practice: Common Matrix Operations

A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)

# Transpose
t(A)
#>      [,1] [,2]
#> [1,]    1    2
#> [2,]    3    4
#> [3,]    5    6

# Matrix multiplication
A %*% B  # Result: 2×2
#>      [,1] [,2]
#> [1,]   22   49
#> [2,]   28   64

# Element-wise operations (same dimensions needed)
C <- matrix(7:12, nrow = 2, ncol = 3)
A + C
#>      [,1] [,2] [,3]
#> [1,]    8   12   16
#> [2,]   10   14   18
A - C
#>      [,1] [,2] [,3]
#> [1,]   -6   -6   -6
#> [2,]   -6   -6   -6
A * C  # Hadamard product (element-wise)
#>      [,1] [,2] [,3]
#> [1,]    7   27   55
#> [2,]   16   40   72
A / C
#>           [,1]      [,2]      [,3]
#> [1,] 0.1428571 0.3333333 0.4545455
#> [2,] 0.2500000 0.4000000 0.5000000

# Cross product
crossprod(A)    # t(A) %*% A
#>      [,1] [,2] [,3]
#> [1,]    5   11   17
#> [2,]   11   25   39
#> [3,]   17   39   61
tcrossprod(A)   # A %*% t(A)
#>      [,1] [,2]
#> [1,]   35   44
#> [2,]   44   56

# Determinant
D <- matrix(c(1, 2, 3, 4), nrow = 2)
det(D)
#> [1] -2

# Inverse (square matrices only)
solve(D)
#>      [,1] [,2]
#> [1,]   -2  1.5
#> [2,]    1 -0.5

# Diagonal
diag(D)          # Extract diagonal
#> [1] 1 4
diag(c(1, 2, 3)) # Create diagonal matrix
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    2    0
#> [3,]    0    0    3

# Eigenvalues and eigenvectors
eigen(D)
#> eigen() decomposition
#> $values
#> [1]  5.3722813 -0.3722813
#> 
#> $vectors
#>            [,1]       [,2]
#> [1,] -0.5657675 -0.9093767
#> [2,] -0.8245648  0.4159736

# Singular value decomposition
svd(A)
#> $d
#> [1] 9.5255181 0.5143006
#> 
#> $u
#>            [,1]       [,2]
#> [1,] -0.6196295 -0.7848945
#> [2,] -0.7848945  0.6196295
#> 
#> $v
#>            [,1]       [,2]
#> [1,] -0.2298477  0.8834610
#> [2,] -0.5247448  0.2407825
#> [3,] -0.8196419 -0.4018960

11.8 Dimension Preservation

⚠️ Common Pitfall: Dropping Dimensions

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

# Extract row (becomes vector!)
row1 <- A[1, ]
dim(row1)  # NULL (it's a vector now)
#> NULL

# Extract column (becomes vector!)
col1 <- A[, 1]
dim(col1)  # NULL
#> NULL

# Preserve matrix structure
row1 <- A[1, , drop = FALSE]
dim(row1)  # 1 3
#> [1] 1 4

col1 <- A[, 1, drop = FALSE]
dim(col1)  # 3 1
#> [1] 3 1

When it matters:

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

# Extract column from A (becomes vector)
A_col <- A[, 1]

# Try to multiply
A_col %*% B  # Error! A_col is vector
#>      [,1]
#> [1,]   14
# Fix: preserve dimensions
A_col <- A[, 1, drop = FALSE]
A_col %*% t(B)  # Works!
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    2    4    6
#> [3,]    3    6    9

11.9 Converting Between Structures

💡 Key Insight: Conversions

# Vector to matrix
vec <- 1:12
mat <- matrix(vec, nrow = 3, ncol = 4)

# Matrix to vector
as.vector(mat)  # Column-major order
#>  [1]  1  2  3  4  5  6  7  8  9 10 11 12

# Matrix to data frame
df <- as.data.frame(mat)
class(df)
#> [1] "data.frame"

# Data frame to matrix
mat2 <- as.matrix(df)
class(mat2)
#> [1] "matrix" "array"

# List to matrix (if all same length)
lst <- list(a = 1:3, b = 4:6, c = 7:9)
mat3 <- do.call(cbind, lst)
mat3
#>      a b c
#> [1,] 1 4 7
#> [2,] 2 5 8
#> [3,] 3 6 9

# Matrix to list (by column)
lst2 <- as.list(as.data.frame(mat))

Warning: Type coercion

# Data frame with mixed types
df_mixed <- data.frame(
  x = 1:3,
  y = c("a", "b", "c")
)

# Converting to matrix coerces to common type
as.matrix(df_mixed)  # All become character!
#>      x   y  
#> [1,] "1" "a"
#> [2,] "2" "b"
#> [3,] "3" "c"

11.10 Summary

Key Takeaways:

  1. Element-wise operations: Need identical dimensions
  2. Matrix multiplication: Inner dimensions must match
  3. Singular matrices: Cannot be inverted (det = 0)
  4. Filling order: Column-major by default (use byrow = TRUE)
  5. drop = FALSE: Preserves matrix structure
  6. Type coercion: Converting mixed-type df to matrix coerces all
  7. Check dimensions: Always verify before operations

Quick Reference:

Error Cause Fix
non-conformable arrays Different dimensions for +,-,*,/ Match dimensions
non-conformable arguments ncol(A) ≠ nrow(B) for %*% Transpose or reshape
computationally singular Matrix not invertible Check det(), use ginv()
incorrect number of dimensions Wrong subscripts Match matrix structure

Matrix Operations Checklist:

# Before operations:
dim(A)                    # Check dimensions
det(A)                    # Check if invertible
ncol(A) == nrow(B)        # Check for multiplication

# Safe extraction:
A[i, , drop = FALSE]      # Preserve row
A[, j, drop = FALSE]      # Preserve column

# Matrix multiplication:
A %*% B                   # Matrix product
A * B                     # Element-wise (Hadamard)
crossprod(A, B)           # t(A) %*% B
tcrossprod(A, B)          # A %*% t(B)

Best Practices:

# ✅ Good
identical(dim(A), dim(B))      # Check before adding
ncol(A) == nrow(B)             # Check before multiply
det(A) != 0                    # Check before invert
A[i, , drop = FALSE]           # Preserve structure

# ❌ Avoid
A + B                          # Without checking dims
solve(A)                       # Without checking singular
A[i, ]                         # Drops to vector unexpectedly

11.11 Exercises

📝 Exercise 1: Matrix Dimension Checker

Write a function that checks if two matrices can be: 1. Added/subtracted 2. Multiplied (A %% B) 3. Multiplied (B %% A)

Return TRUE/FALSE for each operation.

📝 Exercise 2: Safe Matrix Operations

Create matrix_op(A, B, op) that: - Checks dimensions before operation - Supports: “add”, “subtract”, “multiply”, “divide” - Gives clear error messages - Returns result or NULL

📝 Exercise 3: Matrix Inversion Check

Write safe_invert(A) that: 1. Checks if matrix is square 2. Checks if singular 3. Warns if near-singular 4. Returns inverse or NULL 5. Provides diagnostic information

📝 Exercise 4: Matrix Creation Helper

Write make_matrix(...) that: - Takes values and shape (nrow, ncol) - Handles different input formats (vector, list, data frame) - Validates dimensions - Allows byrow specification - Returns matrix with informative errors

11.12 Exercise Answers

Click to see answers

Exercise 1:

check_matrix_ops <- function(A, B) {
  result <- list(
    can_add = identical(dim(A), dim(B)),
    can_multiply_AB = ncol(A) == nrow(B),
    can_multiply_BA = ncol(B) == nrow(A)
  )
  
  # Add details
  result$dim_A <- paste(dim(A), collapse = "×")
  result$dim_B <- paste(dim(B), collapse = "×")
  
  if (result$can_multiply_AB) {
    result$result_dim_AB <- paste(c(nrow(A), ncol(B)), collapse = "×")
  }
  
  if (result$can_multiply_BA) {
    result$result_dim_BA <- paste(c(nrow(B), ncol(A)), collapse = "×")
  }
  
  class(result) <- "matrix_ops_check"
  return(result)
}

print.matrix_ops_check <- function(x, ...) {
  cat("Matrix A:", x$dim_A, "\n")
  cat("Matrix B:", x$dim_B, "\n\n")
  cat("Can add/subtract:", x$can_add, "\n")
  cat("Can multiply A %*% B:", x$can_multiply_AB)
  if (x$can_multiply_AB) {
    cat(" (result:", x$result_dim_AB, ")")
  }
  cat("\n")
  cat("Can multiply B %*% A:", x$can_multiply_BA)
  if (x$can_multiply_BA) {
    cat(" (result:", x$result_dim_BA, ")")
  }
  cat("\n")
}

# Test
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)
check_matrix_ops(A, B)
#> Matrix A: 2×3 
#> Matrix B: 3×2 
#> 
#> Can add/subtract: FALSE 
#> Can multiply A %*% B: TRUE (result: 2×2 )
#> Can multiply B %*% A: TRUE (result: 3×3 )

Exercise 2:

matrix_op <- function(A, B, op = c("add", "subtract", "multiply", "divide")) {
  op <- match.arg(op)
  
  # Validate inputs
  if (!is.matrix(A) || !is.matrix(B)) {
    stop("Both A and B must be matrices")
  }
  
  # Check dimensions based on operation
  if (op %in% c("add", "subtract", "divide")) {
    if (!identical(dim(A), dim(B))) {
      stop("For ", op, ", matrices must have same dimensions. ",
           "A is ", paste(dim(A), collapse = "×"),
           ", B is ", paste(dim(B), collapse = "×"))
    }
    
    result <- switch(op,
      add = A + B,
      subtract = A - B,
      divide = A / B
    )
  } else if (op == "multiply") {
    if (ncol(A) != nrow(B)) {
      stop("For multiplication, ncol(A) must equal nrow(B). ",
           "A is ", paste(dim(A), collapse = "×"),
           ", B is ", paste(dim(B), collapse = "×"))
    }
    
    result <- A %*% B
  }
  
  return(result)
}

# Test
A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(1:6, nrow = 3, ncol = 2)
matrix_op(A, B, "multiply")
#>      [,1] [,2]
#> [1,]   22   49
#> [2,]   28   64

Exercise 3:

safe_invert <- function(A, tol = 1e-10, warn_threshold = 1e-8) {
  # Check if matrix
  if (!is.matrix(A)) {
    message("Input is not a matrix")
    return(NULL)
  }
  
  # Check if square
  if (nrow(A) != ncol(A)) {
    message("Matrix is not square: ", 
            paste(dim(A), collapse = "×"))
    return(NULL)
  }
  
  # Calculate determinant
  d <- det(A)
  
  # Check if singular
  if (abs(d) < tol) {
    message("Matrix is singular (det = ", d, ")")
    message("Consider using MASS::ginv() for generalized inverse")
    return(NULL)
  }
  
  # Warn if near-singular
  if (abs(d) < warn_threshold) {
    warning("Matrix is near-singular (det = ", d, "). ",
            "Results may be numerically unstable.")
  }
  
  # Compute inverse
  A_inv <- solve(A)
  
  # Verify (optional)
  check <- A %*% A_inv
  is_identity <- all(abs(check - diag(nrow(A))) < 1e-10)
  
  if (!is_identity) {
    warning("Inversion may be inaccurate (A %*% A^-1 != I)")
  }
  
  # Return with diagnostics
  attr(A_inv, "determinant") <- d
  attr(A_inv, "condition_number") <- kappa(A)
  
  return(A_inv)
}

# Test
A <- matrix(c(1, 2, 3, 4), nrow = 2)
A_inv <- safe_invert(A)
A_inv
#>      [,1] [,2]
#> [1,]   -2  1.5
#> [2,]    1 -0.5
#> attr(,"determinant")
#> [1] -2
#> attr(,"condition_number")
#> [1] 18.77778

# Singular
B <- matrix(c(1, 2, 2, 4), nrow = 2)
safe_invert(B)
#> Matrix is singular (det = 0)
#> Consider using MASS::ginv() for generalized inverse
#> NULL

Exercise 4:

make_matrix <- function(x, nrow, ncol, byrow = FALSE) {
  # Handle different input types
  if (is.matrix(x)) {
    if (missing(nrow) && missing(ncol)) {
      return(x)
    }
    x <- as.vector(x)
  } else if (is.data.frame(x)) {
    x <- as.matrix(x)
    if (missing(nrow) && missing(ncol)) {
      return(x)
    }
    x <- as.vector(x)
  } else if (is.list(x)) {
    # Check if all elements same length
    lens <- lengths(x)
    if (length(unique(lens)) != 1) {
      stop("List elements have different lengths")
    }
    x <- unlist(x)
  }
  
  # Validate dimensions
  n_elements <- length(x)
  
  if (missing(nrow) && missing(ncol)) {
    stop("Must provide nrow, ncol, or both")
  }
  
  if (missing(nrow)) {
    nrow <- ceiling(n_elements / ncol)
  } else if (missing(ncol)) {
    ncol <- ceiling(n_elements / nrow)
  }
  
  expected_elements <- nrow * ncol
  
  if (n_elements != expected_elements) {
    if (n_elements < expected_elements) {
      warning("Data length (", n_elements, 
              ") is less than matrix size (", expected_elements,
              "). Recycling values.")
    } else {
      warning("Data length (", n_elements,
              ") is greater than matrix size (", expected_elements,
              "). Truncating values.")
      x <- x[1:expected_elements]
    }
  }
  
  # Create matrix
  result <- matrix(x, nrow = nrow, ncol = ncol, byrow = byrow)
  
  return(result)
}

# Test
make_matrix(1:6, nrow = 2, ncol = 3)
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    2    4    6
make_matrix(1:6, nrow = 2, byrow = TRUE)
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    4    5    6
make_matrix(list(a = 1:3, b = 4:6), ncol = 3)
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    2    4    6