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:
where
The simplex algorithm to solve the problem consists in the next steps:
- If
for all , then the current solution is optimal. Basic variables are equal to and non-basic variables are equal to 0. - If
for some , we choose it to enter the base. We chose the variable with the most negative , let’s say that it is . - If
for all , then the problem is unbounded. - If
for some , we choose such that and pivot on , to then go back to step 1.
The coefficients are updated according to:
for for
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:
We need to re-write the problem in canonical form:
The initial tableau for this problem is:
The first row is the cost row, the last column is the right-hand side, and the rest is the matrix
We pivot on row 2, column 2:
Then, we pivot on row 2, column 1:
Here we reached a stopping criterion: the minimum cost is non-negative, so we have found an optimal solution, which is
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())
{
("incompatible dimensions between A and b");
stop}
if (n != c.size())
{
("incompatible dimensions between A and c");
stop}
// Create the tableau
// DON'T FORGET TO USE writable:: TO MAKE IT MODIFIABLE
::doubles_matrix<> T(m + 1, n + m + 1);
writable
// b2 = [0, b]
::doubles b2(m + 1);
writable
// fill b2_i+1 with b_i
for (int i = 0; i < m; i++)
{
[i] = 0.0;
b2}
for (int i = 0; i < m; i++)
{
[i + 1] += b[i];
b2}
// 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++)
{
(i, j) = 0.0;
T}
}
// put b2 in the last column of T
// notice the 0-indexing here!!!
for (int i = 0; i < m + 1; i++)
{
(i, n + m) += b2[i];
T}
// put c in the first row of T
for (int j = 0; j < n; j++)
{
(0, j) += c[j];
T}
// put A in the lower-left corner of T
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
(i + 1, j) += A(i, j);
T}
}
// fill the diagonal of non-basic variables with 1s (for the identity matrix)
for (int i = 0; i < m; i++)
{
(i + 1, n + i) += 1.0;
T}
// Print the tableau
<< "Initial tableau:" << endl;
cout
// 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 }
<< T(i, j) << " ";
cout }
<< endl;
cout }
// 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)
{
= T(0, j);
min_c_j }
}
// Print the minimum cost
<< "Minimum cost: " << min_c_j << endl;
cout
// 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))
{
= j;
pivot_col }
}
// 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) {
= i;
pivot_row }
}
}
// Print the pivot row and column
<< "Pivot row: " << pivot_row << endl;
cout << "Pivot column: " << pivot_col + 1 << endl;
cout
// 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++)
{
(pivot_row, j) /= pivot;
T}
// 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++)
{
(i, j) -= ratio * T(pivot_row, j);
T}
}
}
// Print the new tableau
<< "====\nNew tableau:" << endl;
cout for (int i = 0; i < m + 1; i++)
{
for (int j = 0; j < n + m + 1; j++)
{
if (T(i, j) >= 0.0)
{
<< " ";
cout }
<< T(i, j) << " ";
cout }
<< endl;
cout }
// Cost criteria
// start with the first element in the first row (arbitrary)
= T(0, 0);
min_c_j
// 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)
{
= T(0, j);
min_c_j }
}
// Print success message and stop if the minimum cost is non-negative
if (min_c_j >= 0)
{
<< "Optimal solution found in 1 step!" << endl;
cout }
// Print the minimum cost
<< "Minimum cost: " << min_c_j << endl;
cout
// 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)
{
("too many iterations");
stop}
// Find the pivot column = Find the most negative entry in the first row
= 0;
pivot_col for (int j = 0; j < n + m; j++)
{
if (T(0, j) < T(0, pivot_col))
{
= j;
pivot_col }
}
// Find the pivow row = Find the row with the smallest positive ratio
= 1;
pivot_row 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) {
= i;
pivot_row }
}
}
// Print the pivot row and column
<< "Pivot row: " << pivot_row << endl;
cout << "Pivot column: " << pivot_col + 1 << endl;
cout
// Get the pivot element
= T(pivot_row, pivot_col);
pivot
// Divide the pivot row by the pivot element
for (int j = 0; j < n + m + 1; j++)
{
(pivot_row, j) /= pivot;
T}
// 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++)
{
(i, j) -= ratio * T(pivot_row, j);
T}
}
}
// Print the new tableau
<< "====\nNew tableau:" << endl;
cout for (int i = 0; i < m + 1; i++)
{
for (int j = 0; j < n + m + 1; j++)
{
if (T(i, j) >= 0.0)
{
<< " ";
cout }
<< T(i, j) << " ";
cout }
<< endl;
cout }
// Cost criteria
// start with the first element in the first row (arbitrary)
= T(0, 0);
min_c_j
// 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)
{
= T(0, j);
min_c_j }
}
// Print the minimum cost
<< "Minimum cost: " << min_c_j << endl;
cout
// Increment the counter
++;
iter_count
// Print success message and stop if the minimum cost is non-negative
if (min_c_j >= 0)
{
<< "Optimal solution found in " << iter_count << " steps !" << endl;
cout 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
<- function(c, b, A) {
cpp11_simplex_phase2 cpp11_simplex_phase2_(c, b, A)
}
Build and test
We can test with the previous example and the following R code:
# build
::clean_dll()
devtools::cpp_register()
cpp11::document()
devtools::load_all()
devtools
# test
<- c(-1, -3)
c <- c(3, 2)
b
<- matrix(
A 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
: -3
Minimum cost: 2
Pivot row: 2
Pivot column====
:
New tableau-10 0 0 3 6
4 0 1 -1 1
-3 1 0 1 2
: -10
Minimum cost: 1
Pivot row: 1
Pivot column====
:
New tableau0 0 2.5 0.5 8.5
1 0 0.25 -0.25 0.25
0 1 0.75 0.25 2.75
: 0
Minimum costin 2 steps ! Optimal solution found