Python finding local maxima/minima for multiple polynomials efficiently

user_7942208 0
Viewed 114

I am looking for an efficient way to find local mins for multiple (>1 million) but independent 4th order polynomials in given/ specified ranges/ boundaries.

I have two requirements:

R1: efficient even for 1 million different polynomial equations

R2: the local min is accurate up to 0.01 (i.e. 2dp)

Here is some code I have created using scipy. It's okay but I am wondering if there's any other better packages in performing such a task before I go for parallel programming.

To illustrate my problem, let's start with one polynomial first:

Below I am trying to find the local min of 4x^4 + 6x^3 + 3x^2 + x + 5 within the range (-5, 5).

On my laptop, it takes about 2ms to find the local min (which is at ~ -0.72770502).

The time is alright for one polynomial but I would want something faster as I need to perform this operation over 1 million times regularly.

from scipy import optimize
import numpy as np

# Define a objective and gradient function for 4th order polynomial
# x is the value to be evaluated
# par is a numpy array of len 5 that specifies the polynomial coefficients.
def obj_grad_fun_custom(x,par):
    obj = (np.array([x**4,x**3,x**2,x**1,1]) * par).sum()
    grad = (np.array([4*x**3,3*x**2,2*x,1]) * par[:-1]).sum()
    return obj, grad

# Try minimise an example polynomial of 4x^4 + 6x^3 + 3x^2 + x + 5
# with contrainted bound
res = optimize.minimize(
    fun = obj_grad_fun_custom,
    x0 = 0,
    args=(np.array([4,6,3,1,5])), # polynomial coefficients
    jac=True ,
    bounds=[(-2, 10)],
    tol=1e-10)
print(res.x)

# Timing (this takes about 2 ms for me)
%timeit optimize.minimize(fun = obj_grad_fun_custom, x0 = 0, args=(np.array([4,6,3,1,5])), jac=True, bounds=[(-5, 5)], tol=1e-10)

Below is what I am planning to do regular with 1 million different order 4 polynomials I would want to minimise locally. Hopefully, someone could point me to a more suitable package other than scipy. Or any alternative methods? Thanks!

# Multiple polynomials
result = [] # saving the local minima
poly_sim_no = 1000000 #ideally 1 million or even more
np.random.seed(0)
par_set = np.random.choice(np.arange(10), size=(poly_sim_no, 5), replace=True) #generate some order 4 polynomial coefficients 

for a in par_set:
    res = optimize.minimize(obj_grad_fun_custom, 0,args=(a), jac=True ,bounds=[(-5, 5)], tol=1e-10)
    result.append(res.x)

print(result)
polynomials scipy optimization python
1 Answers

Since you're finding the minimum of a polynomial, you can take advantage of the fact that it's easy to take the derivative of a polynomial, and that there are many good algorithms for finding the roots of a polynomial.

Here's how it works:

  1. First, take the derivative. All of the points which are minimums will have a derivative of zero.
  2. Look for those zeros, aka find the roots of the derivative.
  3. Once we have the list of candidates, check that the solutions are real.
  4. Check that the solutions are within the bounds you set. (I don't know if you added bounds because you actually want the bounds, or to make it go faster. If it's the latter, feel free to remove this step.)
  5. Actually evaluate the candidates with the polynomial and find the smallest one.

Here's the code:

import numpy as np
from numpy.polynomial import Polynomial

def find_extrema(poly, bounds):
    deriv = poly.deriv()
    extrema = deriv.roots()
    # Filter out complex roots
    extrema = extrema[np.isreal(extrema)]
    # Get real part of root
    extrema = np.real(extrema)
    # Apply bounds check
    lb, ub = bounds
    extrema = extrema[(lb <= extrema) & (extrema <= ub)]
    return extrema

def find_minimum(poly, bounds):
    extrema = find_extrema(poly, bounds)
    # Note: initially I tried taking the 2nd derivative to filter out local maxima.
    # This ended up being slower than just evaluating the function.

    # Either bound could end up being the minimum. Check those too.
    extrema = np.concatenate((extrema, bounds))
    # Check every candidate by evaluating the polynomial at each possible minimum,
    # and picking the minimum.
    value_at_extrema = poly(extrema)
    minimum_index = np.argmin(value_at_extrema)
    return extrema[minimum_index]

# Warning: polynomial expects coeffients in the opposite order that you use.
poly = Polynomial([5,1,3,6,4]) 
print(find_minimum(poly, (-5, 5)))

This takes 162 microseconds on my computer, making it about 6x faster than the scipy.optimize solution. (The solution shown in the question takes 1.12 ms on my computer.)

Edit: A faster alternative

Here's a faster approach. However, it abandons bounds checking, uses a deprecated API, and is generally harder to read.

p = np.poly1d([4,6,3,1,5])  # Note: polynomials are opposite order of before

def find_minimum2(poly):
    roots = np.real(np.roots(poly.deriv()))
    return roots[np.argmin(poly(roots))]

print(find_minimum2(p))

This clocks in at 110 microseconds, making it roughly 10x faster than the original.

user_5301600
Related