0

I am filling an Eigen matrix with the following code:

int M = 3;
int N = 4;
MatrixXd A(M, N);

double res = sin(4);

for (int i = 0; i < M; i++) {
    for (int j = 0; j < N; j++) {
        A(i, j) = sin(i+j);
    }
}

In Matlab I only need 1 for loop to do the same thing using vectorization:

M = 3;
N = 4;
N_Vec = 0:(N-1);
A = zeros(M,N);
for i=1:M
    A(i,:) = sin((i-1)+N_Vec);
end

Is it possible to do something similar in C++/Eigen so that I can get rid of one of the for loops? If it is possible to somehow get rid of both for loops that would be even better. Is that possible?

csss
  • 1,937
  • 2
  • 14
  • 24
  • 1
    Related: https://stackoverflow.com/questions/53877801/eigen-indexing-update-all-column-of-a-specific-row – πάντα ῥεῖ Feb 07 '21 at 13:13
  • @πάνταῥεῖ I should mention that I used M=3 and N=4 just for simplicity in my question. In M and N will be hundreds of times larger. So in post you linked the explicit setting of vectors (`vec << i, i + 1, i + 2, i+3, i+4;`) will not be possible as it is not feasible to type `vec << i, i+1, i+2, ...` all the way up to say N=500. – csss Feb 07 '21 at 13:32
  • Well, if you have to set a large number of values individually using a calculation, there's no other way than looping through them. At least none I know. – πάντα ῥεῖ Feb 07 '21 at 13:34
  • Is that not extremely inefficient? Matlab will surely be much, much faster than Eigen if we can't do vectorization in Eigen. https://uk.mathworks.com/help/matlab/matlab_prog/vectorization.html – csss Feb 07 '21 at 13:36
  • No less efficient than specifying a function to be applied to all values of a row. If Eigen doesn't provide such functionality, it's easy to encapsulate the inner loop into another function which takes a row vector, and a function pointer applicable to all individual values, so that you can write a single call to that function instead of the inner loop. In other words: The looping has to be done somewhere, where it's done doesn't have much impact on _"efficiency"_, other than readability of the code. – πάντα ῥεῖ Feb 07 '21 at 13:41
  • 2
    MATLAB can do this without any loops at all. But just because the MATLAB code doesn’t show loops, doesn’t mean they’re not there. MATLAB runs those loops for you behind the scenes. Doing a loop in C++ is not slow (just like it has no longer been slow in MATLAB for many years now, the slowness of loops is something from 15 years ago). – Cris Luengo Feb 07 '21 at 14:52
  • 1
    "Vectorization" = "making a faster programming language handle the loops for you" – beaker Feb 07 '21 at 16:40
  • @CrisLuengo If I run the code in Matlab with two for loops with N,M = 2e3 it takes about 0.086 seconds, whereas the code with one for loop and vectorization takes 0.037 seconds. Maybe for loops are not as slow as they were 15 years ago but they are certainly still slower than vectorization, hence why MathWorks say 'Vectorized code often runs much faster than the corresponding code containing loops' here: https://uk.mathworks.com/help/matlab/matlab_prog/vectorization.html – csss Feb 10 '21 at 11:20
  • 2
    That used to be a 100x difference back in the day. Nowadays there are many cases of vectorized code being slower than the equivalent loop. Your case has just a function call as loop body. Function calls still have a relatively high overhead in MATLAB, so calling it once or many times makes a difference. Anyway, calling `sin` with an array causes the execution of a loop written in C or C++ under the hood. – Cris Luengo Feb 10 '21 at 14:49

1 Answers1

5

Using a NullaryExpr you can do this with zero (manual) loops in Eigen:

Eigen::MatrixXd A = Eigen::MatrixXd::NullaryExpr(M, N,
      [](Eigen::Index i, Eigen::Index j) {return std::sin(i+j);});

When compiled with optimization this is not necessarily faster than the manual two-loop version (and without optimization it could even be slower).

You can write int or long instead of Eigen::Index, if that is more readable ...

chtz
  • 17,329
  • 4
  • 26
  • 56