When migrating financial or scientific simulations from MATLAB to C++, developers often expect an immediate 10x to 100x performance leap. However, upon rewriting vectorized code like randn(N, M) and element-wise array expressions into nested C++ loops, the results are frequently disappointing: verbose syntax and execution speeds that can actually lag behind MATLAB.

Understanding why this happens—and how to apply modern C++ paradigms—allows you to achieve both expressive, maintainable code and performance that easily outpaces MATLAB.

Why Naive C++ Often Lags Behind MATLAB Vectorization

MATLAB is fundamentally backed by highly optimized native libraries, primarily Intel MKL (Math Kernel Library) and Intel VSL (Vector Statistics Library). When you write:

X = randn(N, M);
Y = exp(mu + sigma .* X);
payoff = max(Y - K, 0);
price = mean(payoff);

MATLAB executes this with:

  • Vectorized pseudo-random generation: Using batch generation algorithms (like SIMD-optimized Ziggurat or Box-Muller) tuned for your CPU architecture.
  • Multi-threading: Array operations in MATLAB automatically dispatch across multiple CPU cores via OpenMP.
  • Vectorized transcendental functions: exp() is computed across entire cache lines using AVX-512 or AVX2 vector registers rather than one scalar double at a time.

Conversely, a naive C++ nested loop using standard <random> (std::mt19937 and std::normal_distribution) and std::exp evaluates scalar operations sequentially on a single thread. The scalar overhead of std::exp and random distribution sampling quickly dominates execution time.

Solution 1: Streamed Parallel Execution with OpenMP (Maximum Cache Efficiency)

While MATLAB allocates a massive N x M matrix in memory, C++ can compute this in-place or via streaming, vastly reducing RAM bandwidth and keeping data in L1/L2 cache. Using OpenMP and a fast random number generator (such as PCG or Xoshiro), you can match or exceed MATLAB's multithreading easily:

#include <iostream>
#include <vector>
#include <cmath>
#include <random>
#include <omp.h>

double monte_carlo_option_price(std::size_t N, std::size_t M, double S0, double K, double r, double sigma, double T) {
    const double dt = T / M;
    const double drift = (r - 0.5 * sigma * sigma) * dt;
    const double vol = sigma * std::sqrt(dt);
    
    double global_payoff_sum = 0.0;

    #pragma omp parallel reduction(+:global_payoff_sum)
    {
        // Thread-local RNG to avoid contention
        std::mt19937_64 rng(1337 + omp_get_thread_num());
        std::normal_distribution<double> norm(0.0, 1.0);

        #pragma omp for schedule(static)
        for (std::size_t i = 0; i < N; ++i) {
            double spot = S0;
            for (std::size_t j = 0; j < M; ++j) {
                spot *= std::exp(drift + vol * norm(rng));
            }
            global_payoff_sum += std::max(spot - K, 0.0);
        }
    }

    return std::exp(-r * T) * (global_payoff_sum / N);
}

Why this is better: It removes allocation overhead entirely and achieves linear scaling across physical CPU cores.

Solution 2: Expression Templates via Eigen (MATLAB-Like Syntax)

If you prefer code that mirrors MATLAB's array syntax without manual loops, the Eigen library is the standard choice. Eigen uses C++ expression templates to fuse array operations, avoiding intermediate allocations:

#include <Eigen/Dense>
#include <iostream>
#include <random>

double run_simulation_eigen(int N, int M, double mu, double sigma, double K) {
    std::mt19937_64 gen(42);
    std::normal_distribution<double> dist(0.0, 1.0);

    // Allocate matrix
    Eigen::ArrayXXd X(N, M);
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < M; ++j) {
            X(i, j) = dist(gen);
        }
    }

    // Vectorized arithmetic matching MATLAB syntax
    Eigen::ArrayXXd Y = (mu + sigma * X).exp();
    Eigen::ArrayXXd payoff = (Y - K).max(0.0);
    
    return payoff.mean();
}

Eigen automatically enables SIMD vectorization (SSE, AVX, NEON) for arithmetic and unary math functions when compiled with -O3 -march=native.

Solution 3: The True MATLAB Equivalent (Intel oneMKL)

To replicate MATLAB’s exact internal mechanism, use Intel oneMKL’s VSL (Vector Statistics Library) and VML (Vector Math Library). This provides batch RNG and vector math written in hand-tuned assembly:

#include <mkl.h>
#include <vector>
#include <algorithm>

double run_mkl_simulation(std::size_t total_samples, double mu, double sigma, double K) {
    std::vector<double> X(total_samples);
    std::vector<double> Y(total_samples);

    // 1. Vectorized RNG
    VSLStreamStatePtr stream;
    vslNewStream(&stream, VSL_BRNG_SFMT19937, 42);
    vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, stream, total_samples, X.data(), mu, sigma);
    vslDeleteStream(&stream);

    // 2. Vectorized exp()
    vdExp(total_samples, X.data(), Y.data());

    // 3. Reduction / Payoff
    double sum = 0.0;
    #pragma omp parallel for reduction(+:sum)
    for (std::size_t i = 0; i < total_samples; ++i) {
        sum += std::max(Y[i] - K, 0.0);
    }

    return sum / total_samples;
}

Can C++20 Ranges Replace Vectorization Libraries?

C++20 ranges provide expressive, lazy evaluation pipelines:

auto payoffs = views::iota(0, N)
             | views::transform([&](int) { return compute_path(); });

While ranges enhance readability, standard range pipelines currently do not auto-vectorize transcendental functions (SIMD) as effectively as Eigen or oneMKL. They are best suited for structural clarity rather than raw vector floating-point throughput.

Summary Recommendations

  • For maximum raw speed: Use streamed OpenMP execution with thread-local lightweight RNGs (such as Xoshiro256**) to avoid allocating large arrays altogether.
  • For code readability matching MATLAB: Use Eigen with compiler flags -O3 -march=native.
  • For industrial-grade quantitative finance: Use Intel oneMKL to leverage SIMD-optimized vector random generation and transcendentals.