Problem Set 2

Using the procedure outlined in section 1 of the weekly handout, write a function to compute the inverse of a well-conditioned full-rank square matrix using co-factors. In order to compute the co-factors, you may use built-in commands to compute the determinant. Your function should have the following signature:

\[B = myinverse(A)\]

where A is a matrix and B is its inverse and AxB = I. The off-diagonal elements of I should be close to zero, if not zero. Likewise, the diagonal elements should be close to 1, if not 1. Small numerical precision errors are acceptable but the function myinverse should be correct and must use co-factors and determinant of A to compute the inverse.

Solution

First we’ll define the function myInverse. The approach will be to cycle through our matrix and determine each element of the cofactor matrix, C. We will then compute \(A^{-1} = \frac{C^T}{det(A)}\).

myInverse <- function(A){
  
  n <- nrow(A) # define the dimensions of the matrix
  dimList <- seq(nrow(A)) # create a sequence for the for loop to iterate over
  tempList <- matrix(, nrow = n, ncol = n) # store cofactors here
  
  # for all rows
  for(i in dimList){
  
    # set row indices
    if (i==1){
      rows <- seq((i+1),n)
    } else if (i==n) {
      rows <- seq(1,(i-1))
    } else {
      rowa <- seq(1,(i-1))
      rowb <- seq((i+1),n)
      rows <- c(rowa,rowb)
    }
    
    # for all columns
    for(j in dimList){
  
      # set col indices
      if (j==1){
        cols <- seq((j+1),n)
      } else if (j==n) {
        cols <- seq(1,(j-1))
      } else {
        cola <- seq(1,(j-1))
        colb <- seq((j+1),n)
        cols <- c(cola,colb)
      }
      
      if(n==2){
        finalVal <- ((-1)^(i+j)) * A[rows,cols]
      } else {
        cofacDet <- det(A[rows,cols])
        finalVal <- (-1^(i+j))* cofacDet
      }
      
      # append the cofactor to the appropriate element of the tempList matrix
      tempList[i,j] <- finalVal
    }
  }
  
  Ainv <- t(tempList)/det(A)
  
  return(Ainv)
  
}

We will test our function using the matrix: \[ A = \left[\begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array}\right] \]

First, let’s generate the inverse function:

A <- matrix(c(1,3,2,4),nrow=2)
myInv <- myInverse(A)
myInv
##      [,1] [,2]
## [1,] -2.0  1.0
## [2,]  1.5 -0.5

Next, we will verify that \(AA^{-1}=I\):

A%*%myInv
##              [,1] [,2]
## [1,] 1.000000e+00    0
## [2,] 8.881784e-16    1