Doing mathematics with R: Linear algebra in R
Row reduction and determinant
The R package pracma provides row and column operations to compute reduction processes in matrices numerically, but you can calculate the reduced row echelon form of a matrix with the command rref. This allows you to perform the method of row reduction for calculating an inverse matrix.
Computing the inverse of a matrix
> A <- matrix(c(3, 1, -4, -1, 0, 2, 1, 1, -1), nrow=3, byrow=TRUE)
> library(pracma)
> AI <- cbind(A,diag(3)); AI # augmented matrix
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3 1 -4 1 0 0
[2,] -1 0 2 0 1 0
[3,] 1 1 -1 0 0 1
> rref(AI) # reduced row echelon form
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 2 3 -2
[2,] 0 1 0 -1 -1 2
[3,] 0 0 1 1 2 -1
> .Last.value[,4:6] # submatrix that is equal to inverse of A
[,1] [,2] [,3]
[1,] 2 3 -2
[2,] -1 -1 2
[3,] 1 2 -1
> inv(A) # inverse of a matrix via inv function
[,1] [,2] [,3]
[1,] 2 3 -2
[2,] -1 -1 2
[3,] 1 2 -1
The determinant of the matrix \(A\) in the above example can be calculated with the command det, and the result shows that the matrix is invertible. We give below also an example of a non-invertible matrix because the determinant is equal to 0 and the reduced echelon form is not equal to an identity matrix. If you're still trying to calculate the inverse, you actually get a result with a matching warning. In addition, be always on guard for numerical rounding errors in calculations.
Determinant of a matrix
> A <- matrix(c(3, 1, -4, -1, 0, 2, 1, 1, -1), nrow=3, byrow=TRUE)
> A # invertible matrix
[,1] [,2] [,3]
[1,] 3 1 -4
[2,] -1 0 2
[3,] 1 1 -1
> det(A)
[1] -1
> B <- matrix(c(1, -2, -2, 1, 14, 14, 4, -24, -24), nrow=3, byrow=TRUE)
> B # non-invertible matrix (singular matrix)
[,1] [,2] [,3]
[1,] 1 -2 -2
[2,] 1 14 14
[3,] 4 -24 -24
> det(B) # determinant of a matrix
[1] 0
> qr(B)$rank # rank of matrix B (via QR decomposition)
[1] 2
> library(pracma)
> Rank(B) # rank of matrix B (via pracma package)
> rref(B) # reduced row echelon form of matrix B
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 1
[3,] 0 0 0
> inv(B) # attempt to compute inverse of B
[,1] [,2] [,3]
[1,] Inf Inf Inf
[2,] Inf Inf Inf
[3,] Inf Inf Inf
Warning message:
In inv(B) : Matrix appears to be singular.
