Vectorizing Loops In Matlab

Question

I am trying to find a smooth curve through data points using Beizer Curves.

I have the following loop that computes the Binomial Coefficients for the Beizer Curves:

for i = 1:n; // Here 'n' is total number of Control Points.
    a(i) = (factorial(n-1))/(factorial(i-1)*(factorial(n-i)));
end

How would you 'Vectorize' this code in Matlab ?

If I do the following ... I am getting an error.

i = 1:n;
a(i) = (factorial(n-1))/(factorial(i-1)*(factorial(n-i)));

Another question I have is ... if I have nested loops how would you 'Vectorize' that ?

For example: To find all the other points for the Beizer Curve ... I have the following code:

for j=2:c-1
x=0;
y=0;
for i=1:n
    kx(i)=a(i)*t(j).^(i-1)*(1-t(j)).^(n-i)*px(i);
    ky(i)=a(i)*t(j).^(i-1)*(1-t(j)).^(n-i)*py(i);
    x=x+kx(i);
    y=y+ky(i);
end
bx(j)=x;
by(j)=y;
end

Answer 1

For your first example, there are two issues. The first, as pointed out by @Cris Luengo, is that you need to use element-wise division (and also element-wise multiplication). The other semi-issue is that you do not need to use the syntax a(i), which will try to assign the entire output array into a subset of a (using i as the indices). In this case that will be fine, but it's unnecessary and in some other cases may not do what you want/expect; it's a habit you need to be careful about. This formulation produced the same output as your original form, in my limited testing:

i = 1:n;
a2 = (factorial(n-1))./(factorial(i-1).*(factorial(n-i)));

I should note that there is a built-in function nchoosek() to do essentially this calculation, and it supports a vector input, but only for the "n" term (which does not vary in your case). So if you want your code to be vectorized, it may not be useful to you. Although since I assume faster runtime is your goal, it might be worth testing a non-vectorized version that calls nchoosek(), just in case there is some optimization voodoo at work there. You can measure runtime with various tools, such as tic()/toc() or the profiler.

For your second example, depending on the specifics of your problem it may be easy or difficult. Your example relies on enough pre-defined variables (like c, t, etc.) that it was less trivial for me to test. But in general the strategy would be to define one of the independent variables as a row vector, and the other as a column vector. For example, say that i is of size n x 1 and j is of size 1 x m. Variables that vary only with respect to i would also be size n x 1; those that vary only with respect to j would be 1 x m; and those that vary with respect to both would be n x m. You then need to use only functions and operators that can handle the arrays that you feed them, and you likely need to not use the elementwise operators that you needed in the vectorized example above. Here is a somewhat trivial example:

n = 6;
m = 4;

i = (1:n)'; % size 6 x 1
j = 1:m; % size 1 x 4

q = i * j; % size 6 x 4

All of the basic math operations (multiplication, addition, etc.) support vectors (and often n-dimensional arrays as well), and they will do the "right" linear algebra thing with them. Additionally, depending on your version of Matlab, some of the basic operations also support implicit array expansion (so for example you can add a constant to a vector, and it will just "do what you mean", even though this is improper linear algebra)

For more complex situations, this can get a bit unwieldy to understand. Creating temporary intermediate variables can help, so that your code is doing fewer things at a time.

Lastly as a parting shot, I should note that Mathworks has historically discouraged the use of i and j as variable names (this could have changed in the newer versions). The reason is that Matlab uses i and j to indicate imaginary numbers. So if you use them in your code, then Matlab has to figure out what you mean. This can make your code slower, or in the worst case it can result in the wrong behavior. A common practice is to use ii and jj instead, or just pick different iterator variables entirely.