Matrix Inversion by Gaussian Elimination

Description

Matrix Inversion by Gaussian Elimination

Usage

invGE(A, loc = 1:nrow(A), nopermute = 1, eps = 1e-09, simplify = 1, print = 0)

Arguments

Details

This function performs the Gaussian elimination on A matrix. When loc=1:nrow(A), the elimination starts from column 1 to column n=ncol(A) as:

  B0 <- cbind( A, I )
  B1 <- T1 %*% B0
  B2 <- T2 %*% B1
   :        :
  Bk <- Tk %*% Bk-1
   :        :
  Bn <- Tn %*% Bn-1
 

where Tk matrix makes the k-th column of Bk matrix a unit vector whose k-th element is 1 and 0 elsewhere. B[i,i] is used as the pivot.

At the end of the iteration, we have

   Bn = cbind( I, inv(A) )
 

In other words

  inv(A) = Tn %*% Tn-1 %*% ... %*% T2 %*% T1
 

When loc != 1:nrow(A), at each step, B[loc[i],loc[i]] is used as the pivot.

When length(loc) < nrow(A), the intermediate result will be returned.

Value

cbind(I, inv(A)
The inverse of A is given as cbind(I,inv{A})[,-c(1:nrow(A))].

Examples

# Symbolic matrix A.

n <- 2
Aa <- demomat(n,n,root="a")

set.seed(1701)
An <- matrix(rnorm(n*n),n,n)

res <- invGE(Aa, print=1)
II <- res[,c(1:n)]
invAa <- res[,-c(1:n)]

# numerical checking
Eval(II, Aa=An)
# Eval(Aa%*%invAa, Aa=An)
Eval(matTimes(Aa,invAa), Aa=An)


# Numeric matrix A.

set.seed(1702)
n <- 4
A <- matrix( sample(0:9,n*n,replace=1), n,n )

# result
res <- invGE(A, print=1)
invA <- res[,-(1:n)]
printm(A, invA, A%*%invA, fmt=".2")



# Example of permutation
# needs permutation
set.seed(1701)
n=4
A=matrix( sample(0:9,n*n,replace=1), n,n )

# wrong result
invA=invGE(A, print=1)[,c(1:n)]
printm( A%*%invA )

# correct result by nopermute=0 option
invA=invGE(A, print=1, nopermute=0)[,-c(1:n)]
printm( A%*%invA, fmt="7.3"  )

# correct by specifying loc argument
invA=invGE(A, print=1, loc=c(1,4,3,2))[,-c(1:n)]
printm( A%*%invA, fmt="6.3" )



# non full rank matrix A.
set.seed(1701)
n <- 4
A <- matrix( sample(0:9,n*n,replace=1), n,n )
A[,4] <- A[,1]+A[,2]

res <- invGE(A, print=1, nopermute=0)
invA <- res[,-c(1:n)]
round( A%*%invA ,6 )
chkginv(A,invA)