maximum of a polynomial

Question

I have a polynomial of order N (where N is even). This polynomial is equal to minus infinity for x minus/plus infinity (thus it has a maximum). What I am doing right now is taking the derivative of the polynomial by using polyder then finding the roots of the N-1 th order polynomial by using the roots function in Matlab which returns N-1 solutions. Then I am picking the real root that really maximizes the polynomial. The problem is that I am updating my polynomial a lot and at each time step I am using the above procedure to find the maximizer. Therefore, the roots function takes too much of a computation time making my application slow. Is there a way either in Matlab or a proposed algorithm that does this maximization in a computationally efficient fashion( i.e. just finding one solution instead of N-1 solutions)? Thanks.

Edit: I would also like to know whether there is a routine in Matlab that only returns the real roots instead of roots which returns all real/complex ones.

Rapture

You mean to say that your polynomial coefficients change EACH time step ? – Cabrales

yes, the coefficients change at each time step, also the coefficients are real. – Rapture

I doubt you can do it more efficiently than that. How do your polynomial coefficients change? What is the order of magnitude of N? Can you post some code and the reason why you need to optimize it? (i.e. why you think it's worth to optimize this specific part of the code). – Subspecies

@pedrosorio, N is a design parameter. It is the order of the Taylor approximation in my application. According to the profiler roots takes more than 80% of the computation time. I am using this in the context of a statistical inference application. – Rapture

Unless there is some relation between the polynomials, being able to do this faster would mean that knowing the maximum of any Nth degree polynomial lets you find the maximum of any other (unrelated) Nth degree polynomial faster. – Causal

What about finding a local maxima instead of an global one? Is there a specifically fast way to do that? – Rapture

Answer 1

There is a file exchange submission by Steve Morris which finds all real roots of functions on a given interval. It does so by interpolating the polynomial by a Chebychev polynomial, and finding its roots.

You can modify the eig evaluation of the companion matrix in there, to eigs. This allows you to find only one (or a few) roots and save time (there's a fair chance it's also possible to compute the roots or extrema of a Chebychev analytically, although I could not find a good reference for that (or even a bad one for that matter...)).

Another attempt that you can make in speeding things up, is to note that polyder does nothing more than

Pprime = (numel(P)-1:-1:1) .* P(1:end-1);

for your polynomial P. Also, roots does nothing more than find the eigenvalues of the companion matrix, so you could find these eigenvalues yourself, which prevents a call to roots. This could both be beneficial, because calls to non-builtin functions inside a loop prevent Matlab's JIT compiler from translating the loop to machine language. This could otherwise give you a large speed gain (factors of 100 or more are not uncommon).

Branks

Answer 2

I think that you are probably out of luck. If the coefficients of the polynomial change at every time step in an arbitrary fashion, then ultimately you are faced with a distinct and unrelated optimisation problem at every stage. There is insufficient information available to consider calculating just a subset of roots of the derivative polynomial - how could you know which derivative root provides the maximum stationary point of the polynomial without comparing the function value at ALL of the derivative roots?? If your polynomial coefficients were being perturbed at each step by only a (bounded) small amount or in a predictable manner, then it is conceivable that you would be able to try something iterative to refine the solution at each step (for example something crude such as using your previous roots as starting point of a new set of newton iterations to identify the updated derivative roots), but the question does not suggest that this is in fact the case so I am just guessing. I could be completely wrong here but you might just be out of luck in getting something faster unless you can provide more information of have some kind of relationship between the polynomials generated at each step.

Cabrales

the maximizer probably does not move so much at each time step especially as algorithm continues running (the distribution that I am trying to maximize reaches an equilibrium), so I can use the older root as initial guess as you have stated. But at the beginning the polynomials are arbitrarily changing, so your method makes sense after a point. However, there is still the issue of N solutions instead of 1. The roots return even imaginary roots whereas I am searching for the maximizer on the real line. Thanks for the answer btw. – Rapture

is there a faster routine than the roots in matlab that only returns the real roots instead of all real/complex roots. Thanks. – Rapture

@YBE I am afraid I do not know – Cabrales

@YBE - No. There is no faster routine that returns ONLY real roots. Wanting magic to work is not sufficient for it to happen, even if you want it very much. – Marston