A naive simplex phase 2 implementation with C++ 11 and R
R and Shiny Training: If you find this blog to be interesting, please note that I offer personalized and group-based training sessions that may be reserved through Buy me a Coffee. Additionally, I provide training services in the Spanish language and am available to discuss means by which I may contribute to your Shiny project.
Motivation
I am in the process of learning C++, so I created the Cpp11 + R examples repository.
Without covering the details of the simplex algorithm, I will focus on the implementation of the algorithm in C++ 11 and R.
Final code
You can find it here. In the rest of the post I explain the logic behind the code and how to build and test it.
Simplex (phase 2) algorithm
The simplex algorithm is well described in Introduction to Linear Optimization and there is efficient software to solve this. The goal here is to show how to print the intermediate steps of the algorithm, so that it can be used to practice the simplex algorithm for undergraduate students.
Here I present the recipe and and example in the next section.
A problem written in canonical form is represented by a table such as:
\[ \begin{array}{ccc|c} x_1 & \cdots & x_n & \\ \hline c_1 & \cdots & c_n & -z \\ a_{11} & \cdots & a_{1n} & b_1 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m1} & \cdots & a_{mn} & b_m \end{array} \]
where \(c_1, \ldots, c_n\) are the coefficients of the objective function (i.e., costs), \(a_{11}, \ldots, a_{mn}\) are the coefficients of the constraints, and \(b_1, \ldots, b_m\) are the right-hand side of the constraints.
The simplex algorithm to solve the problem consists in the next steps:
- If \(c_j \geq 0\) for all \(j\), then the current solution is optimal. Basic variables are equal to \(b_i\) and non-basic variables are equal to 0.
- If \(c_j < 0\) for some \(j\), we choose it to enter the base. We chose the variable with the most negative \(c_j\), let’s say that it is \(j = s\).
- If \(a_{is} \leq 0\) for all \(i\), then the problem is unbounded.
- If \(a_{is} > 0\) for some \(i\), we choose \(i = r\) such that \(\frac{b_r}{a_{rs}} = \min(\frac{b_i}{a_is},\: a_{is} 0)\) and pivot on \(a_{rs}\), to then go back to step 1.
The coefficients are updated according to:
- \(a_{ij} \leftarrow a_{ij} - \frac{a_{is} a_{rj}}{a_{rs}}\) for \(j \neq s\)
- \(a_{rj} \leftarrow \frac{a_{rj}}{a_{rs}}\)
- \(b_i \leftarrow b_i - \frac{a_{is} b_r}{a_{rs}}\) for \(i \neq r\)
- \(b_r \leftarrow \frac{b_r}{a_{rs}}\)
- \(c_j \leftarrow c_j - \frac{c_s a_{rj}}{a_{rs}}\)
- \(-z \leftarrow -z - \frac{c_s b_r}{a_{rs}}\)
This algorithm is equivalent to Gauss method to solve linear systems.
Simplex example
Let’s say we want to solve the following linear programming problem:
\[ \begin{aligned} \text{min} \quad & -x_1 - 3x_2 \\ \text{subject to} \quad & x_1 + x_2 \geq 3 \\ & -3x_1 + x_2 \geq 2 \\ & x_1, x_2 \geq 0 \end{aligned} \]
We need to re-write the problem in canonical form:
\[ \begin{aligned} \text{min} \quad & -x_1 - 3x_2 + 0x_3 + 0x_4 \\ \text{subject to} \quad & x_1 + x_2 + x_3 = 3 \\ & -3x_1 + x_2 + x_4 = 2 \\ & x_1, x_2,x_3,x_4 \geq 0 \end{aligned} \]
The initial tableau for this problem is:
\[ \begin{array}{cccc|c} x_1 & x_2 & x_3 & x_4 & -z \\ \hline -1 & -3 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 3 \\ -3 & 1 & 0 & 1 & 2 \end{array} \]
The first row is the cost row, the last column is the right-hand side, and the rest is the matrix \(A\).
We pivot on row 2, column 2:
\[ \begin{array}{cccc|c} x_1 & x_2 & x_3 & x_4 & -z \\ \hline -10 & 0 & 0 & 3 & 6 \\ 4 & 0 & 1 & -1 & 1 \\ -3 & 1 & 0 & 1 & 2 \end{array} \]
Then, we pivot on row 2, column 1:
\[ \begin{array}{cccc|c} x_1 & x_2 & x_3 & x_4 & -z \\ \hline 0 & 0 & 5/2 & 1/2 & 17/2 \\ 1 & 0 & 1/4 & -1/4 & 1/4 \\ 0 & 1 & 3/4 & 1/4 & 11/4 \end{array} \]
Here we reached a stopping criterion: the minimum cost is non-negative, so we have found an optimal solution, which is \(x^* = (\frac{1}{4}, \frac{11}{4}, 0 , 0)\) and the optimal value is \(z^* = -\frac{17}{2}\).
C++ 11
We will use the following C++ 11 code to implement the simplex algorithm. The code consists in writing the correct steps for the first pivot and then using a while loop to repeat the steps until the minimum cost is non-negative.
// doubles and matrix to do the math
// iostream to print the results
#include <cpp11/doubles.hpp>
#include <cpp11/matrix.hpp>
#include <iostream>
using namespace cpp11;
using namespace std;
// this function has no output
// i.e., it won't send a matrix or vector back to R
// instead, it will print the results to the console
[[cpp11::register]] void cpp11_simplex_phase2_(doubles c, doubles b, doubles_matrix<> A) {
// Naive implementation of the simplex algorithm (phase 2)
// Check dimensions
int m = A.nrow();
int n = A.ncol();
if (m != b.size())
{
stop("incompatible dimensions between A and b");
}
if (n != c.size())
{
stop("incompatible dimensions between A and c");
}
// Create the tableau
// DON'T FORGET TO USE writable:: TO MAKE IT MODIFIABLE
writable::doubles_matrix<> T(m + 1, n + m + 1);
// b2 = [0, b]
writable::doubles b2(m + 1);
// fill b2_i+1 with b_i
for (int i = 0; i < m; i++)
{
b2[i] = 0.0;
}
for (int i = 0; i < m; i++)
{
b2[i + 1] += b[i];
}
// Fill the tableau with the data
// fill with 0s
// we don't really need this but I don't like to see 1.0e-400 in the output
for (int i = 0; i < m + 1; i++)
{
for (int j = 0; j < n + m + 1; j++)
{
T(i, j) = 0.0;
}
}
// put b2 in the last column of T
// notice the 0-indexing here!!!
for (int i = 0; i < m + 1; i++)
{
T(i, n + m) += b2[i];
}
// put c in the first row of T
for (int j = 0; j < n; j++)
{
T(0, j) += c[j];
}
// put A in the lower-left corner of T
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
T(i + 1, j) += A(i, j);
}
}
// fill the diagonal of non-basic variables with 1s (for the identity matrix)
for (int i = 0; i < m; i++)
{
T(i + 1, n + i) += 1.0;
}
// Print the tableau
cout << "Initial tableau:" << endl;
// this won't work
// it prints a funky 0x556da90020e8 or similar
// cout << T << endl;
// this will work
// it prints the matrix in a nice format
for (int i = 0; i < m + 1; i++)
{
for (int j = 0; j < n + m + 1; j++)
{
// cout << T(i, j) << " ";
// if positive, print a space before or it will look misaligned
if (T(i, j) >= 0.0)
{
cout << " ";
}
cout << T(i, j) << " ";
}
cout << endl;
}
// STEP 1 OF THE SIMPLEX
// Cost criteria
// start with the first element in the first row (arbitrary)
double min_c_j = T(0, 0);
// replace min_c_j with the smallest number in the first row
for (int j = 0; j < n + m; j++)
{
if (T(0, j) < min_c_j)
{
min_c_j = T(0, j);
}
}
// Print the minimum cost
cout << "Minimum cost: " << min_c_j << endl;
// Find the pivot column = Find the most negative entry in the first row
int pivot_col = 0;
for (int j = 0; j < n + m; j++)
{
if (T(0, j) < T(0, pivot_col))
{
pivot_col = j;
}
}
// Find the pivow row = Find the row with the smallest positive ratio
// again, 0-indexing
int pivot_row = 1;
for (int i = 1; i < m + 1; i++)
{
if (T(i, n + m) / T(i, pivot_col) < T(pivot_row, n + m) / T(pivot_row, pivot_col))
{
// here we need a 2nd if as this can take a_rs < 0 for the pivot
// pivot_row = i;
// if the ratio is positive, then we can use it
if (T(i, n + m) / T(i, pivot_col) > 0) {
pivot_row = i;
}
}
}
// Print the pivot row and column
cout << "Pivot row: " << pivot_row << endl;
cout << "Pivot column: " << pivot_col + 1 << endl;
// Get the pivot element
double pivot = T(pivot_row, pivot_col);
// Divide the pivot row by the pivot element
for (int j = 0; j < n + m + 1; j++)
{
T(pivot_row, j) /= pivot;
}
// Subtract multiples of the pivot row from each row to make all other entries in the pivot column zero
for (int i = 0; i < m + 1; i++)
{
if (i != pivot_row)
{
double ratio = T(i, pivot_col);
for (int j = 0; j < n + m + 1; j++)
{
T(i, j) -= ratio * T(pivot_row, j);
}
}
}
// Print the new tableau
cout << "====\nNew tableau:" << endl;
for (int i = 0; i < m + 1; i++)
{
for (int j = 0; j < n + m + 1; j++)
{
if (T(i, j) >= 0.0)
{
cout << " ";
}
cout << T(i, j) << " ";
}
cout << endl;
}
// Cost criteria
// start with the first element in the first row (arbitrary)
min_c_j = T(0, 0);
// replace min_c_j with the smallest number in the first row
for (int j = 0; j < n + m; j++)
{
if (T(0, j) < min_c_j)
{
min_c_j = T(0, j);
}
}
// Print success message and stop if the minimum cost is non-negative
if (min_c_j >= 0)
{
cout << "Optimal solution found in 1 step!" << endl;
}
// Print the minimum cost
cout << "Minimum cost: " << min_c_j << endl;
// STEP 2, 3, ... OF THE SIMPLEX
// Repeat until the first row contains no negative numbers
// add a counter to keep track of the number of iterations
// we previously completed 1 iteration
int iter_count = 1;
// the computational complexity of this algorithm is O(mn^2)
// so we can use the dimensions of the tableau to set an upper bound on the number of iterations
double iter_bound = m * n * n;
while (min_c_j < 0)
{
if (iter_count > iter_bound)
{
stop("too many iterations");
}
// Find the pivot column = Find the most negative entry in the first row
pivot_col = 0;
for (int j = 0; j < n + m; j++)
{
if (T(0, j) < T(0, pivot_col))
{
pivot_col = j;
}
}
// Find the pivow row = Find the row with the smallest positive ratio
pivot_row = 1;
for (int i = 1; i < m + 1; i++)
{
if (T(i, n + m) / T(i, pivot_col) < T(pivot_row, n + m) / T(pivot_row, pivot_col))
{
if (T(i, n + m) / T(i, pivot_col) > 0) {
pivot_row = i;
}
}
}
// Print the pivot row and column
cout << "Pivot row: " << pivot_row << endl;
cout << "Pivot column: " << pivot_col + 1 << endl;
// Get the pivot element
pivot = T(pivot_row, pivot_col);
// Divide the pivot row by the pivot element
for (int j = 0; j < n + m + 1; j++)
{
T(pivot_row, j) /= pivot;
}
// Subtract multiples of the pivot row from each row to make all other entries in the pivot column zero
for (int i = 0; i < m + 1; i++)
{
if (i != pivot_row)
{
double ratio = T(i, pivot_col);
for (int j = 0; j < n + m + 1; j++)
{
T(i, j) -= ratio * T(pivot_row, j);
}
}
}
// Print the new tableau
cout << "====\nNew tableau:" << endl;
for (int i = 0; i < m + 1; i++)
{
for (int j = 0; j < n + m + 1; j++)
{
if (T(i, j) >= 0.0)
{
cout << " ";
}
cout << T(i, j) << " ";
}
cout << endl;
}
// Cost criteria
// start with the first element in the first row (arbitrary)
min_c_j = T(0, 0);
// replace min_c_j with the smallest number in the first row
for (int j = 0; j < n + m; j++)
{
if (T(0, j) < min_c_j)
{
min_c_j = T(0, j);
}
}
// Print the minimum cost
cout << "Minimum cost: " << min_c_j << endl;
// Increment the counter
iter_count++;
// Print success message and stop if the minimum cost is non-negative
if (min_c_j >= 0)
{
cout << "Optimal solution found in " << iter_count << " steps !" << endl;
break;
}
}
}R
We will use the following R code to call the C++ 11 function and print the results.
## usethis namespace: start
#' @useDynLib cpp11simplexphase2, .registration = TRUE
## usethis namespace: end
NULL
#' Print the table of the simplex algorithm
#' @param c vector of coefficients of the objective function
#' @param b vector of the right hand side of the constraints
#' @param A matrix of the coefficients of the constraints
#' @export
cpp11_simplex_phase2 <- function(c, b, A) {
cpp11_simplex_phase2_(c, b, A)
}Build and test
We can test with the previous example and the following R code:
# build
devtools::clean_dll()
cpp11::cpp_register()
devtools::document()
devtools::load_all()
# test
c <- c(-1, -3)
b <- c(3, 2)
A <- matrix(
c(1, -3, 1, 1),
nrow = 2,
ncol = 2,
byrow = FALSE
)
cpp11_simplex_phase2(c, b, A)The output shows the correct solution:
Initial tableau:
-1 -3 0 0 0
1 1 1 0 3
-3 1 0 1 2
Minimum cost: -3
Pivot row: 2
Pivot column: 2
====
New tableau:
-10 0 0 3 6
4 0 1 -1 1
-3 1 0 1 2
Minimum cost: -10
Pivot row: 1
Pivot column: 1
====
New tableau:
0 0 2.5 0.5 8.5
1 0 0.25 -0.25 0.25
0 1 0.75 0.25 2.75
Minimum cost: 0
Optimal solution found in 2 steps !