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 are the coefficients of the objective function (i.e., costs), are the coefficients of the constraints, and are the right-hand side of the constraints.

The simplex algorithm to solve the problem consists in the next steps:

1. If for all , then the current solution is optimal. Basic variables are equal to and non-basic variables are equal to 0.
2. If for some , we choose it to enter the base. We chose the variable with the most negative , let’s say that it is .
3. If for all , then the problem is unbounded.
4. If for some , we choose such that and pivot on , to then go back to step 1.

The coefficients are updated according to:

1. for
2. 
3. for
4. 
5. 
6. 

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 and the optimal value 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())
    {
        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;
        }
    }
}

## 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

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)

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 !