C++ Numerical Computing

Table of contents

Background

For accuracy, all the results are the execution time after poisoning the L1 cache before each run.

Also, all possible system soft interrupts were bound to core 0 using taskset on the test machine, and these tests were done on core 1. (Rocky Linux 9.2, AMD 7950X, same configuration as the live trading machine we use for quantitative trading)

Interpolation (linear or non-linear)

Armadillo 1-D Linear Interpolation

void armadillo_1d_test() {
  // arma::vec x = arma::linspace<arma::vec>(0, 3, 20);
  arma::vec x = {0, 1, 2, 3, 4, 5};
  arma::vec y =
      arma::square(x);  // arma::vec xx = arma::linspace<arma::vec>(0, 3, 100);
  arma::vec xx = {2.5};
  arma::vec yy;
  // arma::interp1(x, y, xx, yy); // use linear interpolation by default
  arma::interp1(x, y, xx, yy,
                "*linear");  // faster than "linear"
  std::cout << "Interpolated 1d value at 2.5: " << yy << std::endl;
}

boost.Math 1-D Linear Interpolation

void boost_1d_test() {
  std::vector<double> x = {0, 1, 2, 3, 4, 5};
  std::vector<double> y = {0, 1, 4, 9, 16, 25};
  boost::math::interpolators::barycentric_rational<double> interpolator(
      x.begin(), x.end(), y.begin());
  std::cout << "Interpolated 1d value at 2.5: " << interpolator(2.5)
            << std::endl;
}

GSL 1-D Linear Interpolation

void gsl_1d_test() {
  double x[] = {0, 1, 2, 3, 4, 5};
  double y[] = {0, 1, 4, 9, 16, 25};
  int n = sizeof(x) / sizeof(x[0]);
  gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, n);
  gsl_interp_init(interp, x, y, n);
  double yi = gsl_interp_eval(interp, x, y, 2.5, NULL);
  std::cout << "Interpolated 1d value at 2.5: " << yi << std::endl;
  gsl_interp_free(interp);
}

Armadillo 2-D Linear Interpolation

void armadillo_2d_test() {
  arma::mat Z = {{0, 1, 4}, {1, 2, 5}, {4, 5, 8}};
  arma::vec X = arma::regspace(0, Z.n_cols - 1);  // X=horizontal spacing
  arma::vec Y = arma::regspace(0, Z.n_rows - 1);  // Y=vertical spacing
  arma::vec XI = arma::vec{1.5};                    // magnify by approx 2
  arma::vec YI = arma::vec{1.5};                    // magnify by approx 3
  arma::vec ZI;
  arma::interp2(X, Y, Z, XI, YI, ZI);  // use linear interpolation by default
  std::cout << "Interpolated value at (" << 1.5 << ", " << 1.5 << "): " << ZI<< std::endl;
}

boost.Math 2-D Linear Interpolation

void boost_2d_test() {
  std::vector<double> v{0, 1, 4, 1, 2, 5, 4, 5, 8};
  using boost::math::interpolators::bilinear_uniform;
  int rows = 3;
  int cols = 3;
  double dx = 1;
  double dy = 1;
  auto bu = bilinear _ uniform(std::move(v), rows, cols, dx, dy);
  // evaluate at a point:
  double x = 1.5, y = 1.5;
  double z = bu(x, y);
  std::cout << "Interpolated value at (" << x << "," << y << "): " << z<< std::endl;
}

GSL 2-D Linear Interpolation

void gsl_2d_test() {
  const size_t xSize = 3;
  const size_t ySize = 3;
  double xArr[xSize] = {0, 1, 2};
  double yArr[ySize] = {0, 1, 2};
  double zArr[xSize * ySize] = {0, 1, 4, 1, 2, 5, 4, 5, 8};
  gsl_interp_accel *xAcc = gsl_interp_accel_alloc();
  gsl_interp_accel *yAcc = gsl_interp_accel_alloc();
  gsl_spline2d *spline =
      gsl_ spline2d_alloc(gsl_interp2d_bilinear, xSize, ySize);
  gsl_spline2d_init(spline, xArr, yArr, zArr, xSize, ySize);
  double x = 1.5, y = 1.5;
  double result = gsl_spline2d_eval(spline, x, y, xAcc, yAcc);
  std::cout << "Interpolated value at (" << x << "," << y << "): " << result<< std::endl;
  gsl_spline2d_free(spline);
  gsl_interp_accel_free(xAcc);
  gsl_interp_accel_free(yAcc);
}

Interpolation Performance Comparisons

Common mathematical functions (such as natural logarithm, standard normal distribution probability density and cumulative density functions, etc.)

Armadillo Common mathematical functions

void armadillo_test() {
  arma::vec x = {1.0};
  std::cout << "ln(" << x << "): " << arma::log(x) << std::endl;
  std::cout << "PDF: " << arma::normpdf(x) << std::endl;
  std::cout << "CDF: " << arma::normcdf(x) << std::endl;
}

boost.MATH Common mathematical functions

void boost_test() {
  double x = 1.0;
  boost::math::normal_distribution<> dist;
  double ln = boost::math::log1p(x);
  double pdf = boost::math::pdf(dist, x);
}

GSL Common mathematical functions

void gsl_test() {
  double x = 1.0;
  double ln = gsl_sf_log(x);
  double pdf = gsl_ran_gaussian_pdf(x, 1.0);
  double cdf = gsl_cdf_gaussian_P(x, 1.0);
}

Common mathematical functions Performance Comparisons

Matrix Operations

Matrix Operations with Eigen

Matrix Add:

void eigen_add_test() {
  Eigen::MatrixXd A(2, 2);
  A << 1, 2, 3, 4;
  Eigen::MatrixXd B(2, 2);
  B << 2, 3, 4, 5;
  Eigen::MatrixXd C = A + B;
}

Matrix Multiplication

void eigen_mul_test() {
  Eigen::MatrixXd A(100, 100);
  Eigen::MatrixXd B(100, 100);
  // Fill matrices A and B with data
  int count = 0;
  for (int i = 0; i < 100; ++i) {
    for (int j = 0; j < 100; ++j) {
      A(i, j) = count;
      B(i, j) = count + 1;  // Just to differentiate from A
      count++;
    }
  }
  Eigen::MatrixXd C = A * B;
}

Matrix Operations with Armadillo

Matrix Add:

void armadillo_add_test() {
  arma::mat A = {{1, 2}, {3, 4}};
  arma::mat B = {{2, 3}, {4, 5}};
  arma::mat C = A + B;
}

Matrix Multiplication

void armadillo_mul_test() {
  arma::mat A(100, 100);
  arma::mat B(100, 100);
  // Fill matrices A and B with data
  int count = 0;
  for (size_t i = 0; i < 100; ++i) {
    for (size_t j = 0; j < 100; ++j) {
      A(i, j) = count;
      B(i, j) = count + 1;  // Just to differentiate from A
      count++;
    }
  }
  arma::mat C = A * B;
}

Matrix Operations with GSL

Matrix Add:

void gsl_add_test() {
  gsl_matrix *A = gsl_matrix_alloc(2, 2);
  gsl_matrix_set(A, 0, 0, 1);
  gsl_matrix_set(A, 0, 1, 2);
  gsl_matrix_set(A, 1, 0, 3);
  gsl_matrix_set(A, 1, 1, 4);
  gsl_matrix *B = gsl_matrix_alloc(2, 2);
  gsl_matrix_set(B, 0, 0, 2);
  gsl_matrix_set(B, 0, 1, 3);
  gsl_matrix_set(B, 1, 0, 4);
  gsl_matrix_set(B, 1, 1, 5);
  gsl_matrix_add(A, B);  // Result is stored in A
  for (size_t i = 0; i < 2; ++i) {
    for (size_t j = 0; j < 2; ++j) {
      std::cout << gsl_matrix_get(A, i, j) << " ";
    }
    std::cout << std::endl;
  }
  gsl_matrix_free(A);
  gsl_matrix_free(B);
}

Matrix Multiplication

void gsl_mul_test() {
  gsl_matrix *A = gsl_matrix_alloc(100, 100);
  gsl_matrix *B = gsl_matrix_alloc(100, 100);
  gsl_matrix *C = gsl_matrix_alloc(100, 100);
  // Fill matrices A and B with data
  int count = 0;
  for (size_t i = 0; i < 100; ++i) {
    for (size_t j = 0; j < 100; ++j) {
      gsl_matrix_set(A, i, j, count);
      gsl_matrix_set(B, i, j, count + 1);  // Just to differentiate from
      A count++;
    }
  }
  gsl_matrix_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, B, 0.0, C);
  // ... use matrix C ...
  gsl_matrix_free(A);
  gsl_matrix_free(B);
  gsl_matrix_free(C);
}

Matrix Operations Comaprisons

Summary

GSL is generally a bit tricky to write, it’s C-style, and you have to take care of your own to be aware of the memory release.

In terms of speed:

Armadillo is faster for 1D linear interpolation, boost is faster for 2D linear interpolation, and then Armadillo is faster than GSL, but there is not much difference between boost and Armadillo are not much different.

For common math functions Armadillo is fastest and GSL is faster than boost.

For matrix addition Armadillo is the fastest, and matrix multiplication Armadillo is also the fastest.

References

• Armadillo: C++ library for linear algebra & scientific computing
• Boost Library Documentation - Math and numerics
• GSL - GNU Scientific Library
• Eigen: a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms