Implementing the Simplex Method in the Matrix Way

Sandipan Dey

Introduction

The goal of this article is to explain, implement and demonstrate how the popular simplex algorithm for solving a linear program can be revised to form a matrix-based algorithm (instead of tableau-based), so that the algorithm is known as revised simplex algorithm. This algorithm has been taught in the coursera course Operations Research Theory by the National Taiwan University. The theory part will be described as is from the course itself, along with an attempt to implement the matrix-based simplex method with R (from scratch) and demonstrate with examples.

Theory

Let’s start with the theory part. The simplex method is a very popular algorithm to solve constrained linear optimization (maximization or minimization) problems. It can be thought of as optimizing a linear objective function of the decision variables, so that the (linear) constraints remain satisfied, by traveling in the direction of increase in the value of the objective function (or equivalently decrease in the value of objective, in case of a minimization problem), from one corner to the another corner of a convex polytope formed by the constraints (each corner represents a basic feasible solution), until the optimal solution is obtained (assuming the problem is bounded and feasible) and stopping as soon as the optimal solution is obtained.

The iterative algorithm starts with a basic feasible solution and with each iteration it improves (e.g., for a maximization problem increases) the objective function value by changing a nonbasic variable to a basic one. First the linear program needs to be converted to the standard form. At each iteration the nonbasic variable corresponding to the most negative value in the reduced costs (allowing the highest possible increase) enters the basis, which constitutes the pivot column, Now, the basic variable which is first supposed to violate a constraint leaves the basis, determined by the minimum non-negative finite value in the ratio test. The algorithm terminates when all values are non-negative in the reduced costs, since there is no further scope of increase in value for the objective function.

The next figure describes how the simplex method can be thought of in matrix way. Note that the classic tableau can be computed using the matrix method at any step, although it’s not required.

Implementation

The next R code snippet implements the above algorithm in matrix way, at every iterative step of the simplex method, the function simplex.matrix.iter() needs to get called. Note that we only need to consider the pivot column to determine which basic variable leaves the basis, using the ratio test and not the entier matrix.

library(tidyverse)
library(knitr)

print.matrices <- function(A, b, c, X_B, X_N, A_B, A_N, c_B, c_N, x, z, r_c, opt) {

  X <- 1:(ncol(A))
  dimnames(A) <- list(NULL, paste0('x', X)) #[[2]]
  print(A %>%  kable(caption='**A**'))
  print(b %>% kable(caption='**b**')) 
  c_t <- matrix(c, nrow=1)
  dimnames(c_t)[[2]] <- paste0('x', X)
  print(c_t %>% kable(caption='$$c^T$$'))
  X_B_t <- matrix(X_B, nrow=1)
  dimnames(X_B_t)[[2]] <- paste0('x', X_B)
  print(X_B_t %>% kable(caption='$$X_B$$'))
  X_N_t <- matrix(X_N, nrow=1)
  dimnames(X_N_t)[[2]] <- paste0('x', X_N)
  print(X_N_t %>% kable(caption='$$X_N$$'))
  c_B_t <- matrix(c_B, nrow=1)
  dimnames(c_B_t)[[2]] <- paste0('x', X_B)
  print(c_B_t %>% kable(caption='$$c_B^T$$'))
  c_N_t <- matrix(c_N, nrow=1)
  dimnames(c_N_t)[[2]] <- paste0('x', X_N)
  print(c_N_t %>% kable(caption='$$c_N^T$$'))
  dimnames(A_B) <- list(NULL, paste0('x', X_B))
  dimnames(A_N) <- list(NULL, paste0('x', X_N))
  print(A_B %>% kable(caption='$$A_B$$'))
  print(A_N %>% kable(caption='$$A_N$$'))
  x_t <- matrix(x, nrow=1)
  dimnames(x_t)[[2]] <- paste0('x', X)
  print(x_t %>% kable(caption='**x (solution)**'))
  print(z %>% kable(caption='**z (objective value)**'))
  print(opt %>% kable(caption='**Optimal**'))
  r_c <- matrix(r_c, nrow=1)
  dimnames(r_c) <- list(NULL, paste0('x', X_N))
  print(r_c %>% kable(caption='**reduced cost** = $$c_B^TA_B^{-1}A_N-c_N^T$$'))
}

simplex.matrix.iter <- function(A, b, c, X_B) {
  X <- 1:(ncol(A))
  X_N <- setdiff(X, X_B)
  A_B <- A[,X_B]
  A_N <- A[,X_N]
  c_B <- c[X_B]
  c_N <- c[X_N]
  A_B_1 <- solve(A_B)
  #A_B_1A <- A_B_1 %*% A
  A_B_1b <- solve(A_B, b)
  c_BA_B_1 <- c_B %*% A_B_1
  c_BA_B_1A_N_c_N <- c_BA_B_1%*%A_N - c_N
  x <- rep(0, ncol(A))
  x[X_B] <- A_B_1b
  z <- c_BA_B_1 %*% b
  opt <- all(c_BA_B_1A_N_c_N >= 0) #all(c(c_BB_1, c_BB_1A_c) >= 0)
  print.matrices(A, b, c, X_B, X_N, A_B, A_N, c_B, c_N, x, z, c_BA_B_1A_N_c_N, opt)

  res <- list()
  m <- min(c_BA_B_1A_N_c_N)
  if (m < 0) {
    e_ix <- X_N[which.min(c_BA_B_1A_N_c_N)]
    A_B_1A_e <- A_B_1 %*% A[,e_ix]
    m <- Inf
    l_ix <- -1
    for (i in 1:length(A_B_1b)) {
       if ((A_B_1A_e[i] > 0) & (A_B_1b[i] / A_B_1A_e[i] < m)) {
         m <- A_B_1b[i] / A_B_1A_e[i]
         l_ix <- i
      }
    }
    if (l_ix < 0) {
      res$unbounded <- TRUE
      return(res)
    }
    print(data.frame(step=paste0('x', e_ix, ' enters the basis')) %>% kable())
    print(A_B_1A_e %>%  kable(caption=paste0('$$A_B^{-1}A$$', e_ix)))
    print(A_B_1b %>%  kable(caption='$$A_B^{-1}b$$'))
    r <- A_B_1b / A_B_1A_e
    dimnames(r) <- list(paste0('x', X_B), NULL)
    print(r %>% kable(caption='**ratio**'))
    print(data.frame(step=paste0('x', X_B[l_ix], ' leaves the basis')) %>% kable())
    X_B[l_ix] <- e_ix
  } else {
    res$opt <- TRUE
  }
  res$X_B <- X_B  # vars
  res$x <- x     # values
  res$z <- z    # objective
  res$opt <- opt
  res$unbounded <- FALSE
  return(res)
}

Examples

Now, let’s consider a few examples to demonstrate the above algorithm implementation step by step.

Example 1

Using the above theory, here’s how the simplex method works, starting from the basis \[X_B=(x_1,x_4,x_5)\]:

Now, let’s obtain the same result using the implementation above:

A <- matrix(c(2,2,0,-1,1,1), ncol=2) 
m <- nrow(A) 
n <- ncol(A) 
b <- matrix(c(4,8,3), ncol=1) 
c <- c(1,0)

A <- cbind(A, diag(m)) # add slack columns
c <- c(c, rep(0, m))

X_B <- c(1,4,5)
opt <- FALSE
i <- 1
while (!opt) { 
  print(data.frame(Iteration=i) %>% kable()) 
  res <- simplex.matrix.iter(A, b, c, X_B) 
  X_B <- res$X_B 
  opt <- res$opt
  #print(res)
  i <- i + 1
}

Iteration
1

A

x1	x2	x3	x4	x5
2	-1	1	0	0
2	1	0	1	0
0	1	0	0	1

b

4

8

3

c^T

x1	x2	x3	x4	x5
1	0	0	0	0

X_B

x1	x4	x5
1	4	5

X_N

x2	x3
2	3

c_B^T

x1	x4	x5
1	0	0

c_N^T

x2	x3
0	0

A_B

x1	x4	x5
2	0	0
2	1	0
0	0	1

A_N

x2	x3
-1	1
1	0
1	0

x (solution)

x1	x2	x3	x4	x5
2	0	0	4	3

z (objective value)

2

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x2	x3
-0.5	0.5

step
x2 enters the basis

A_B⁻¹A2

-0.5

2.0

1.0

A_B⁻¹b

2

4

3

ratio

x1	-4
x4	2
x5	3

step
x4 leaves the basis

Iteration
2

A

x1	x2	x3	x4	x5
2	-1	1	0	0
2	1	0	1	0
0	1	0	0	1

b

4

8

3

c^T

x1	x2	x3	x4	x5
1	0	0	0	0

X_B

x1	x2	x5
1	2	5

X_N

x3	x4
3	4

c_B^T

x1	x2	x5
1	0	0

c_N^T

x3	x4
0	0

A_B

x1	x2	x5
2	-1	0
2	1	0
0	1	1

A_N

x3	x4
1	0
0	1
0	0

x (solution)

x1	x2	x3	x4	x5
3	2	0	0	1

z (objective value)

3

Optimal

x
TRUE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x3	x4
0.25	0.25

Example 2

This time let’s start with \[X_B=(x_3,x_4)\], i.e. the slack variables \[(s_1,s_2)\] as the basic variables.

Now, lets demonstrate with our implementation with the simplex matrix method:

A <- matrix(c(1,1,1,2), ncol=2) 
m <- nrow(A) 
n <- ncol(A) 
b <- matrix(c(4,6), ncol=1) 
c <- c(2,3)

A <- cbind(A, diag(m)) # add slack columns
c <- c(c, rep(0, m))

X_B <- (n+1):(n+m)
opt <- FALSE
i <- 1
while (!opt) { 
  print(data.frame(Iteration=i) %>% kable()) 
  res <- simplex.matrix.iter(A, b, c, X_B) 
  X_B <- res$X_B 
  opt <- res$opt
  #print(res)
  i <- i + 1
}

Iteration
1

A

x1	x2	x3	x4
1	1	1	0
1	2	0	1

b

4

6

c^T

x1	x2	x3	x4
2	3	0	0

X_B

x3	x4
3	4

X_N

x1	x2
1	2

c_B^T

x3	x4
0	0

c_N^T

x1	x2
2	3

A_B

x3	x4
1	0
0	1

A_N

x1	x2
1	1
1	2

x (solution)

x1	x2	x3	x4
0	0	4	6

z (objective value)

0

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x1	x2
-2	-3

step
x2 enters the basis

A_B⁻¹A2

1

2

A_B⁻¹b

4

6

ratio

x3	4
x4	3

step
x4 leaves the basis

Iteration
2

A

x1	x2	x3	x4
1	1	1	0
1	2	0	1

b

4

6

c^T

x1	x2	x3	x4
2	3	0	0

X_B

x3	x2
3	2

X_N

x1	x4
1	4

c_B^T

x3	x2
0	3

c_N^T

x1	x4
2	0

A_B

x3	x2
1	1
0	2

A_N

x1	x4
1	0
1	1

x (solution)

x1	x2	x3	x4
0	3	1	0

z (objective value)

9

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x1	x4
-0.5	1.5

step
x1 enters the basis

A_B⁻¹A1

0.5

0.5

A_B⁻¹b

1

3

ratio

x3	2
x2	6

step
x3 leaves the basis

Iteration
3

A

x1	x2	x3	x4
1	1	1	0
1	2	0	1

b

4

6

c^T

x1	x2	x3	x4
2	3	0	0

X_B

x1	x2
1	2

X_N

x3	x4
3	4

c_B^T

x1	x2
2	3

c_N^T

x3	x4
0	0

A_B

x1	x2
1	1
1	2

A_N

x3	x4
1	0
0	1

x (solution)

x1	x2	x3	x4
2	2	0	0

z (objective value)

10

Optimal

x
TRUE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x3	x4
1	1

Verify the soltion obtained with R package’s implementation:

library(lpSolve)
objective.fn <- c
const.mat <- A
const.dir <- rep("<=", nrow(A))
const.rhs <- b
lp.solution <- lp("max", objective.fn, const.mat, 
                  const.dir, const.rhs, compute.sens=TRUE)
#lp.solution$objective
print(lp.solution$solution)

[1] 2 2 0 0

print(lp.solution$objval)

[1] 10

Example 3

Again, let’s start with \[X_B=(x_3,x_4,x_5)\], i.e. the slack variables as the basic variables.

A <- matrix(c(-1,-1,3,1,2,1), ncol=2) 
m <- nrow(A) 
n <- ncol(A) 
b <- matrix(c(3,8,18), ncol=1) 
c <- c(1,3)

A <- cbind(A, diag(m)) # add slack columns
c <- c(c, rep(0, m))

X_B <- (n+1):(n+m)
opt <- FALSE
i <- 1
while (!opt) { 
  print(data.frame(Iteration=i) %>% kable()) 
  res <- simplex.matrix.iter(A, b, c, X_B) 
  X_B <- res$X_B 
  opt <- res$opt
  #print(res)
  i <- i + 1
}

Iteration
1

A

x1	x2	x3	x4	x5
-1	1	1	0	0
-1	2	0	1	0
3	1	0	0	1

b

3

8

18

c^T

x1	x2	x3	x4	x5
1	3	0	0	0

X_B

x3	x4	x5
3	4	5

X_N

x1	x2
1	2

c_B^T

x3	x4	x5
0	0	0

c_N^T

x1	x2
1	3

A_B

x3	x4	x5
1	0	0
0	1	0
0	0	1

A_N

x1	x2
-1	1
-1	2
3	1

x (solution)

x1	x2	x3	x4	x5
0	0	3	8	18

z (objective value)

0

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x1	x2
-1	-3

step
x2 enters the basis

A_B⁻¹A2

1

2

1

A_B⁻¹b

3

8

18

ratio

x3	3
x4	4
x5	18

step
x3 leaves the basis

Iteration
2

A

x1	x2	x3	x4	x5
-1	1	1	0	0
-1	2	0	1	0
3	1	0	0	1

b

3

8

18

c^T

x1	x2	x3	x4	x5
1	3	0	0	0

X_B

x2	x4	x5
2	4	5

X_N

x1	x3
1	3

c_B^T

x2	x4	x5
3	0	0

c_N^T

x1	x3
1	0

A_B

x2	x4	x5
1	0	0
2	1	0
1	0	1

A_N

x1	x3
-1	1
-1	0
3	0

x (solution)

x1	x2	x3	x4	x5
0	3	0	2	15

z (objective value)

9

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x1	x3
-4	3

step
x1 enters the basis

A_B⁻¹A1

-1

1

4

A_B⁻¹b

3

2

15

ratio

x2	-3.00
x4	2.00
x5	3.75

step
x4 leaves the basis

Iteration
3

A

x1	x2	x3	x4	x5
-1	1	1	0	0
-1	2	0	1	0
3	1	0	0	1

b

3

8

18

c^T

x1	x2	x3	x4	x5
1	3	0	0	0

X_B

x2	x1	x5
2	1	5

X_N

x3	x4
3	4

c_B^T

x2	x1	x5
3	1	0

c_N^T

x3	x4
0	0

A_B

x2	x1	x5
1	-1	0
2	-1	0
1	3	1

A_N

x3	x4
1	0
0	1
0	0

x (solution)

x1	x2	x3	x4	x5
2	5	0	0	7

z (objective value)

17

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x3	x4
-5	4

step
x3 enters the basis

A_B⁻¹A3

-1

-2

7

A_B⁻¹b

5

2

7

ratio

x2	-5
x1	-1
x5	1

step
x5 leaves the basis

Iteration
4

A

x1	x2	x3	x4	x5
-1	1	1	0	0
-1	2	0	1	0
3	1	0	0	1

b

3

8

18

c^T

x1	x2	x3	x4	x5
1	3	0	0	0

X_B

x2	x1	x3
2	1	3

X_N

x4	x5
4	5

c_B^T

x2	x1	x3
3	1	0

c_N^T

x4	x5
0	0

A_B

x2	x1	x3
1	-1	1
2	-1	0
1	3	0

A_N

x4	x5
0	0
1	0
0	1

x (solution)

x1	x2	x3	x4	x5
4	6	1	0	0

z (objective value)

22

Optimal

x
TRUE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x4	x5
1.142857	0.7142857

Let’s compare the above solution obtained using the function lp() from R package lpSolve. Note that exactly same solution is obtained.

library(lpSolve)
objective.fn <- c
const.mat <- A
const.dir <- rep("<=", nrow(A))
const.rhs <- b
lp.solution <- lp("max", objective.fn, const.mat, 
                  const.dir, const.rhs, compute.sens=TRUE)
print(lp.solution$solution)

[1] 4 6 0 0 0

print(lp.solution$objval)

[1] 22

Example 4

A <- matrix(c(1,2,0,0,2,1,3,0,5,1,-3,2), nrow=3) 
m <- nrow(A) 
n <- ncol(A) 
b <- matrix(c(15,18,20), ncol=1) 
c <- c(4,3,2,3) 

A <- cbind(A, diag(m)) # add slack columns
c <- c(c, rep(0, m))
X_B <- (n+1):(n+m)

opt <- FALSE
i <- 1
while (!opt) { 
  print(data.frame(Iteration=i) %>% kable()) 
  res <- simplex.matrix.iter(A, b, c, X_B) 
  X_B <- res$X_B 
  opt <- res$opt
  #print(res)
  i <- i + 1
}

Iteration
1

A

x1	x2	x3	x4	x5	x6	x7
1	0	3	1	1	0	0
2	2	0	-3	0	1	0
0	1	5	2	0	0	1

b

15

18

20

c^T

x1	x2	x3	x4	x5	x6	x7
4	3	2	3	0	0	0

X_B

x5	x6	x7
5	6	7

X_N

x1	x2	x3	x4
1	2	3	4

c_B^T

x5	x6	x7
0	0	0

c_N^T

x1	x2	x3	x4
4	3	2	3

A_B

x5	x6	x7
1	0	0
0	1	0
0	0	1

A_N

x1	x2	x3	x4
1	0	3	1
2	2	0	-3
0	1	5	2

x (solution)

x1	x2	x3	x4	x5	x6	x7
0	0	0	0	15	18	20

z (objective value)

0

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x1	x2	x3	x4
-4	-3	-2	-3

step
x1 enters the basis

A_B⁻¹A1

1

2

0

A_B⁻¹b

15

18

20

ratio

x5	15
x6	9
x7	Inf

step
x6 leaves the basis

Iteration
2

A

x1	x2	x3	x4	x5	x6	x7
1	0	3	1	1	0	0
2	2	0	-3	0	1	0
0	1	5	2	0	0	1

b

15

18

20

c^T

x1	x2	x3	x4	x5	x6	x7
4	3	2	3	0	0	0

X_B

x5	x1	x7
5	1	7

X_N

x2	x3	x4	x6
2	3	4	6

c_B^T

x5	x1	x7
0	4	0

c_N^T

x2	x3	x4	x6
3	2	3	0

A_B

x5	x1	x7
1	1	0
0	2	0
0	0	1

A_N

x2	x3	x4	x6
0	3	1	0
2	0	-3	1
1	5	2	0

x (solution)

x1	x2	x3	x4	x5	x6	x7
9	0	0	0	6	0	20

z (objective value)

36

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x2	x3	x4	x6
1	-2	-9	2

step
x4 enters the basis

A_B⁻¹A4

2.5

-1.5

2.0

A_B⁻¹b

6

9

20

ratio

x5	2.4
x1	-6.0
x7	10.0

step
x5 leaves the basis

Iteration
3

A

x1	x2	x3	x4	x5	x6	x7
1	0	3	1	1	0	0
2	2	0	-3	0	1	0
0	1	5	2	0	0	1

b

15

18

20

c^T

x1	x2	x3	x4	x5	x6	x7
4	3	2	3	0	0	0

X_B

x4	x1	x7
4	1	7

X_N

x2	x3	x5	x6
2	3	5	6

c_B^T

x4	x1	x7
3	4	0

c_N^T

x2	x3	x5	x6
3	2	0	0

A_B

x4	x1	x7
1	1	0
-3	2	0
2	0	1

A_N

x2	x3	x5	x6
0	3	1	0
2	0	0	1
1	5	0	0

x (solution)

x1	x2	x3	x4	x5	x6	x7
12.6	0	0	2.4	0	0	15.2

z (objective value)

57.6

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x2	x3	x5	x6
-2.6	8.8	3.6	0.2

step
x2 enters the basis

A_B⁻¹A2

-0.4

0.4

1.8

A_B⁻¹b

2.4

12.6

15.2

ratio

x4	-6.000000
x1	31.500000
x7	8.444444

step
x7 leaves the basis

Iteration
4

A

x1	x2	x3	x4	x5	x6	x7
1	0	3	1	1	0	0
2	2	0	-3	0	1	0
0	1	5	2	0	0	1

b

15

18

20

c^T

x1	x2	x3	x4	x5	x6	x7
4	3	2	3	0	0	0

X_B

x4	x1	x2
4	1	2

X_N

x3	x5	x6	x7
3	5	6	7

c_B^T

x4	x1	x2
3	4	3

c_N^T

x3	x5	x6	x7
2	0	0	0

A_B

x4	x1	x2
1	1	0
-3	2	2
2	0	1

A_N

x3	x5	x6	x7
3	1	0	0
0	0	1	0
5	0	0	1

x (solution)

x1	x2	x3	x4	x5	x6	x7
9.222222	8.444444	0	5.777778	0	0	0

z (objective value)

79.55556

Optimal

x
TRUE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x3	x5	x6	x7
12.55556	2.444444	0.7777778	1.444444

Verification with lpSolve:

library(lpSolve)
objective.fn <- c
const.mat <- A
const.dir <- rep("<=", nrow(A))
const.rhs <- b
lp.solution <- lp("max", objective.fn, const.mat, 
                  const.dir, const.rhs, compute.sens=TRUE)
#lp.solution$objective
print(lp.solution$solution)

[1] 9.222222 8.444444 0.000000 5.777778 0.000000 0.000000 0.000000

print(lp.solution$objval)

[1] 79.55556

Example 5

A <- matrix(c(2,1,2,1,2,2,1,3,1), nrow=3) 
m <- nrow(A) 
n <- ncol(A) 
b <- matrix(c(2,5,6), ncol=1) 
c <- c(3,1,3) 

A <- cbind(A, diag(m)) # add slack columns
c <- c(c, rep(0, m))
X_B <- (n+1):(n+m)

opt <- FALSE
i <- 1
while (!opt) { 
  print(data.frame(Iteration=i) %>% kable()) 
  res <- simplex.matrix.iter(A, b, c, X_B) 
  X_B <- res$X_B 
  opt <- res$opt
  #print(res)
  i <- i + 1
}

Iteration
1

A

x1	x2	x3	x4	x5	x6
2	1	1	1	0	0
1	2	3	0	1	0
2	2	1	0	0	1

b

2

5

6

c^T

x1	x2	x3	x4	x5	x6
3	1	3	0	0	0

X_B

x4	x5	x6
4	5	6

X_N

x1	x2	x3
1	2	3

c_B^T

x4	x5	x6
0	0	0

c_N^T

x1	x2	x3
3	1	3

A_B

x4	x5	x6
1	0	0
0	1	0
0	0	1

A_N

x1	x2	x3
2	1	1
1	2	3
2	2	1

x (solution)

x1	x2	x3	x4	x5	x6
0	0	0	2	5	6

z (objective value)

0

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x1	x2	x3
-3	-1	-3

step
x1 enters the basis

A_B⁻¹A1

2

1

2

A_B⁻¹b

2

5

6

ratio

x4	1
x5	5
x6	3

step
x4 leaves the basis

Iteration
2

A

x1	x2	x3	x4	x5	x6
2	1	1	1	0	0
1	2	3	0	1	0
2	2	1	0	0	1

b

2

5

6

c^T

x1	x2	x3	x4	x5	x6
3	1	3	0	0	0

X_B

x1	x5	x6
1	5	6

X_N

x2	x3	x4
2	3	4

c_B^T

x1	x5	x6
3	0	0

c_N^T

x2	x3	x4
1	3	0

A_B

x1	x5	x6
2	0	0
1	1	0
2	0	1

A_N

x2	x3	x4
1	1	1
2	3	0
2	1	0

x (solution)

x1	x2	x3	x4	x5	x6
1	0	0	0	4	4

z (objective value)

3

Optimal

x
FALSE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x2	x3	x4
0.5	-1.5	1.5

step
x3 enters the basis

A_B⁻¹A3

0.5

2.5

0.0

A_B⁻¹b

1

4

4

ratio

x1	2.0
x5	1.6
x6	Inf

step
x5 leaves the basis

Iteration
3

A

x1	x2	x3	x4	x5	x6
2	1	1	1	0	0
1	2	3	0	1	0
2	2	1	0	0	1

b

2

5

6

c^T

x1	x2	x3	x4	x5	x6
3	1	3	0	0	0

X_B

x1	x3	x6
1	3	6

X_N

x2	x4	x5
2	4	5

c_B^T

x1	x3	x6
3	3	0

c_N^T

x2	x4	x5
1	0	0

A_B

x1	x3	x6
2	1	0
1	3	0
2	1	1

A_N

x2	x4	x5
1	1	0
2	0	1
2	0	0

x (solution)

x1	x2	x3	x4	x5	x6
0.2	0	1.6	0	0	4

z (objective value)

5.4

Optimal

x
TRUE

reduced cost = c_B^TA_B⁻¹A_N − c_N^T

x2	x4	x5
1.4	1.2	0.6

Verification with lpSolve:

library(lpSolve)
objective.fn <- c
const.mat <- A
const.dir <- rep("<=", nrow(A))
const.rhs <- b
lp.solution <- lp("max", objective.fn, const.mat, 
                  const.dir, const.rhs, compute.sens=TRUE)
#lp.solution$objective
print(lp.solution$solution)

[1] 0.2 0.0 1.6 0.0 0.0 0.0

print(lp.solution$objval)

[1] 5.4