Cramer’s Rule in R

Cramer’s Rule

In linear algebra, Cramer’s rule is an explicit formula for the solution of a system of linear equations with as many equations as unknowns, valid whenever the system has a unique solution. It expresses the solution in terms of the determinants of the (square) coefficient matrix and of matrices obtained from it by replacing one column by the column vector of right-hand-sides of the equations. It is named after Gabriel Cramer (1704–1752).

Cramer’s rule can also be numerically unstable even for 2×2 systems. However, it has recently been shown that Cramer’s rule can be implemented in \(O(n^3)\) time, which is comparable to more common methods of solving systems of linear equations, such as Gaussian elimination (consistently requiring 2.5 times as many arithmetic operations for all matrix sizes), while exhibiting comparable numeric stability in most cases.

Calculation

Consider a system of n linear equations for n unknowns, represented in matrix multiplication form as follows:

\(Ax=B\)

where the n × n matrix \(A\) has a nonzero determinant, and the vector \(x=(x_1,x_2,...,x_n)^T\) is the column vector of the variables. Then the theorem states that in this case the system has a unique solution, whose individual values for the unknowns are given by:

\(x_i=\frac{det(A_i)}{det(A)},~i=1,2,...,n\) (equation no.1)

where \(A_i\) is the matrix formed by replacing the i-th column of \(A\) by the column vector \(b\).

Implementation of a Sample

Lets create a random n × n matrix, \(A\) and a vector \(b\). To do this, we assign 2 to n and generate \(n^2\) numbers of random values using uniform distribution on the interval \([0,1]\). Then we create matrix \(A\) by the matrix function. For visualization simplicity we round the numbers to 2 digit decimal places using round() function. For construction of the vector \(b\) only assigning a random values and transpose it, would be sufficient.

n <- 2
values <- round(runif(n*n), 2)
A <-  matrix(ncol = n, nrow = n, data = values, byrow = T)
b <-  t(round(runif(n), 2))

# Some naming sugars
colnames(A) <- paste0(rep("x", n), 1:n)
row.names(A) <- paste0(rep("Eq", n), 1:n)
colnames(b) <- paste0(rep("b", n), 1:n)
row.names(b) <- ""

Thus the matrix \(A\) would be:

##       x1   x2
## Eq1 0.79 0.19
## Eq2 0.07 0.60

and the vector \(b\) would be:

##    b1   b2
##  0.16 0.85

For solving the equation \(Ax=b\) construction of the matrices \(A_1\) and \(A_2\) are necessary. We write the code in a way that makes us easy to change the size of matrix \(A\) (the value n) for further Samples. Thus we crate a variable Ai as list. Its clear that for constructing \(A_i\)s for \(i=1,2,...,n\) we need n loops.

Ai <- list()

for (i in 1:n) {
  A_i <- A
  A_i[, i] <- b
  Ai[[paste0("A", i)]] <- A_i
}

Thus the list \(A\) would be:

## $A1
##       x1   x2
## Eq1 0.16 0.19
## Eq2 0.85 0.60
## 
## $A2
##       x1   x2
## Eq1 0.79 0.16
## Eq2 0.07 0.85

For solving the system, we need to calculate equation no.1 for every value of i. Thus we need n loops. But as the denominator stays the same for every value i in the formula, we can calculate it once and use it throw the loops.

x <- list()

det_A <- det(A)

for (i in 1:n) {
  x[[paste0("x", i)]] <- det(Ai[[paste0("A", i)]]) / det_A
}

Consequently the list x contains the results:

## $x1
## [1] -0.142175
## 
## $x2
## [1] 1.433254

Generalization of the algorithm

We can compact the written codes in a function as bellow:

CR <- function(Ab){
  
  A <- Ab$A
  b <- Ab$b
  
  # input values A and b dimension checking
  if (!((dim(A)[1] == dim(A)[2]) && dim(A)[1] == length(b))) {
    stop("Input value dimension error.")
  }
  
  n <- length(b)
  
  Ai <- list()
  
  for (i in 1:n) {
    A_i <- A
    A_i[, i] <- b
    Ai[[paste0("A", i)]] <- A_i
  }
  
  x <- list()

  for (i in 1:n) {
    x[[paste0("x", i)]] <- det(Ai[[paste0("A", i)]]) / det_A
  }
  
  return(x)
}

Checking the function CR() by giving the value A and b as inputs should return the previews calculation result x.

CR(list(A=A, b=b))

## $x1
## [1] -0.142175
## 
## $x2
## [1] 1.433254

Benchmarking and modeling

In order to evaluate the effect of different matrix sizes on CR() function, first lets create a function random.A.b() which generates random variables A and b by given a size n.

random.A.b <- function(n){
  
  values <- round(runif(n*n), 2)
  A <-  matrix(ncol = n, nrow = n, data = values, byrow = T)
  b <-  t(round(runif(n), 2))
  
  return(list(A=A, b=b))
}

We declare variable sizes, in order to apply the mentioned function CR() and record the time consumed for calculation in milliseconds. For increasing statistical validity, every measurement will be repeated 10 times for each benchmarks. The simple average of every sample for fixed size matrix will be kept in vector means.
Then we use lm() function to create a linear model of the result. The formula for this model is assumed to be a polynomial by degree of 4.

library(microbenchmark)
means <- c()
sizes <- (1:25)*10
i <- 1

for (size in sizes) {
  data <- random.A.b(size)
  means[i] <- mean(microbenchmark(CR(data), times = 10)$time) / 10^6
  i <- i+1
} 


model <- lm(means ~ sizes + I(sizes^2) + I(sizes^3) + I(sizes^4))
c <- as.numeric(coefficients(model))
y <- c[1] + c[2]*sizes + c[3]*sizes^2 + c[4]*sizes^3 + c[5]*sizes^4

library(ggplot2)
ggplot() +
  geom_point(aes(sizes, means), color="black", na.rm = T) +
  geom_line(aes(sizes, y), col = "red") +
  labs(x="Matrix A's Size (n)", y="Time (ms)", 
       title = paste("Benchmarking of", length(sizes), "samples"))

Complexity order

The the most expensive computational operation in the implemented CR() function is the det() function. R-base uses Bareiss algorithm for calculating the determinant of a square matrices. Which has the time complexity order of \(O(n^3)\).
Determinant operation is done n times in order to calculate \(det(A_i);i=1,2,...,n\) and only one time for \(det(A)\). Thus the CR() algorithm time complexity function would be:

\(T_1(n)=(n+1)*n^3=n^4+n^3\)
\(O(T_1(n)) = O(n^4)\)

However, there are some other methods of determinant calculation. The simplest one is the Laplace expansion (aks cofactor expansion) which has the time complexity order of \(O(n!)\). Therefore if we had used this method in the CR() function for calculation, the time complexity order would be:

\(T_2(n)=(n+1)*n!\)
\(O(T_2(n)) = O(n!)\)

---

resources:

[1] https://en.wikipedia.org/wiki/Cramer%27s_rule#cite_note-11
[2] https://en.wikipedia.org/wiki/Laplace_expansion
[3] https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations