Integrate function that calls a function

I'm having trouble integrating a function that contains another function. To give an example:

   value = X*func(y)
   prob = lognpdf(y, muy, stdy)
   int_func = @(y)(value(y)*prob(y))

i.e. I have a function value that depends on a variable y, which is log-normally distributed. The problem is that the multiplication in the first line does not work when done inside the integral function:

>> integral(int_func, 0, 5)
Error using _*_: Innder matrix dimensions must agree

As func(y) creates a matrix, the size of which does not line up with the size of X when it is not passed a scalar as input. It appears that integral() internally passes vectors to the function that is being integrated, which then causes the error. Is there any way around this?

Answer 1

[Mike Hosea]

Except when computing array-valued integrals, the INTEGRAL function requires that the integrand function be coded to accept arrays as input and return an array of the same size as output. We call this "vectorization". It is critical to performance in MATLAB. One rarely needs matrix multiplication * when defining an integrand function, rather element-wise multiplication .* .

If I understand you correctly, this is one of those cases where an intermediate quantity is an array, and yet the final result of evaluating the integrand function with a scalar is a scalar. We just need to understand that INTEGRAL is going to pass you an array, and it expects you to evaluate your integrand function at each element of the array. Given that, and the simplicity of the problem definition, we can go one of three ways with this.

One is to use 'ArrayValued',true. This is for calculating the integral of array-valued functions, usually vector-valued functions. Of course a scalar is a vector of length 1. It's the same computation, but in this mode INTEGRAND just passes scalars to the integrand function because it is expecting vectors (or any other fixed shape) back. So using 'ArrayValued',true with a scalar-valued problem essentially has the effect of turning off "vectorization" inside of the integral function. As a result, it will be slower than if we fixed the integrand to accept and return arrays, but the results should be valid.

Two is to use arrayfun. In this approach, you define your integrand f(y) however is natural for you, but you don't integrate f(y) directly, rather you pass this

fv = @(y)arrayfun(f,y);

to INTEGRAL. That is the syntax when f is a function_handle. If, on the other hand, you had an m-file named myfun.m, you would write

fv = @(y)arrayfun(@myfun,y);

These ARRAYFUN definitions mean "Take this function (the first input) and evaluate it element-by-element on entries of the following array (the second input).". It's like writing a loop. This is usually faster than 'ArrayValued',true on scalar-valued integrals.

Finally, in this case you might be able to write func in such a way that it would accept a row (or column) vector input and return an array of size 128-by-numel(y). Then X*func(y) would return a row-vector, which we would then reshape back to size y, i.e.

value = @(y)reshape(X*func(y(:).'),size(y))
prob = @(y)lognpdf(y, muy, stdy)
int_func = @(y)(value(y).*prob(y))

Note the use of elementwise multiplication .* rather than * in the definition of int_fun.

[Nils]

Hi Mike,

thanks for the detailed answer, which gives me a much better understanding of what's going on. One more thing I'd like to ask, if you could shed any light on this: I'll have to perform this integration a couple million times, what is the performance loss I'll likely experience when working without the built-in vectorization?

Thanks!

[Mike Hosea]

Varies, but let me check real quick on something like this...

[Mike Hosea]

I did a little example and got 18x speedup by rewriting func to accept a vector input and return a 128-by-numel(y) array. I don't know what's typical, but I'm not surprised by that number. It's significant.

Answer 2

[Honglei Chen]

What is the dimensions of value and prob? What are the operation you want to do?

If the function is vector valued, integral can deal with it by setting 'ArrayValued' to true

http://www.mathworks.com/help/matlab/ref/integral.html

[Nils]

Hi Honglei,

thanks for your answer. Using the ArrayValued argument does give me an answer, although I'm not entirely sure it's the correct one (I don't have any analytical expression for the function I'm integrating, so verifying is not straightforward). Can you give an intuitive explanation of what ArrayValued does in this case?

Before I answer your question about the dimensions of my objects, let me note that I have made a mistake in my original post. The second line should have read:

prob = @(y)(lognpdf(y, muY, stdY))

so that both value and prob are functions of y. The dimensions are as follows:

>> size(X)
1 128
>> size(func(y))  % with y scalar
128 1
>> size(value)
1 1
>> size prob(y)  % with y scalar
1 1

so that in the case of scalar y I simply evaluate value at y, which gives a scalar, and multiply by the probability of y occuring, which is also scalar. However, when I integrate using integral, it seems that a vector is passed to the func function in value, so that the multiplication X*func(y) fails. This is as the size of X is fixed, but func(y) changes shape when a vector is passed as an argument rather than a scalar.