multi-threaded random search optimizer

random_search.cpp

// Copyright 2023 Taylor Howell

#include <iostream>
#include <random>
#include <thread>
#include <vector>

// derivative-free random search
template <typename T>
class RandomSearch {
 public:
  // constructor
  RandomSearch(float (*objective)(const std::vector<T>&, int), int dim,
               int n_sample, int n_iter, float stddev) {
    // set attributes
    objective_ = objective;
    dim_ = dim;
    n_sample_ = n_sample;
    n_iter_ = n_iter;

    // allocate memory
    for (int i = 0; i <= n_sample_; i++) {
      parameters_.push_back(std::vector<T>(dim_));
    }
    values_.resize(n_sample_ + 1);

    // normal distribution standard deviation
    stddev_ = stddev;

    // history
    history_.resize(n_iter_);
  }

  // sample parameters
  void SampleParameters(int i) {
    // copy nominal parameters
    parameters_[i] = parameters_[n_sample_];

    // create distribution
    std::random_device rd;
    std::mt19937 gen(rd());
    std::normal_distribution<float> dist(0.0, stddev_);

    // add Gaussian noise
    for (int j = 0; j < dim_; j++) {
      // sample noise
      float noise = dist(gen);

      // add to nominal parameters
      parameters_[i][j] += noise;
    }
  }

  // single random search iteration
  void Iteration() {
    // sample parameters and evaluate objective
    // for (int i = 0; i < n_sample_; i++) {
    //     // noisy parameters
    //     SampleParameters(i);

    //     // objective value for noisy parameters
    //     values_[i] = objective_(parameters_[i], dim_);
    // }

    // sample and evaluate parameters in parallel
    for (int i = 0; i < n_sample_; i++) {
      workers_.push_back(std::thread([&rs = this, i]() {
        // noisy parameters
        rs->SampleParameters(i);

        // objective values
        rs->values_[i] = rs->objective_(rs->parameters_[i], rs->dim_);
      }));
    }

    // wait for all workers
    std::for_each(workers_.begin(), workers_.end(),
                  [](std::thread& t) { t.join(); });

    // reset workers
    workers_.clear();

    // find minimum value
    int idx = std::distance(values_.begin(),
                            std::min_element(values_.begin(), values_.end()));

    // set nominal parameters
    if (idx < n_sample_) {
      parameters_[n_sample_] = parameters_[idx];
      values_[n_sample_] = values_[idx];
    }
  }

  // solve via random search
  void Solve() {
    for (int i = 0; i < n_iter_; i++) {
      // single random search iteration
      Iteration();

      // cache value for history
      history_[i] = values_[n_sample_];
    }
  }

  // initialize parameters
  void Initialize(std::vector<T> parameters) {
    parameters_[n_sample_] = parameters;
    values_[n_sample_] = objective_(parameters, dim_);
  }

  // solve status
  void Status() {
    std::cout << "HISTORY" << std::endl;
    for (int i = 0; i < n_iter_; i++) {
      std::cout << "  value (" << i << "): " << history_[i] << std::endl;
    }

    std::cout << "SOLUTION" << std::endl;

    std::cout << "parameters: " << std::endl;
    for (int i = 0; i < dim_; i++) {
      std::cout << "  [" << i << "] = " << parameters_[n_sample_][i]
                << std::endl;
    }

    std::cout << "value = " << values_[n_sample_] << std::endl;
  }

 private:
  // objective function
  float (*objective_)(const std::vector<T>&, int);

  // parameter dimension
  int dim_;

  // number of random samples
  int n_sample_;

  // number of iterations
  int n_iter_;

  // parameters
  std::vector<std::vector<T>> parameters_;

  // objective values
  std::vector<T> values_;

  // normal distribution standard deviation
  T stddev_;

  // history of best iteration values
  std::vector<T> history_;

  // workers
  std::vector<std::thread> workers_;
};

// quadratic objective
float quadratic(const std::vector<float>& x, int dim) {
  float obj = 0.0f;

  for (int i = 0; i < dim; i++) {
    obj += x[i] * x[i];
  }

  return 0.5 * obj;
}

// rosenbrock
float rosenbrock(const std::vector<float>& x, int dim) {
  float a = 1.0f;
  float b = 100.0f;

  float obj = 0.0;

  float t0 = (a - x[0]);
  float t1 = (x[1] - x[0] * x[0]);

  obj += t0 * t0;
  obj += b * t1 * t1;

  return obj;
}

int main() {
  std::cout << "random search" << std::endl;

  // // constructor
  // RandomSearch optimizer(quadratic, 2, 8, 1000, 0.01);

  // // initial parameters
  // std::vector<float> x = {1.0, 1.0};

  // constructor
  RandomSearch<float> optimizer(rosenbrock, 2, 32, 100, 0.1);

  // initial parameters
  std::vector<float> x = {-1.0, 2.0};

  // initialize
  optimizer.Initialize(x);

  // solve
  optimizer.Solve();

  // status
  optimizer.Status();

  return 0;
}