0

I'm basicaly trying to find the product of an expression that goes like this:

(x-(N-1)/2).....(x+(N-1)/2) for even value of N

x is a value that I will set at the beginning that changes too but that is a different problem...

let's say for the sake of argument that for now x is a constant (ex x=1)

example for N=6

(x-5/2)(x-3/2)(x-1/2)(x+1/2)(x+3/2)*(x+5/2)

the idea was to create a row vector every element of which is each individual term (P(1)=x-5/2) (P(2)=x-3/2)...etc and then calculate its product

N=6;
x=1;


P=ones(1,N);
for k=(-N-1)/2:(N-1)/2

for n=1:N
            P(n)=(x-k);

end
end
y=prod(P);

instead this creates a vector that takes only the first value of the epxression and then repeats the same value at each cell.

there is obviously a fundamental problem with my loop but I just can't see it.

So if anyone can help with that OR suggest a better way to calculate the product I would be grateful.

0

3 Answers 3

2

Use vectorized commands

Why use a loop when you can use vectorized commands like prod?

y = prod(2 * x + [-N + 1 : 2 : N - 1]) / 2;

For convenience, you may want to define an anonymous function for it:

f = @(N,x) reshape(prod(bsxfun(@plus, 2 * x(:), -N + 1 : 2 : N - 1) / 2, 2), size(x));

Note that the function is compatible with a (row or column) vector input x.

Tests in MATLAB's Command Window

>> f(6, [2,2]')

ans =
  -14.7656
    4.9219
   -3.5156
    4.9219
  -14.7656

>> f(6, [2,2])

ans =
  -14.7656    4.9219   -3.5156    4.9219   -14.7656

Benchmark

Here is a comparison of rayreng's approach versus mine. The former emerges as the clear winner... :'( ...at least as N increases.

Varying N, fixed x

enter image description here

Fixed N (= 10), vector x of varying length

enter image description here

Fixed N (= 100), vector x of varying length

enter image description here

Benchmark code

function benchmark
% varying N, fixed x

    clear all
    n = logspace(2,4,20)';
    x = rand(1000,1);
    tr = zeros(size(n));
    tj = tr;

    for k = 1 : numel(n)
        % rayreng's approach (poly/polyval)
        fr = @() rayreng(n(k), x);
        tr(k) = timeit(fr);

        % Jubobs's approach (prod/reshape/bsxfun)
        fj = @() jubobs(n(k), x);
        tj(k) = timeit(fj);
    end

    figure
    hold on
    plot(n, tr, 'bo')
    plot(n, tj, 'ro')
    hold off
    xlabel('N')
    ylabel('time (s)')
    legend('rayreng', 'jubobs')
end

function y = jubobs(N,x)
    y = reshape(prod(bsxfun(@plus,...
                            2 * x(:),...
                            -N + 1 : 2 : N - 1) / 2,...
                     2),...
                size(x));
end

function y = rayreng(N, x)
    p = poly(linspace(-(N-1)/2, (N-1)/2, N));
    y = polyval(p, x);
end

function benchmark2
% fixed N, varying x    

    clear all
    n = 100;
    nx = round(logspace(2,4,20));
    tr = zeros(size(n));
    tj = tr;

    for k = 1 : numel(nx)
        disp(k)
        x = rand(nx(k), 1);

        % rayreng's approach (poly/polyval)
        fr = @() rayreng(n, x);
        tr(k) = timeit(fr);

        % Jubobs's approach (prod/reshape/bsxfun)
        fj = @() jubobs(n, x);
        tj(k) = timeit(fj);
    end

    figure
    hold on
    plot(nx, tr, 'bo')
    plot(nx, tj, 'ro')
    hold off
    xlabel('number of elements in vector x')
    ylabel('time (s)')
    legend('rayreng', 'jubobs')
    title(['n = ' num2str(n)])
end

function y = jubobs(N,x)
    y = reshape(prod(bsxfun(@plus,...
                            2 * x(:),...
                            -N + 1 : 2 : N - 1) / 2,...
                     2),...
                size(x));
end

function y = rayreng(N, x)
    p = poly(linspace(-(N-1)/2, (N-1)/2, N));
    y = polyval(p, x);
end

An alternative

Alternatively, because the terms in your product form an arithmetic progression (each term is greater than the previous one by 1/2), you can use the formula for the product of an arithmetic progression.

1
  • Wow, I'm surprised that poly and polyval beat out bsxfun.... that's very surprising indeed. Either way, thank you for the benchmarking :)
    – rayryeng
    Dec 14, 2014 at 3:04
1

I agree with @Jubobs in that you should avoid using for loops for this kind of computation. There are cases where for loops perform fast, but for something as simple as this, avoid using loops if possible.

An alternative approach to what Jubobs has suggested is that you can consider that polynomial equation to be in factored form where each factor denotes a root located at that particular location. You can use poly to convert these factors into a polynomial equation, then use polyval to evaluate the expression at the point you want. First, generate your roots by linspace where the points vary from -(N-1)/2 to (N-1)/2 and there are N of them, then plug this into poly. Finally, for any values of x, put this into polyval with the output of poly. The advantage of this approach is that you can evaluate multiple points of x in a single sweep.

Going with what you have, you would simply do this:

p = poly(linspace(-(N-1)/2, (N-1)/2, N));
out = polyval(p, x);

With your example, supposing that N = 6, this would be the output of the first line:

p =

1.0000         0   -8.7500         0   16.1875         0   -3.5156

As such, this is saying that when we expand out (x-5/2)(x-3/2)(x-1/2)(x+1/2)(x+3/2)(x+5/2), we get:

x^6 - 8.75x^4 + 16.1875x^2 - 3.5156

If we take a look at the roots of this equation, this is what we get:

r = roots(p)

r =

   -2.5000
    2.5000
   -1.5000
    1.5000
   -0.5000
    0.5000

As you can see, each term corresponds to one factor in your polynomial equation, so we do have the right mindset here. Now, all you have to do is use p with your values of x into polyval to obtain your results. For example, if I wanted to evaluate that polynomial from -2 <= x <= 2 where x is an integer, this is the result I get:

polyval(p, -2:2)

ans =

  -14.7656    4.9219   -3.5156    4.9219  -14.7656

Therefore, when x = -2, the result is -14.7656 and so on.

7
  • Is there a benefit in favouring your alternative?
    – jub0bs
    Dec 8, 2014 at 21:16
  • @Jubobs - The only two things that someone may favour in this method over yours is that it's more readable. Also, this can take in multiple values of x and generate those corresponding output values. With yours, you can only take in one value of x at any given time. You could, however, place multiple values of x if you make prod work over a given dimension, but I'm just nitpicking.
    – rayryeng
    Dec 8, 2014 at 21:22
  • I guess I could make mine more readable/user-friendly by defining an anonymous function for it.
    – jub0bs
    Dec 8, 2014 at 21:23
  • 1
    @Jubobs - That's true!... either way, I did give you a +1 as it's more in tune with what the OP wants.
    – rayryeng
    Dec 8, 2014 at 21:24
  • thank you both for your answers and for your time. rayryeng answer was a starting point for discovering new methods of working in matlab. If I could I would vote both up (no rep yet)
    – Sobeit
    Dec 8, 2014 at 21:26
1

Though I would recommend the solution by @Jubobs, it is also good to check what the issue is with your loop.

The first indication that something is wrong, is that you have a nested loop over 2 variables, and only index with one of them to store the result. Probably you just need a single loop.

Here is a loop that you may be interested in that should do roughly what you need:

N=6;
x=1;
k=(-N-1)/2:(N-1)/2
P = ones(size(k));

for n=1:numel(k)
            P(n)=(x-k(n));
end
y=prod(P);

I tried to keep the code close to the original, so hopefully it is easy to understand.

1
  • 1
    I tried to keep the code close to the original [...] including the awful code formatting :) +1 anyway.
    – jub0bs
    Dec 10, 2014 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.