'ArrayValued' Integral Extremely Slow: Slower than For-Loop Through Each Element of Array

The title says it all. I first wrote out some code for creating a vector of simulated data through integrating a function for many different values of Q:

%%For-loop implementation

S_ani_int_func = @(Q,r)(r.^3 .* gam(r) .* diff_g_r_func(r) .* sphbesselj(2,Q*r));
S_ani = @(Q)(C * integral(@(r)S_ani_tanh_int_func(Q,r),0,Inf));
my_S_ani_Q = zeros(Q_len,1);
for Q_i = 1:Q_len
    my_S_ani_Q(Q_i) = S_ani(Q_vec(Q_i));
end

I then realized the integral function had a built-in array-valued flag, which I then tried to use to my advantage:

%%'ArrayValued' implementation

S_ani_int_func = @(r)(r.^3 .* gam(r) .* diff_g_r_func(r) .* sphbesselj(2,r.*Q_vec));
my_S_ani_Q = C * integral(S_ani_int_func,0,Inf,'ArrayValued',true);

The two implementations return the exact same vectors, which is good, but the stock array-valued way is for some reason A LOT slower. I run this piece of code in a larger loop several hundred times, and it becomes very obvious: the latter takes at least three times as long. Why is this the case?

[Chad O'Brien]

I observe the same thing for small length vectors. However if you increase the number of points in the array you'll see it switches. I'd like to know why as well. There must be some non-scalable function integral calls when it is in arrayvalued mode.

Answer 1

[Mike Hosea]

Normally the integral function makes vectorized calls over the interval of integration, approximating the subintegrals over each active subinterval simultaneously. You'll see that the integrand f(x) is called with a vector of x values that range across the interval of integration. The length of the vector is typically in the hundreds or thousands.

With ArrayValued=true, the subintervals are handled serially while vectorization of the arithmetic is performed over the elements of the vector-valued function. You'll see that the integrand function f(x) is called only with scalar values of x. The vectors returned by f(x) need to be long enough before the vectorization of the arithmetic using these vectors represents a significant time savings.

Also, different components of an array-valued problem may require smaller subintervals in different places along the interval of integration. Where one component could have a long subinterval, another may require a short one. Solving the problem as an array-valued one implies that the subintervals will generally be shorter, since they must be sized to meet the error requirement implied by the most restrictive component at every point along the interval of integration. This is one reason why we haven't added an ArrayValued option for integral2 or integral3. The vectorization issue is also in play, as the crossover would likely be moved even further out in 2D and 3D problems, making it usually better just to call these functions in a loop.