Gaussian elimination figures prominently in the matrixDemos package, and is implemented directly in the GaussianElimination
function. Gaussian elimination is also used by matrixInverse
to compute the inverse of a nonsingular square matrix; by Ginv
to compute a generalized inverse of a singular square matrix or a rectangular matrix; and by RREF
to compute the reduced row-echelon form of a matrix.
GaussianElimination
has one required argument, the matrix A
, which is transformed to reduced row-echelon form; for example:
> (A <- matrix(c(1,2,3,4,5,6,7,8,10), 3, 3)) # a nonsingular matrix
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 10
> GaussianElimination(A)
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> (B <- matrix(1:9, 3, 3)) # a singular matrix
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> GaussianElimination(B)
[,1] [,2] [,3]
[1,] 1 0 -1
[2,] 0 1 2
[3,] 0 0 0
In these applications, GaussianElimination
is equivalent to RREF
.
A second optional argument to GaussianElimination
is B
, which may be a matrix or a vector, and which is transformed in the same manner as A
. For example,
> A # nonsingular
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 10
> GaussianElimination(A, diag(3)) # find inverse of A
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 -0.6666667 -0.6666667 1
[2,] 0 1 0 -1.3333333 3.6666667 -2
[3,] 0 0 1 1.0000000 -2.0000000 1
> b <- 1:3
> GaussianElimination(A, b) # solve the matrix equation Ax = b
[,1] [,2] [,3] [,4]
[1,] 1 0 0 1
[2,] 0 1 0 0
[3,] 0 0 1 0
For the last example, the solution is \(x_1 = 1, x_2 = x_3 = 0\).
Other arguments to GaussianElimination
:
tol
, allows you to control the tolerance for checking that a number is (effectively) 0 within rounding error
verbose
, if TRUE
(default is FALSE
) shows the process of reducing A
to RREF
fractions
, if TRUE
(default is FALSE
) tries to express the result as rational (rather than decimal) fractions
Some further examples:
> matrixInverse(A, fractions=TRUE)
[,1] [,2] [,3]
[1,] -2/3 -2/3 1
[2,] -4/3 11/3 -2
[3,] 1 -2 1
> RREF(A, verbose=TRUE)
[,1] [,2] [,3]
[1,] 1 2 3.333333
[2,] 0 1 1.333333
[3,] 0 2 3.666667
[,1] [,2] [,3]
[1,] 1 0 -0.3333333
[2,] 0 1 1.8333333
[3,] 0 0 -0.5000000
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> Ginv(B, fractions=TRUE) # a generalized inverse of singular B
[,1] [,2] [,3]
[1,] -3/4 0 7/12
[2,] 0 0 0
[3,] 1/4 0 -1/12
> B %*% Ginv(B) %*% B # check
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> (C <- matrix(1:12, 3, 4)) # rectangular matrix
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> Ginv(C, fractions=TRUE)
[,1] [,2] [,3]
[1,] -2/3 0 5/9
[2,] 0 0 0
[3,] 0 0 0
[4,] 1/6 0 -1/18
> round(C %*% Ginv(C) %*% C, 5) # check
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12