Finding a local max/min of a function in Python

[Eclair_de_XII]

Okay, so my algorithm looks something like this:
====
1. Locate mid-point of the interval . This is our estimate of the local max.
2. Evaluate .
3. Divide the main interval into two subintervals: a left and right of equal length.
4. Check to see if there is a point in the interval such that . If there is, check, because the maximum is there.
5. Do the same for .
6. If there is not a point in either interval such that greater than the current iteration of , then set . This is the maximum.
====

from numpy import vectorize, linspace

def find_max(x0,xf,h,N=100):
    xm=(xf+x0)/2; f=vectorize(h)
    xleft=linspace(x0,xm,N); xright=linspace(xm,xf,N)
    if (f(xleft)-f(xm)>1E-16).any:
        return find_max(x0,xm,h,N)
    if (f(xright)-f(xm)>1E-16).any:
        return find_max(xm,xf,h,N)
    else:
        return xm

t=(-1,1,lambda x: x**2)
find_max(*t)

Edit: I saw what was wrong with my code. I forgot to put the () after the .any method.

[CODE lang="python" title="Final version."]from numpy import vectorize, linspace, pi,sin

def diff(f,x1,x2):
return f(x1)-f(x2)

def find_crt(x0,xf,h,N=100,local='max'):
local.lower()
assert local=='min' or local=='max'
xm=(xf+x0)/2; f=vectorize(h)
xleft=linspace(x0,xm,N); xright=linspace(xm,xf,N)
t1=[xleft,xm]; t2=[xright,xm]
if local=='min':
t1,t2=t1[::-1],t2[::-1]
if (diff(f,*t1)>=1E-16).any():
return find_crt(x0,xm,h,N)
if (diff(f,*t2)>=1E-16).any():
return find_crt(xm,xf,h,N)
else:
return xm[/CODE]

It's odd, though, that whenever I use this function for this somewhat messy, amateur module:

https://github.com/gchang12/scientific_programming/blob/master/inheritance

with the parameters (x0=-pi/2,xf=pi/2,h=cos,N=50), the module I made implementing it returns no minimum points, while the module here returns a point exponentially close to zero (which I should maybe work on fixing).

 

[Dr.AbeNikIanEdL]

Your post does not really contain a question (anymore?), so I am not sure if you are still looking for help with this. Anyway, what strikes me a bit odd about your algorithm is that you only divide the interval into two at each step, but to do that, you compute the function on all 2N points. After doing that, you could just look what point actually had the maximum value and then limit the search to the neighborhood of that point. Or, if you really want to split the interval into two step by step, you should be able to achieve this with just two evaluations of the function per step (at the neighboring points of you initial guess). Is there any particular reason for the algorithm you chose?

 

[Eclair_de_XII]

I am not sure if you are still looking for help with this.

At this point, I am just wondering why it does not work in the module that imported the function.

After doing that, you could just look what point actually had the maximum value and then limit the search to the neighborhood of that point.

The problem is that a function evaluated on an array generated with the numpy.linspace function might not necessarily contain the maximum value of the function. The array is not "dense" in the interval $[x0,xf]$, to use the term very, very lightly.

Or, if you really want to split the interval into two step by step, you should be able to achieve this with just two evaluations of the function per step (at the neighboring points of you initial guess).

That, I admit would be a much simpler solution than checking the arrays if there is a maximum value. Checking if the mid-point of the interval is greater than the mid-point of either the left or right sub-interval, then checking it, if is not, would be much more straight-forward.

 

[Eclair_de_XII]

Odd; I just reread what I wrote, and it sounds a lot like a proof I read for that Bolzano-Weierstrauss Theorem.
===
Anyway, I wrote some code again, without using the numpy module. Yeah, your way is better.

def find_crt2(x0,xf,f,local='max'):
    """This function finds the maximum of an interval [x0,xf] that contains exactly one minimum or maximum."""
    # Want to make sure end-user specifies valid response to type of critical point.
    # local.lower(); assert local=='min' or local=='max' # Nulled code because it does not distinguish between min/max.
    # Want to initialize mid-point of main interval, and mid-point of left and right sub-intervals.
    xm=(xf+x0)/2; xl=(x0+xm)/2; xr=(xm+xf)/2
    # Crucial to returning min or max, per request.
    arg1=[xm,xl]; arg2=[xm,xr]
    # If user seeks minimum, then order of terms in difference is reversed.
    #if local=='min':
    #    arg1=arg1[::-1]; arg2=arg2[::-1]
    # Check if mid-point of main interval is greater than mid-point of left sub-interval.
    if diff(f,*arg1)>1E-16:
        # If it is, then check right sub-interval [xm,xf] because the left sub-interval does not contain max.
        return find_crt2(xm,xf,f)
    elif diff(f,*arg2)>1E-16:
        return find_crt2(x0,xm,f)
    else:
        # Return mid-point of current interval. This ought to be the maximum of the function.
        return xm

 

[Dr.AbeNikIanEdL]

The problem is that a function evaluated on an array generated with the numpy.linspace function might not necessarily contain the maximum value of the function. The array is not "dense" in the interval [x0,xf], to use the term very, very lightly.

Thats not what I said, I said, with new emphasis

you could just look what point actually had the maximum value and then limit the search to the neighborhood of that point.

i.e. if the maximum of the discrete version is at $x_i$ (at least if that is not the start or end point), then $\max{f}\in (x_{i-1}, x_{i+1})$ (otherwise, there would be more then one local max in that interval).

Of course, you could always have a constant function plus a very thin peak that would probably be hard to resolve by either method...

 

[Eclair_de_XII]

i.e. if the maximum of the discrete version is at (at least if that is not the start or end point), then $\max f\in (x_{i-1},x_{i+1})$ (otherwise, there would be more then one local max in that interval).

So something like this, I think?

def find_crt3(x0,xf,f,N=51):
    """This function searches an interval [x0,xf] containing exactly one minimum or maximum of a function f,
    and returns the critical point by recursively checking neighborhoods of estimated maximum."""
    f=vectorize(f); xinc=(xf-x0)/N # This will be the radius of the interval if recursion is necessary.
    xarr=linspace(x0,xf,N); yarr=f(xarr)
    ymax=max(yarr) # Setting the max y-value
    for x in xarr:
        if f(x)==ymax:
            xmax=x # Setting the corresponding x-value for max y-value
            break
    if abs(xmax-(xf+x0)/2)<=1E-16: # Assumed that local max occurs at mid-point of interval at final recursion.
        return xmax
    else: # Search in a neighborhood about the current point at which f is max.
        return find_crt3(xmax-xinc,xmax+xinc,f)

constant function plus a very thin peak

So like a function that is doing a terrible impression of the dirac-delta relation.

 

[Dr.AbeNikIanEdL]

So something like this, I think?

Kind of yeah. I think you can make better use of numpy functions to find the index there, but that's what I had in mind.With the "constant plus very thin peak" I was, in an extreme case, thinking of something like

def f(x):
    if x == 0.468: return 2.
    else: return 1.

 

[pbuk]

Hopefully you have now had more lessons where you have learned about ternary (three-point) and Fibonnaci (golden section) search.