More on expansion: arrayfun

Since there seems to be so much interest in the topic of going beyond scalar expansion, it seems worth talking about yet some other ways to achieve the same goals. I am very appreciative of the thoughtful comments we've received so far on this topic and I know we will be weighing this input carefully as we design how to move forward.

Suppose we want to raise all the values in x to all the powers in powers, where both x and powers are vectors, like this:

powers = 1:3;
x = (0:.1:.4)';

We know from the last blog article, we can achieve this by creating larger intermediate arrays like this, using Tony's trick (used for most of the work in repmat:

X = x(:,ones(size(powers,2),1))
P = powers(ones(size(x,1),1),:)

X =

         0         0         0
    0.1000    0.1000    0.1000
    0.2000    0.2000    0.2000
    0.3000    0.3000    0.3000
    0.4000    0.4000    0.4000

P =

     1     2     3
     1     2     3
     1     2     3
     1     2     3
     1     2     3

and then we are ready to perform the desired calculation:

q = X.^P

q =

         0         0         0
    0.1000    0.0100    0.0010
    0.2000    0.0400    0.0080
    0.3000    0.0900    0.0270
    0.4000    0.1600    0.0640

but at the cost of 2 larger arrays, X and P, that we might not even need going forward.

In the absence of new operators in MATLAB 7, how might we solve this today without creating these temporary arrays and without using a for loop? A new feature in R14sp3, arrayfun is a good candidate. Here are two ways to achieve the same results using arrayfun, the first iterating calculating results for each value of x and the second, for each value of powers. I am taking advantage of anonymous functions, new in MATLAB 7, to define the function I want to calculate, first as a function of powers, and second, as a function of the x values.

pq = cell2mat(arrayfun(@(p)x.^p,powers,'uniformoutput',false));
xq = cell2mat(arrayfun(@(x)x.^powers,x,'uniformoutput',false));

Before proceeding, let's just make sure we're getting the right answer!

isequal(xq,pq,q)
ans =

     1

Now that we've taken care of checking for correctness, let's decipher what the statements mean. Taking the final arguments to arrayfun first, I chose nonuniform output since for each value of my function, I calculate a vector of output values, not a single value. This places outputs for each input value (in this case, the values in powers), into cells in a cell array.

The first argument is an anonyous function, a function defined "inline" that takes its non-argument values from the workspace. So the first function, @(p)x.^p, takes x from the workspace. The function itself takes the values in x and raises them each to some value p, the function argument. Then with the second input to arrayfun, we set the argument to be iterated over to powers from our workspace.

Finally, since our output is now in a cell array, in this case with each cell containing equal numbers of elements, we convert the cell array to a numeric matrix using cell2mat.

Which of the 2 methods using arrayfun is preferable? Well, as is so often, that depends. Since there's an intermediate cell array, we might want to minimize its size. How large are typical sizes for the vectors x and powers? If one is typically substantially larger than the other, we might want to profile the code, either from the command line or the desktop. The profiler can help us make these tradeoffs.

Once again I ask for you to add your thoughts here. Thanks.