Sample R function demonstrating how a 2x2 or 3x3 matrix inverse is computed.

matrix_inverse.R

# Sample function demonstrating how a 2x2 or 3x3 matrix inverse is computed. 
# Requires an object with no more than 3 columns and an equal number of rows that can be coerced into a matrix.
# Used in post on computing the inverse of a matrix if one exists: http://wp.me/p4aZEo-5Ln


matrix.inverse <- function(mat) {
  A <- as.matrix(mat)
  
  # If there are more than four columns in the supplied matrix, stop routine
  if ((ncol(A) >= 4) | (nrow(A) >= 4)) {
    stop('Matrix is not 2x2 or 3x3')
  }
  
  # Stop if matrix is a single column vector
  if (ncol(A) == 1) {
    stop('Matrix is a vector')
  }
  
  # 2x2 matrix inverse
  if(ncol(A) == 2) {
    # Determinant
    a <- A[1]
    b <- A[3]
    c <- A[2]
    d <- A[4]
    det <- a * d - b * c
    # Check to see if matrix is singular
    if (det == 0) {
      stop('Determinant of matrix equals 0, no inverse exists')
    }
    # Compute matrix inverse elements
    a.inv <- d / det
    b.inv <- -b / det
    c.inv <- -c / det
    d.inv <- a / det
    # Collect the results into a new matrix
    inv.mat <- as.matrix(cbind(c(a.inv,c.inv), c(b.inv,d.inv)))
  }
  
  # 3x3 matrix inverse
  if (ncol(A) == 3) {
    # Extract the entries from the matrix
    a <- A[1]
    b <- A[4]
    c <- A[7]
    d <- A[2]
    e <- A[5]
    f <- A[8]
    g <- A[3]
    h <- A[6]
    k <- A[9]
    
    # Compute the determinant and check that it is not 0
    det <- a * (e * k - f * h) - b * (d * k - f * g) + c * (d * h - e * g)
    if (det == 0) {
      stop('Determinant of matrix equals 0, no inverse exists')
    }
    
    # Using the equations defined above, calculate the matrix inverse entries.
    A.inv <- (e * k - f * h) / det
    B.inv <- -(b * k - c * h) / det
    C.inv <- (b * f - c * e) / det
    D.inv <- -(d * k - f * g) / det
    E.inv <- (a * k - c * g) / det
    F.inv <- -(a * f - c * d) / det
    G.inv <- (d * h - e * g) / det
    H.inv <- -(a * h - b * g) / det
    K.inv <- (a * e - b * d) / det
    
    # Collect the results into a new matrix
    inv.mat <- as.matrix(cbind(c(A.inv,D.inv,G.inv), c(B.inv,E.inv,H.inv), c(C.inv,F.inv,K.inv)))
  }
  
  return(inv.mat)
}