we can use the rref() function from the pracma package to compute the inverse of an invertible matrix. For example, suppose we have a matrix
library(matlib)
A <- matrix(c(1, 3, 2, 1, 2, 5, -2, -3, 2, 2, -3, -3,
-4, -2, -1, 4), nrow = 4, ncol = 4, byrow=TRUE)Then we create a matrix M with the function diag():
M <- cbind(A, diag(4))Now we upload the pracma package.
library(pracma)##
## Attaching package: 'pracma'## The following objects are masked from 'package:matlib':
##
## angle, invThen call the rref() function.
R <- rref(M)If you want to obtain the inverse of the matrix A, then we take the fifth, sixth, seventh and eighth columns of R.
Ainv <- R[, 5:8]
Ainv## [,1] [,2] [,3] [,4]
## [1,] 0.58024691 -0.7037037 0.8641975 -0.02469136
## [2,] 0.04938272 0.2592593 -0.1604938 0.06172840
## [3,] -0.14814815 0.2222222 -0.5185185 -0.18518519
## [4,] 0.56790123 -0.5185185 0.6543210 0.20987654To check if it is the inverse of the matrix A, we can use the inv() function to obtain the inverse of A:
inv(A)## [,1] [,2] [,3] [,4]
## [1,] 0.58024691 -0.7037037 0.8641975 -0.02469136
## [2,] 0.04938272 0.2592593 -0.1604938 0.06172840
## [3,] -0.14814815 0.2222222 -0.5185185 -0.18518519
## [4,] 0.56790123 -0.5185185 0.6543210 0.20987654====Practical Application===== Now we solve the system of linear equations for the working example by using the rref function to compute the inverse of the coefficient matrix. Recall from the practical application in Section 1.4, we define the coefficient matrix for the original problem
A <- matrix(c(1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ,
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ,
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 0 , 1), nrow = 11, ncol = 18, byrow = TRUE)and we write the right hand side in R
b <- c(1298, 1948, 465, 605, 451, 338, 260, 183, 282, 127, 535)First we create the augmented matrix
C <- cbind(A, b)To eliminate one row from this system, as we did in the previous section, we apply the rref() function from the pracma package. Then we have:
library(pracma)
E <- rref(C)
E## b
## [1,] 1 0 0 0 0 0 0 0 0 0 -1 -1 -1 -1 -1 -1 -1 -1 -1483
## [2,] 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 605
## [3,] 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 451
## [4,] 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 338
## [5,] 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 260
## [6,] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 183
## [7,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 282
## [8,] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 127
## [9,] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 535
## [10,] 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1948
## [11,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0E <- E[-11,]G1 <- eye(8)
G2 <- matrix(rep(0, 80), 8, 10)
b2 <- c(266, 223, 140, 264, 137, 67, 130, 24)
G <- cbind(G1, G2, b2)
M <- rbind(E, G)Note that the eye() function outputs the identity matrix. The command “M[, -19]” creates the coefficient matrix for the system of linear equations
M2 <- M[,-19]Now we want to obtain the inverse of this matrix using the rref() function. Since this matrix is an 18 × 18 matrix, we add the identity matrix of size 18.
M3 <- cbind(M2, diag(18))
M4 <- rref(M3)
M4##
## [1,] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## [2,] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## [3,] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
## [4,] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 -1 -1 -1
## [10,] 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 -1 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 -1 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 -1
## [13,] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## [14,] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## [15,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## [16,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## [17,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## [18,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 -1 -1 -1 -1 -1 -1 -1 0 0 1 1 1
##
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 1 0 0 0 0
## [5,] 0 1 0 0 0
## [6,] 0 0 1 0 0
## [7,] 0 0 0 1 0
## [8,] 0 0 0 0 1
## [9,] -1 -1 -1 -1 -1
## [10,] 0 0 0 0 0
## [11,] 0 0 0 0 0
## [12,] 0 0 0 0 0
## [13,] -1 0 0 0 0
## [14,] 0 -1 0 0 0
## [15,] 0 0 -1 0 0
## [16,] 0 0 0 -1 0
## [17,] 0 0 0 0 -1
## [18,] 1 1 1 1 1Note that the 19th column vector to the 36th column vector form the inverse of the matrix M2. Thus we have
Minv <- M4[, 19:36]
Minv %*% M[,19]## [,1]
## [1,] 266
## [2,] 223
## [3,] 140
## [4,] 264
## [5,] 137
## [6,] 67
## [7,] 130
## [8,] 24
## [9,] 47
## [10,] 199
## [11,] 382
## [12,] 311
## [13,] 74
## [14,] 123
## [15,] 116
## [16,] 152
## [17,] 103
## [18,] 488