I wrote a program which iteratively transforms data using matrix multiplications. To minimize the number of large memory allocations, I use two roughly equal-sized std::vector<double>'s data and temp_data to store the current data and the new data respectively. For the matrix operations, I wrap each buffer into an Eigen::Map<Eigen::MatrixXd>, called data_matrix and temp_data_matrix respectively. Then, I need to compute the following operation:

temp_data_matrix.noalias() = data_matrix.transpose()*factor;

where factor is simply an Eigen::MatrixXd which is much smaller than data_matrix or temp_data_matrix. Thanks to .noalias(), this line does not allocate any large temporaries normally. However, when multi-threading the program suddenly allocates a temporary (roughly) equal to the size of temp_data! In practice this can increase memory usage by around 50% so I'm curious about a few things:

  • Why does this only happen when multi-threading? What algorithmic aspect requires a temporary buffer on top of the output buffer?
  • Why does this not happen with the line temp_data_matrix.noalias() = factor*data_matrix? After all, my operation seems more cache-efficient: because we transpose data_matrix, we simply need to compute the scalar product of each pair of columns from data_matrix and factor, which should be fast considering both matrices are column-major. However, I'm not considering the effects of blocking here, which I know can be very important for performance.
  • What would be the fastest way to implement the aforementioned line without allocating any temporaries of the order of data or temp_data? In case the parallel algorithm needs some "scratch space", can't it use temp_data_matrix for that? I don't care about copying factor because it's relatively small in practice. Right now the fastest solution seems to be to switch back to single-threading for this line, but I hope this can be improved. I tried using lazyProduct but this turned out to be much slower than the single-threaded solution.

Appendix: The observations above were made using Massif (valgrind) both in my original program and with the example code below. I compiled the example program using the command g++ -O3 -DNDEBUG -I /usr/include/eigen3 -fopenmp eigen_memory_test.cpp -o eigen_memory_test.

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

int main(int argc, char** argv) {
    bool parallel = (argc > 1 && strcmp(argv[1], "parallel") == 0);
    std::cout << "parallel: " << (parallel  ? "true" : "false") << std::endl;
    bool lazyproduct = (argc > 2 && strcmp(argv[2], "lazyproduct") == 0);
    std::cout << "lazyproduct: " << (lazyproduct  ? "true" : "false") << std::endl;
    bool test = true;
    // Initialize parallellism
    int threads = 4;
    // Initialization and memory allocation
    size_t columns = 100000;
    std::vector<size_t> rows_per_step = {100, 90, 80, 70, 60, 50};
    size_t steps = rows_per_step.size() - 1;
    std::vector<double> data(rows_per_step[0]*columns);
    std::vector<double> temp_data(rows_per_step[0]*columns);
    // Generate random data
    std::random_device rnd_device;
    std::mt19937 mersenne_engine {rnd_device()};
    std::uniform_real_distribution<double> distribution(-1.0, 1.0);
    auto generator = std::bind(distribution, std::ref(mersenne_engine));
    std::generate(data.begin(), data.end(), generator);
    // If testing, print first elements to verify later
    if (test)
        std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;
    // Process steps
    for (size_t step = 0; step < steps; step++) {
        size_t initial_rows = rows_per_step[step];
        size_t new_rows = rows_per_step[step + 1];
        // Choose factor
        Eigen::MatrixXd factor = (test ? 
            Eigen::MatrixXd(Eigen::MatrixXd::Identity(initial_rows, new_rows)) : 
            Eigen::MatrixXd(Eigen::MatrixXd::Random(initial_rows, new_rows))
        // Transform data
        Eigen::Map<Eigen::MatrixXd> data_matrix(data.data(), initial_rows, columns);
        Eigen::Map<Eigen::MatrixXd> temp_data_matrix(temp_data.data(), columns, new_rows);
        if (!parallel)
        // In my real problem setting, the transpose in the following matrix multiplication is actually necessary
        if (lazyproduct)
            temp_data_matrix.noalias() = data_matrix.transpose().lazyProduct(factor);
            temp_data_matrix.noalias() = data_matrix.transpose()*factor;
        // Transpose back to original data
        // This doesn't happen in the real problem, but I put this here to keep this example simpler
        Eigen::Map<Eigen::MatrixXd> transposed_data_matrix(data.data(), new_rows, columns);
        transposed_data_matrix.noalias() = temp_data_matrix.transpose();
    // Print first elements to make sure matrix operations aren't optimized away
    std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;
    return 0;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.