determining when calc_vec can be used

[paciorek]

We need to determine from a given declaration if use of calc_vec, which will rely on Eigen's vectorized expressions, is valid.

Here are my somewhat tentative thoughts.

First, I'm not seeing a lot of documentation regarding how Eigen handles vectorized expressions such as y[1:5] = exp(x[1:5]). AI suggests that it sets up a for loop that avoids creating any temporary arrays and also relies on the processors vectorized/SIMD instructions (e.g., AVX) to compute multiple values in parallel. And that for best efficiency values should be contiguous in memory (but I think this would generally hold for calc_one too and have more to do with the order in which the loops are executed).

So I think we want to consider the following:

• When there are multiple loops, choose an indexing variable to vectorize over. I suppose consider the index associated with the loop with the most iterations as the focal index and vectorize only with respect to that index (though I'm not clear on the mechanism by which for(i in 1:8) exp(x[1:000,i]) is better than for (i in 1:1000) exp(x[i,1:8])).
• Check that all indexing involves scalars to avoid things like x[i,2:5] or x[i:i+3], which presumably couldn't make use of vectorized instructions. That said, what would Eigen do with something like y[1:3, 5:10] = exp(x[1:3,5:10]). Perhaps it would run fine.
• Regarding nested indexing, it's not clear to me that this poses a problem, though it may not provide for any additional efficiency than calc_one.
• Similarly, it's not clear there is a problem with something like y[i+j], if we only have one index vectorized and the other index is looped over explicitly.

I guess one question in all this is if our determination is trying to determine when vectorized code will work versus when it will be more efficient than calc_one.

[perrydv]

Hi @paciorek , thanks for starting this.

• Have you looked at test-math.R? That somewhat documents what operations we support (through Eigen) for vectorization, although it may be a bit thick to get oriented on the test specification objects.
• Be sure when reading on Eigen that you are asking about Eigen's Tensor module, as that is fairly different than standard Eigen, which only handles matrix and array types.
• In some (many element-wise?) cases I think we can vectorize across multiple dimensions in one expression, so exp(x[1:1000, 1:8]).
• In theory I think x[i:(i+3)] could be vectorized.
• Nested indexing: I think x[k[i] : l[i] ] should be fine but x[ k[ i:j] ] would not (easily).
• I'm not clear about e.g. x[(1:10) + j] in terms of the low-level steps and whether it is easy or worthwhile.

[paciorek]

@perrydv and I agreed to punt for now on the question of determining if vectorization is possible/desirable.

Ideally in setting up the calculate stuff, we will allow use of calc_vec and that will give us the ability to easily do some benchmarking and see how much it helps.

I did do a bit of benchmarking of manual looping in nFunction code versus vectorization for exp(x), sin(x), x+4, dnorm(x,0,1). For the first three, using vectorized code can be about 2x as fast. For dnorm it's rather slower (Linux) or the same (MacOS). I'll note that sin gets done via unaryExpr and yet was still faster with vectorization.