One year ago, I wrote an article for the now-defunct Microsoft MSDN Magazine titled “Roach Infestation Optimization”. The C# code associated with that article vanished into the void, so I figured I’d rehydrate the demo program.

Roach infestation optimization (RIO) is a numerical optimization algorithm that’s loosely based on the behavior of common cockroaches such as Periplaneta americana. Let me explain.

In machine learning, a numerical optimization algorithm is often used to find a set of values for variables (usually called weights) that minimize some measure of error. The process of determining the values of the weights is called training the model. The idea is to use a collection of training data that has known correct output values. A numerical optimization algorithm is used to find the values of the weights that minimize the error between computed output values and known correct output values.

The most common training algorithms are based on calculus derivatives (stochastic gradient descent), but there are also algorithms that are based on the behaviors of natural systems. These are sometimes called bio-inspired algorithms. These bio-inspired algorithms are usually not practical, but they’re interesting.

The goal of the demo program is to use RIO to find the minimum value of the Rastrigin function with eight input variables. The Rastrigin function is a standard benchmark function used to evaluate the effectiveness of numerical optimization algorithms. The function has a known minimum value of 0.0 located at x = (0, 0, . . 0) where the number of 0 values is equal to the number of input values.

The Rastrigin function is extremely difficult for most optimization algorithms because it has many peaks and valleys that create local minimum values that can trap the algorithms. It’s not possible to easily visualize the Rastrigin function with eight input values, but you can get a good idea of the function’s characteristics by examining a graph of the function for two input values:

The output of the demo is:

Begin roach infestation optimization demo

Goal is minimize Rastrigin's function in 8 dimensions
Problem known min value = 0.0 at (0, 0, .. 0)

Setting number of roaches  = 50
Setting maximum time = 10000

Starting roach optimization

time =    500 best error = 104.76530
time =   1000 best error = 67.36795
time =   1500 best error = 12.11031
time =   2000 best error = 8.90949
time =   2500 best error = 3.99195
Mass extinction at t =   2500
time =   3000 best error = 3.99195
time =   3500 best error = 1.80840
time =   4000 best error = 1.02166
time =   4500 best error = 0.99382
time =   5000 best error = 0.00000
Mass extinction at t =   5000
time =   5500 best error = 0.00000
time =   6000 best error = 0.00000
time =   6500 best error = 0.00000
time =   7000 best error = 0.00000
time =   7500 best error = 0.00000
Mass extinction at t =   7500
time =   8000 best error = 0.00000
time =   8500 best error = 0.00000
time =   9000 best error = 0.00000
time =   9500 best error = 0.00000

Roach algorithm completed

Best error found = 0.000000 at:

0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000

End roach infestation optimization demo

The demo program sets the number of roaches to 50. Each simulated roach has a position that represents a possible solution to the minimization problem. More roaches increase the chance of finding the true optimal solution at the expense of performance.

RIO is an iterative process and requires a maximum loop counter value. The demo sets the maximum value to 10,000 iterations. The maximum number of iterations will vary from problem to problem.

In the demo, the best (smallest) error associated with the best roach position found so far was displayed every 500 time units. After the algorithm finished, the best position found for any roach was x = (0, 0, 0, 0, 0, 0, 0, 0), which is, in fact, the correct answer. But notice if the maximum number of iterations had been set to 3,000 instead of 10,000, RIO would not have found the one global minimum. RIO is not guaranteed to find an optimal solution in practical scenarios.

RIO is essentially a variation of the similar particle swarm optimization (PSO). The main difference is that in RIO, roaches that are close to each other will sometimes exchange information with each other.

Bio-inspired algorithms such as roach infestation optimization and particle swarm optimization are usually too slow to be practical. But if/when quantum computing becomes a reality, these algorithms could become essential techniques.

“Terra Formars” (2016) is an absolutely bonkers Japanese science fiction movie. In the 21st century, Earth plans to colonize Mars. To terraform the planet, the surface of Mars is seeded with moss and cockroaches to spread the moss.

This movie is one of many things Japanese that absolutely baffle me.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor chokes on symbols).

using System;

// "Roach Infestation Optimization" (2008)
// IEEE Swarm Intelligence Symposium
// T. Havens, C. Spain, N. Salmon, J. Keller

namespace RoachOptimization
{
  class RoachOptimizationProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin roach infestation " +
        "optimization demo\n");

      int dim = 8;
      int numRoaches = 50;
      int tMax = 10000;
      Random rnd = new Random(1);
      
      Console.WriteLine("Goal is minimize Rastrigin's " +
        "function in " + dim + " dimensions");
      Console.WriteLine("Problem known min value = 0.0 " +
        "at (0, 0, .. 0) \n");
      Console.WriteLine("Setting number of roaches  = " +
        numRoaches);
      Console.WriteLine("Setting maximum time = " +
        tMax);
      
      Console.WriteLine("\nStarting roach optimization \n");
      double[] answer = 
        SolveRastrigin(dim, numRoaches, tMax, rnd);
      Console.WriteLine("\nRoach algorithm completed ");

      double err = Error(answer);
      Console.WriteLine("\nBest error found = " +
        err.ToString("F6") + " at: \n");

      for (int i = 0; i "lt" dim; ++i)
        Console.Write(answer[i].ToString("F4") + " ");
      Console.WriteLine("");

      Console.WriteLine("\nEnd roach infestation " +
        "optimization demo ");
      Console.ReadLine();
    } // Main

    // -----------------------------------------------------

    static double[] SolveRastrigin(int dim, int numRoaches,
      int tMax, Random rnd)
    {
      // estimate solution to Rastrigin's function
      // using roach infestation optimization
      double C0 = 0.7;  // RIO is based on PSO
      double C1 = 1.43;
      // probs of exchanging info based on number neighbors
      // actual roach behavior: 0.49, 0.63, 0.65
      // double[] A = new double[] { 0.49, 0.63, 0.65 }; 
      double[] A = new double[] { 0.2, 0.3, 0.4 }; // better 

      // smaller divisor (e.g., 5 instead of 10)
      // roaches stay longer
      int tHunger = tMax / 10;  
      double minX = -10.0;
      double maxX = 10.0;
      int extinction = tMax / 4;  // mass extinction restart
      //Random rnd = new Random(seed);  // random order

      Roach[] herd = new Roach[numRoaches];  // 'intrusion'
      for (int i = 0; i "lt" numRoaches; ++i)
        herd[i] = new Roach(dim, minX, maxX, tHunger, i);

      int t = 0;  // loop counter (time)
      int[] indices = new int[numRoaches];
      for (int i = 0; i "lt" numRoaches; ++i)
        indices[i] = i;

      double bestError = double.MaxValue;  // by any roach
      double[] bestPosition = new double[numRoaches];

      int displayInterval = tMax / 20; // 20 times

      // distances between all pairs of roaches
      double[][] distances = new double[numRoaches][];
      for (int i = 0; i "lt" numRoaches; ++i)
        distances[i] = new double[numRoaches];

      while (t "lt" tMax) // main loop
      {
        if (t "gt" 0 && t % displayInterval == 0)
        {
          Console.Write("time = " +
            t.ToString().PadLeft(6));
          Console.WriteLine(" best error = " +
            bestError.ToString("F5"));
        }

        // compute distances between all pairs
        for (int i = 0; i "lt" numRoaches - 1; ++i)
        {
          for (int j = i + 1; j "lt" numRoaches; ++j)
          {
            double d = 
              Distance(herd[i].position, herd[j].position);
            distances[i][j] = distances[j][i] = d;
          }
        }

        // find 'quarter' median distance
        double[] sortedDists = 
          new double[numRoaches * (numRoaches - 1) / 2];
        int k = 0;
        for (int i = 0; i "lt" numRoaches - 1; ++i)
        {
          for (int j = i + 1; j "lt" numRoaches; ++j)
          {
            sortedDists[k++] = distances[i][j];
          }
        }

        Array.Sort(sortedDists);
        double medianDist = 
          sortedDists[(int)(sortedDists.Length / 4)];

        Shuffle(indices, rnd); // process roaches rnd order

        for (int i = 0; i "lt" numRoaches; ++i)  // each roach
        {
          int idx = indices[i]; // roach index

          Roach curr = herd[idx]; // ref to current roach
          int numNeighbors = 0;

          // calculate number of neighbors of curr roach
          for (int j = 0; j "lt" numRoaches; ++j)  
          {
            if (j == idx) continue; // not a neighbor self
            double d = distances[idx][j];
            if (d "lt" medianDist) // is a neighbor to roach[j]
              ++numNeighbors;
          }

          // at most 3 neighbors are relevant
          int effectiveNeighbors = numNeighbors;
          if (effectiveNeighbors "gte" 3)  
            effectiveNeighbors = 3;

          // possibly swap group best information
          for (int j = 0; j "lt" numRoaches; ++j)
          {
            if (j == idx) continue; // not neighbor of self
            if (effectiveNeighbors == 0) continue;
            double prob = rnd.NextDouble();
            if (prob "gt" A[effectiveNeighbors - 1])
              continue; // don't always swap

            double d = distances[idx][j];
            if (d "lt" medianDist) // is a neighbor
            {
              // exchange info if curr roach better
              if (curr.error "lt" herd[j].error) // 
              {
                for (int p = 0; p "lt" dim; ++p)
                {
                  // use curr roach's personal best
                  // as the group best for both roaches
                  herd[j].groupBestPos[p] = 
                    curr.personalBestPos[p];  
                  curr.groupBestPos[p] = 
                    curr.personalBestPos[p];
                }
              }
              else // roach[j] is better than curr
              {
                for (int p = 0; p "lt" dim; ++p)
                {
                  curr.groupBestPos[p] = 
                    herd[j].personalBestPos[p];
                  herd[j].groupBestPos[p] = 
                    herd[j].personalBestPos[p];
                }
              }

            } // if a neighbor
          } // j, each neighbor

          if (curr.hunger "lt" tHunger) // move roach
          {
            // new velocity
            for (int p = 0; p "lt" dim; ++p)
              curr.velocity[p] = 
                (C0 * curr.velocity[p]) +
                (C1 * rnd.NextDouble() *
                (curr.personalBestPos[p] - 
                curr.position[p])) +
                (C1 * rnd.NextDouble() *
                (curr.groupBestPos[p] - 
                curr.position[p]));

            // update position
            for (int p = 0; p "lt" dim; ++p)
              curr.position[p] = 
                curr.position[p] + curr.velocity[p];

            double e = Error(curr.position); // Rastrigin
            curr.error = e;

            // new personal best position?
            if (curr.error "lt" curr.personalBestErr)
            {
              curr.personalBestErr = curr.error;
              for (int p = 0; p "lt" dim; ++p)
                curr.personalBestPos[p] = curr.position[p];
            }

            // new global best position?
            if (curr.error "lt" bestError)
            {
              bestError = curr.error;
              for (int p = 0; p "lt" dim; ++p)
                bestPosition[p] = curr.position[p];
            }

            ++curr.hunger;
          }
          else // hungry. leave group to random position
          {
            herd[idx] = 
              new Roach(dim, minX, maxX, tHunger, t);
          }

        } // each roach

        // mass extinction?
        if (t "gt" 0 && t % extinction == 0)  
        {
          Console.WriteLine("Mass extinction at t = " +
            t.ToString().PadLeft(6));
          for (int i = 0; i "lt" numRoaches; ++i)
            herd[i] = new Roach(dim, minX, maxX, tHunger, i);
        }

        ++t;
      } // main while loop

      return bestPosition;

    } // SolveRastrigin

    // ------------------------------------------------------

    public static double Error(double[] x)
    {
      // Rastrigin function has min value of 0.0
      double trueMin = 0.0;
      double rastrigin = 0.0;
      for (int i = 0; i "lt" x.Length; ++i)
      {
        double xi = x[i];
        rastrigin += (xi * xi) -
          (10 * Math.Cos(2 * Math.PI * xi)) + 10;
      }

      return (rastrigin - trueMin) * (rastrigin - trueMin);
    }

    // ------------------------------------------------------

    static double Distance(double[] pos1, double[] pos2)
    {
      double sum = 0.0;
      for (int i = 0; i "lt" pos1.Length; ++i)
        sum += (pos1[i] - pos2[i]) * (pos1[i] - pos2[i]);
      return Math.Sqrt(sum);
    }

    // ------------------------------------------------------

    static void Shuffle(int[] indices, Random rnd)
    {
      //Random rnd = new Random(seed);
      for (int i = 0; i "lt" indices.Length; ++i)
      {
        int r = rnd.Next(i, indices.Length);
        int tmp = indices[r];
        indices[r] = indices[i];
        indices[i] = tmp;
      }
    }

  } // Program

  // --------------------------------------------------------

  public class Roach
  {
    public int dim;
    public double[] position;
    public double[] velocity;
    public double error; // at curr position

    public double[] personalBestPos; // best position found
    public double personalBestErr; // 

    public double[] groupBestPos; // by any roach
 
    public int hunger;
    private Random rnd;

    // ------------------------------------------------------

    public Roach(int dim, double minX, double maxX,
      int tHunger, int rndSeed)
    {
      this.dim = dim;
      this.position = new double[dim];
      this.velocity = new double[dim];
      this.personalBestPos = new double[dim];
      this.groupBestPos = new double[dim];

      this.rnd = new Random(rndSeed);
      this.hunger = this.rnd.Next(0, tHunger);

      for (int i = 0; i "lt" dim; ++i)
      {
        this.position[i] = 
          (maxX - minX) * rnd.NextDouble() + minX;
        this.velocity[i] = 
          (maxX - minX) * rnd.NextDouble() + minX;
        this.personalBestPos[i] = this.position[i];
        this.groupBestPos[i] = this.position[i];
      }

      this.error = 
        RoachOptimizationProgram.Error(this.position);
      this.personalBestErr = this.error;
    }

  } // class Roach

} // ns