4. Matrix inversion by elementary row operations

Michael Friendly

2024-10-05

The following examples illustrate the steps in finding the inverse of a matrix using elementary row operations (EROs):

• Add a multiple of one row to another (rowadd())
• Multiply one row by a constant (rowmult())
• Interchange two rows (rowswap())

These have the properties that they do not change the inverse. The method used here is sometimes called the Gauss-Jordan method, a form of Gaussian elimination. Another term is (row-reduced) echelon form.

Steps:

1. Adjoin the identity matrix to the right side of A, to give the matrix $[A | I]$
2. Apply row operations to this matrix until the left ($A$) side is reduced to $I$
3. The inverse matrix appears in the right ($I$) side

Why this works: The series of row operations transforms $$[A | I] \Rightarrow [A^{-1} A | A^{-1} I] = [I | A^{-1}]$$

If the matrix is does not have an inverse (is singular) a row of all zeros will appear in the left ($A$) side.

Load the matlib package

library(matlib)

Create a 3 x 3 matrix

   A <- matrix( c(1, 2, 3,
                  2, 3, 0,
                  0, 1,-2), nrow=3, byrow=TRUE)

Join an identity matrix to A

   (AI <-  cbind(A, diag(3)))

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    2    3    1    0    0
## [2,]    2    3    0    0    1    0
## [3,]    0    1   -2    0    0    1

Apply elementary row operations to reduce A to an identity matrix.

The right three cols will then contain inv(A). We will do this three ways:

1. first, just using R arithmetic on the rows of AI
2. using the ERO functions in the matlib package
3. using the echelon() function

1. Using R arithmetic

    (AI[2,] <- AI[2,] - 2*AI[1,])     # row 2 <- row 2 - 2 * row 1

## [1]  0 -1 -6 -2  1  0

    (AI[3,] <- AI[3,] + AI[2,])       # row 3 <- row 3 + row 2

## [1]  0  0 -8 -2  1  1

    (AI[2,] <- -1 * AI[2,])           # row 2 <- -1 * row 2

## [1]  0  1  6  2 -1  0

    (AI[3,] <-  -(1/8) * AI[3,])        # row 3 <- -.25 * row 3

## [1]  0.000  0.000  1.000  0.250 -0.125 -0.125

Now, all elements below the diagonal are zero

    AI

##      [,1] [,2] [,3] [,4]   [,5]   [,6]
## [1,]    1    2    3 1.00  0.000  0.000
## [2,]    0    1    6 2.00 -1.000  0.000
## [3,]    0    0    1 0.25 -0.125 -0.125

      #--continue, making above diagonal == 0
    AI[2,] <- AI[2,] - 6 * AI[3,]     # row 2 <- row 2 - 6 * row 3
    AI[1,] <- AI[1,] - 3 * AI[3,]     # row 1 <- row 1 - 3 * row 3
    AI[1,] <- AI[1,] - 2 * AI[2,]     # row 1 <- row 1 - 2 * row 2

    AI

##      [,1] [,2] [,3]  [,4]   [,5]   [,6]
## [1,]    1    0    0 -0.75  0.875 -1.125
## [2,]    0    1    0  0.50 -0.250  0.750
## [3,]    0    0    1  0.25 -0.125 -0.125

   #-- last three cols are the inverse
  (AInv <- AI[,-(1:3)])

##       [,1]   [,2]   [,3]
## [1,] -0.75  0.875 -1.125
## [2,]  0.50 -0.250  0.750
## [3,]  0.25 -0.125 -0.125

   #-- compare with inv()
  inv(A)

##       [,1]   [,2]   [,3]
## [1,] -0.75  0.875 -1.125
## [2,]  0.50 -0.250  0.750
## [3,]  0.25 -0.125 -0.125

2. Do the same, using matlib functions rowadd(), rowmult() and rowswap()

   AI <-  cbind(A, diag(3))

   AI <- rowadd(AI, 1, 2, -2)        # row 2 <- row 2 - 2 * row 1
   AI <- rowadd(AI, 2, 3, 1)         # row 3 <- row 3 + row 2
   AI <- rowmult(AI, 2, -1)          # row 1 <- -1 * row 2
   AI <- rowmult(AI, 3, -1/8)        # row 3 <- -.25 * row 3

   # show result so far
   AI

##      [,1] [,2] [,3] [,4]   [,5]   [,6]
## [1,]    1    2    3 1.00  0.000  0.000
## [2,]    0    1    6 2.00 -1.000  0.000
## [3,]    0    0    1 0.25 -0.125 -0.125

    #--continue, making above-diagonal == 0
   AI <- rowadd(AI, 3, 2, -6)        # row 2 <- row 2 - 6 * row 3
   AI <- rowadd(AI, 2, 1, -2)        # row 1 <- row 1 - 2 * row 2
   AI <- rowadd(AI, 3, 1, -3)        # row 1 <- row 1 - 3 * row 3
   AI

##      [,1] [,2] [,3]  [,4]   [,5]   [,6]
## [1,]    1    0    0 -0.75  0.875 -1.125
## [2,]    0    1    0  0.50 -0.250  0.750
## [3,]    0    0    1  0.25 -0.125 -0.125

3. Using echelon()

echelon() does all these steps row by row, and returns the result

   echelon( cbind(A, diag(3)))

##      [,1] [,2] [,3]  [,4]   [,5]   [,6]
## [1,]    1    0    0 -0.75  0.875 -1.125
## [2,]    0    1    0  0.50 -0.250  0.750
## [3,]    0    0    1  0.25 -0.125 -0.125

It is more interesting to see the steps, using the argument verbose=TRUE. In many cases, it is informative to see the numbers printed as fractions.

   echelon( cbind(A, diag(3)), verbose=TRUE, fractions=TRUE)

## 
## Initial matrix:

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  1    2    3    1    0    0  
## [2,]  2    3    0    0    1    0  
## [3,]  0    1   -2    0    0    1  
## 
## row: 1 
## 
##  exchange rows 1 and 2

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  2    3    0    0    1    0  
## [2,]  1    2    3    1    0    0  
## [3,]  0    1   -2    0    0    1  
## 
##  multiply row 1 by 1/2

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]   1  3/2    0    0  1/2    0 
## [2,]   1    2    3    1    0    0 
## [3,]   0    1   -2    0    0    1 
## 
##  subtract row 1 from row 2

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1  3/2    0    0  1/2    0
## [2,]    0  1/2    3    1 -1/2    0
## [3,]    0    1   -2    0    0    1
## 
## row: 2 
## 
##  exchange rows 2 and 3

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1  3/2    0    0  1/2    0
## [2,]    0    1   -2    0    0    1
## [3,]    0  1/2    3    1 -1/2    0
## 
##  multiply row 2 by 3/2 and subtract from row 1

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    3    0  1/2 -3/2
## [2,]    0    1   -2    0    0    1
## [3,]    0  1/2    3    1 -1/2    0
## 
##  multiply row 2 by 1/2 and subtract from row 3

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    3    0  1/2 -3/2
## [2,]    0    1   -2    0    0    1
## [3,]    0    0    4    1 -1/2 -1/2
## 
## row: 3 
## 
##  multiply row 3 by 1/4

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    3    0  1/2 -3/2
## [2,]    0    1   -2    0    0    1
## [3,]    0    0    1  1/4 -1/8 -1/8
## 
##  multiply row 3 by 3 and subtract from row 1

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    0 -3/4  7/8 -9/8
## [2,]    0    1   -2    0    0    1
## [3,]    0    0    1  1/4 -1/8 -1/8
## 
##  multiply row 3 by 2 and add to row 2

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    0 -3/4  7/8 -9/8
## [2,]    0    1    0  1/2 -1/4  3/4
## [3,]    0    0    1  1/4 -1/8 -1/8