In [1]:

    

using Optim


    

In [3]:

    

using Plots
plotlyjs()
ezcontour(x, y, f) = begin
    X = repmat(x', length(y), 1)
    Y = repmat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> log(f([x,y])), X, Y)
    plot(x, y, Z, st=:contour)
end


    

    




WARNING: Method definition ezcontour(Any, Any, Any) in module Main at In[2]:3 overwritten at In[3]:3.






    
Out[3]:





ezcontour (generic function with 1 method)

In [4]:

    

f = x -> (x[1].^2 + x[2] - 11).^2 + (x[1] + x[2].^2 - 7).^2
sol = optimize(f, [0.0;0.0], NelderMead())
xs = sol.minimizer


    

    
Out[4]:





2-element Array{Float64,1}:
 3.00001
 1.99999

In [5]:

    

ezcontour(-5:0.02:5, -5:0.02:5, z -> f(z)[1])
scatter!([xs[1]],[xs[2]],label=false)


    

    
Out[5]:

In [7]:

    

""" Take one step of the Nelder-Mead method to optimize a 
continuous function with no derivatives.

    `nelder_mead_step!(f,X) -> X,type`

This modifies X
"""
function nelder_mead_step!(f,X)
    # just do the simplest thing that could be dramatically improved
    npts = size(X,2) # number of points
    n = size(X,1)    # dim
    assert(npts == n+1) # need a simplex
    fs = zeros(npts)
    for i in 1:npts # compute all function vals
        fs[i] = f(vec(X[:,i]))
    end
    
    p = sortperm(fs) # sort in increasing order
    
    Xmid = @view X[:,p[1:end-1]]
    xworst = @view X[:,p[end]]
    centroid = vec(mean(Xmid,2))
    
    xbar = alpha -> centroid + alpha*(xworst - centroid)
    
    fm1 = f(xbar(-1)) # unit step away from centroid
    step = 0.0
    if fs[p[1]] <= fm1 && fm1 <= fs[p[end-1]]
        # accept
        xworst[:] = xbar(-1)
        step = -1.0
        return X,step
    elseif fm1 < fs[p[1]] # we found a much better point 
        # try and improve
        fm2 = f(xbar(-2))
        if fm2 < fm1
            # we made things better, accept!
            xworst[:] = xbar(-2)
            step = -2.0
            return X,step
        else
            xworst[:] = xbar(-1)
            step = -1.0
            return X,step
        end
    else
        # fm1 > second worst point
        if fs[p[end-1]] <= fm1 && fm1 <= fs[p[end]]
            fmhalf = f(xbar(-1/2))
            if fmhalf <= fm1
                xworst[:] = xbar(-1/2)
                step = -1/2
                return X,step
            end
        else
            fhalf = f(xbar(1/2))
            if fhalf < fs[p[end]] # we got slighlty better
                xworst[:] = xbar(1/2)
                step = 1/2
                return X,step
            end
        end
    end
       
    assert(step == 0.0)
    # shrink towards best point
    xbest = @view X[:,p[1]]
    for i=1:npts
        if i != p[1]
            X[:,i] = 1/2*(X[:,i] + xbest)
        end
    end
    return X,step
end


    

    
Out[7]:





nelder_mead_step!

In [8]:

    

X = 2.5*randn(2,3)
nelder_mead_step!(f,X)


    

    
Out[8]:





(
[-0.880937 -1.98496 -1.10526; 2.4914 0.454804 3.14517],

0.0)

In [10]:

    

using Interact


    

In [11]:

    

nstart = 10
X0s = collect(map( x -> 2.5*randn(2,3), 1:nstart )) # generate 10 sets of starts
@manipulate for start=1:nstart, nstep=0:20
        ezcontour(-5:0.02:5, -5:0.02:5, z -> f(z)[1])
    X = copy(X0s[start])
    for j=1:nstep
        nelder_mead_step!(f, X)
        plot!(Shape([ (X[1,i], X[2,i]) for i in 1:size(X,2) ]), 
            fillalpha=0.2, fillcolor="grey")
    end
    
    plot!(Shape([ (X[1,i], X[2,i]) for i in 1:size(X,2) ]))
    title!(@sprintf("start = %i, step = %i", start, nstep))
    plot!(legend=false)
end


    

    












    












    
Out[11]:

In [48]:

    

X0


    

    
Out[48]:





2×3 Array{Float64,2}:
 -3.60791  -5.06751  -5.82868
 -2.19398  -3.47215  -3.31813