Mapslices and cumsum


#1

I’m finding that the following code grinds my machine to a halt:

n=10^6;
X = randn(1,n);
function V(x)
    return 0.5 * dot(x,x);
end
avgX = cumsum(mapslices(V,X,1))./(1:n);

And by halt, I mean, I had to do a hard reset on one computer.

In contrast, the following code works just fine:

n=10^6;
X = [randn(1) for i=1:n]
function V(x)
    return 0.5 * dot(x,x);
end
avgX = cumsum(V.(X))./(1:n);

#2

The problem is that cumsum(mapslices(V,X,1)) generates a row vector, while (1:n) is interpreted as a column vector. The ./ generates a 1e6 x 1e6 dense matrix. What you want is

avgX = cumsum(mapslices(V, X, 1))./((1:n)')

HOWEVER, this is still much slower than your second version, which is counter-intuitive for me. Maybe somebody else can explain.


#3

Shameless advertisement: You can use UnsafeArrays.@uviews instead of mapslices:

using UnsafeArrays
X = randn(1,n)
@uviews X cumsum([V(view(X, :, i)) for i in indices(X, 2)])

On my system, this is as fast as the @gideonsimpson’s second version, while using the
X-matrix of the first version.