4 ms·

I think you are getting at this above and what I've written below is somewhat redundant with srean's followup, but to clarify, writing vectorized code is mostly

I think you are getting at this above and what I've written below is somewhat redundant with srean's followup, but to clarify, writing vectorized code is mostly orthogonal to using SIMD vector instructions to compute the result. A sufficiently smart compiler can turn a loop into SIMD instructions, and these days most compilers (including LLVM, which Julia uses) are pretty good at that. The problem that Julia solves here is that many scientific programming environments require you to write vectorized code, because they are interpreted and so looping through the elements in a vector incurs massive interpreter overhead. In some cases, this can make your code much more difficult to understand than if you wrote an explicit loop, but even when it does not, it often forces you to give up performance that could be obtained if you could write the code in a devectorized manner. For example, consider the problem of calculating the variance of a sample. To calculate the variance in MATLAB, one might write: mu = mean(x) sum((x - mu).^2) / (length(x) - 1) Although I don't know the details of how MATLAB optimizes this code, I do know that it has roughly equivalent performance characteristics to Julia code that first computes x - mu, then computes abs2(x - mu), then sums along the result. With an interpreter, you can't do much better, and this is actually how the MATLAB var function does it. In Julia, we compute the variance as: mu = mean(x) v = 0. for i = 1:length(x) v += abs2(x[i] - mu) end v / (length(x) - 1) This is several times faster than MATLAB's included var function, because the computations performed on each element are cheap relative to the cost of memory access. If you wrap the loop in the @simd macro and add @inbounds to eliminate the bounds check when accessing x, the LLVM loop vectorizer will actually take that loop and translate it into SIMD instructions, making it even faster. (The @simd macro is necessary to tell LLVM that it's okay to vectorize, since accumulating into n variables and summing them at the end gives different [but usually more accurate] results compared to accumulating into a single variable due to floating point rounding.) The promise of Devectorize.jl is that you will be able to write the first code fragment and have it transformed to the second, which would be neat indeed.

srean 12y ago
Thanks for providing more detail and context. I wonder how the clash in terminology came about ! Who would have imagined naming things uniquely would be so hard (pun intended). I understand that your example is entirely pedagogic, but just a cautionary note for the unwary (a) although it is tempting to fold the variance calculation into a single loop (accumulate the totals of x and x^2), (b) neither that obvious single loop version nor the code above are good ways to compute variance if one cares about preserving precision, more so if x has a wide dynamic range. In large scale problems it does raise its ugly head and these bugs are difficult to catch because the code is mathematically correct. Using double mitigates the problem to an extent (but then floats are faster for SIMD vectorization). Another side note, I have gradually come to realize and appreciate the unique position that Fortran holds. It is not often that you have compiler writers with a background in numerical analysis or vice versa. I BTW have background in none and sorely miss that.
- simonster 12y ago
  The way I compute the variance above is, to my knowledge, the standard algorithm implemented by most software packages (apparently including MATLAB). (The single pass computation you mention is subject to catastrophic cancellation, and thus pretty terrible unless the mean of your data is very small relative to the variance.) However, the "real" implementation in Julia standard library is indeed a bit better: it performs pairwise summation, which has O(log n) error growth instead of O(n) error growth at negligible performance cost (see http://en.wikipedia.org/wiki/Pairwise_summation http://en.wikipedia.org/wiki/Pairwise_summation).
ViralBShah 12y ago
Also, the hope is that in a future version of Julia, we can have a devectorization pass within the compiler at least to handle the most common cases. This is common in Fortran compilers.

srean 12y ago

Thanks for providing more detail and context. I wonder how the clash in terminology came about ! Who would have imagined naming things uniquely would be so hard (pun intended). I understand that your example is entirely pedagogic, but just a cautionary note for the unwary (a) although it is tempting to fold the variance calculation into a single loop (accumulate the totals of x and x^2), (b) neither that obvious single loop version nor the code above are good ways to compute variance if one cares about preserving precision, more so if x has a wide dynamic range. In large scale problems it does raise its ugly head and these bugs are difficult to catch because the code is mathematically correct. Using double mitigates the problem to an extent (but then floats are faster for SIMD vectorization). Another side note, I have gradually come to realize and appreciate the unique position that Fortran holds. It is not often that you have compiler writers with a background in numerical analysis or vice versa. I BTW have background in none and sorely miss that.

simonster 12y ago
The way I compute the variance above is, to my knowledge, the standard algorithm implemented by most software packages (apparently including MATLAB). (The single pass computation you mention is subject to catastrophic cancellation, and thus pretty terrible unless the mean of your data is very small relative to the variance.) However, the "real" implementation in Julia standard library is indeed a bit better: it performs pairwise summation, which has O(log n) error growth instead of O(n) error growth at negligible performance cost (see http://en.wikipedia.org/wiki/Pairwise_summation http://en.wikipedia.org/wiki/Pairwise_summation).

ViralBShah 12y ago

Also, the hope is that in a future version of Julia, we can have a devectorization pass within the compiler at least to handle the most common cases. This is common in Fortran compilers.