Performance improvement for multiple uses of polyval function

Question

I have a simple performance question on using polyval function with Matlab.

Currently, I have a vector of x that can be quite long (>1000 scalars). I want to apply a different polynomial form to each of the x.

The polynomial forms are stored in a 2d array and applied in a loop like the code below. The code is relatively fast as polyval is optimized but the loop can be lengthy and performance is paramount as it is an objective function that can be computed thousands of times in a process.

Any idea on how to improve the performance?

Thanks

% ---------- Objective Function ------------------
function [obj] = obj(x, poly_objective)
    polyvalue = zeros(length(x),1);
    for index = 1: length(x)
        polyvalue (index) = polyval(poly_objective(index,:), x(index));
    end
    obj= -sum(polyvalue );
end
% -------------------------------------------------

What is the size of `poly_objective`? Is it assured that `length(x)` equals `size(poly_objective,1)`? Also, do you need the results of the individual polynomial of only their sum as per your code? — Luis Mendo, Oct 23 '18 at 14:19
Thank you for your response below. You are correct, I only need the sum. Would this change the question (and the answer)? — blouvard, Oct 23 '18 at 15:10
No, I don't see how to speed up my answer to compute the sum. Just obtain `polyvalue` and then `sum(polyvalue)` — Luis Mendo, Oct 23 '18 at 16:12

obchardon · Accepted Answer · 2018-10-23T15:00:13.570

1

You can linearize your for loop manually, here is an example:

p = [3,2,1;  
     5,1,3];           %polynomial coeff
x = [5,6].';           %the query points
d = size(p,2)-1:-1:0;  %the power factors
res = sum(x.^d.*p,2);  %process all the polynome without for loop.

with

Also if you would evaluate each x value for each polynome you could use:

res = x.^d*p.'; %using only matrix multiplication

with

res =
       p1    p2
x1     86    133
x2     121   189

edited Oct 23 '18 at 15:00

answered Oct 23 '18 at 14:32

obchardon

10,614
1
17
33

Thanks! Extremely helpful. 3 really good answers with different perspectives but converging views! – blouvard Oct 23 '18 at 15:06

Luis Mendo · Answer 2 · 2018-10-23T16:11:45.473

I find your question a little confusing, but I think this does what you want:

polyvalue = sum(poly_objective .* x(:).^(numel(x)-1:-1:0), 2);

Note that the above uses implicit expansion. For Matlab vesions before R2016b, use bsxfun:

polyvalue = sum(poly_objective .* bsxfun(@power, x(:), (numel(x)-1:-1:0)), 2);

Example

Random data:

>> x = rand(1,4);
>> poly_objective = randi(9,4,4);

Your code:

>> polyvalue = zeros(length(x),1);
   for index = 1: length(x)
       polyvalue (index) = polyval(poly_objective(index,:), x(index));
   end
>> polyvalue
polyvalue =
  13.545710504297881
  16.286929525147158
  13.289183623920710
   5.777980886766799

My code:

>> polyvalue = sum(poly_objective .* x(:).^(numel(x)-1:-1:0), 2)
polyvalue =
  13.545710504297881
  16.286929525147158
  13.289183623920710
   5.777980886766799

score 1 · Answer 3 · answered Oct 23 '18 at 14:36

The quickest way is likely to evaluate the different polynomials directly, removing the loop (as shown by obchardon or Luis). However, here's a note on polyval performance...

If you type edit polyval in the command window, you can see the source for the polyval function. In particular there is the following conditional evaluation near the top:

nc = length(p);
if isscalar(x) && (nargin < 3) && nc>0 && isfinite(x) && all(isfinite(p(:)))
    % Make it scream for scalar x.  Polynomial evaluation can be
    % implemented as a recursive digital filter.
    y = filter(1,[1 -x],p);
    y = y(nc);
    return
end

I think the "Make it scream" comment is the developer(s) telling us this is a very quick route through the function! Aside; it's also the best comment I've found in a MATLAB built-in.

So let's try to satisfy the conditions for this if statement...

✓ isscalar(x)
✓ nargin < 3 
✓ length(p) > 0
✓ isfinite(x)
✓ all(isfinite(p(:)))

Brilliant, so this is always the evaluation you're using. You might find speed improvements in removing these 5 checks, and simply doing this instead of polyval. In terms of your variables, this looks like so:

y = filter(1,[1 -x(index)],poly_objective(index,:)); 
polyvalue (index) = y(size(poly_objective,2)); 
% Note you should get size(poly_objective,2) outside your loop

_It's also the best comment I've found in a MATLAB built-in_ I agree :-) — Luis Mendo, Oct 23 '18 at 16:08

Performance improvement for multiple uses of polyval function

3 Answers3

Example