Skip to main content

Nicolás' blog

Numerical methods challenge: Day 7

Nicolás Guarín-Zapata

Comments

Also available in:

Español

During October (2017) I will write a program per day for some well-known numerical methods in both Python and Julia. It is intended to be an exercise then don't expect the code to be good enough for real use. Also, I should mention that I have almost no experience with Julia, so it probably won't be idiomatic Julia but more Python-like Julia.

Nelder-Mead method

Today we have Nelder-Mead method. This method uses a simplex over \(\mathbb{R}^n\), e.g., a triangle in \(\mathbb{R}^2\) or a tetrahedron over \(\mathbb{R}^3\).

Animation showing Nelder-Mead method.

In a single step of the method, we remove the vertex with the worst function value and replace it with another one with a better value. The new point is obtained by reflecting, expanding, or contracting the simplex along the line joining the worst vertex with the centroid of the remaining vertices. If we cannot find a better point in this manner, we retain only the vertex with the best function value, and we shrink the simplex by moving all other vertices towards this value. Nocedal and Wright (2006), Numerical Optimization.

Following are the codes.

Python

from __future__ import division, print_function
import numpy as np
from numpy import array
from numpy.linalg import norm


def nelder_mead_step(fun, verts, alpha=1, gamma=2, rho=0.5,
                     sigma=0.5):
    """Nelder-Mead iteration according to Wikipedia _[1]


    References
    ----------
     .. [1] Wikipedia contributors. "Nelder–Mead method." Wikipedia,
         The Free Encyclopedia. Wikipedia, The Free Encyclopedia,
         1 Sep. 2016. Web. 20 Sep. 2016.
    """
    nverts, _ = verts.shape
    f = np.apply_along_axis(fun, 1, verts)
    # 1. Order
    order = np.argsort(f)
    verts = verts[order, :]
    f = f[order]
    # 2. Calculate xo, the centroid"
    xo = verts[:-1, :].mean(axis=0)
    # 3. Reflection
    xr = xo + alpha*(xo - verts[-1, :])
    fr = fun(xr)
    if f[0]<=fr and fr<f[-2]:
        new_verts = np.vstack((verts[:-1, :], xr))
    # 4. Expansion
    elif fr<f[0]:
        xe = xo + gamma*(xr - xo)
        fe = fun(xe)
        if fe < fr:
            new_verts = np.vstack((verts[:-1, :], xe))
        else:
            new_verts = np.vstack((verts[:-1, :], xr))
    # 5. Contraction
    else:
        xc = xo + rho*(verts[-1, :] - xo)
        fc = fun(xc)
        if fc < f[-1]:
            new_verts = np.vstack((verts[:-1, :], xc))
    # 6. Shrink
        else:
            new_verts = np.zeros_like(verts)
            new_verts[0, :] = verts[0, :]
            for k in range(1, nverts):
                new_verts[k, :] = sigma*(verts[k,:] - verts[0,:])

    return new_verts


def nelder_mead(fun, x, niter=200, atol=1e-8, verbose=False):
    msg = "Maximum number of iterations reached."
    f_old = fun(x.mean(0))
    for cont in range(niter):
        if verbose:
            print("n: {}, x: {}".format(cont, x.mean(0)))
        x = nelder_mead_step(fun, x)
        df = fun(x.mean(0)) - f_old
        f_old = fun(x.mean(0))
        if norm(df) < atol:
            msg = "Extremum found with desired accuracy."
            break
    return x.mean(0), f_old, msg


def rosen(x):
    return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2


x = array([[1, 0],
           [1, 1],
           [2, 0]])
print(nelder_mead(rosen, x))

with result

(array([0.99994674, 0.99987728]), 2.90769311467343e-08, 'Extremum found with desired accuracy.')

Julia

function nelder_mead_step(fun, verts; alpha=1, gamma=2, rho=0.5,
                     sigma=0.5)
    nverts, _ = size(verts)
    f = [fun(verts[row, :]) for row in 1:nverts]
    # 1. Order
    order = sortperm(f)
    verts = verts[order, :]
    f = f[order]
    # 2. Calculate xo, the centroid
    xo = mean(verts[1:end - 1, :], 1)
    # 3. Reflection
    xr = xo + alpha*(xo - verts[end, :]')
    fr = fun(xr)
    if f[1]<=fr && fr<f[end - 1]
        new_verts = [verts[1:end-1, :]; xr]
    # 4. Expansion
    elseif fr<f[1]
        xe = xo + gamma*(xr - xo)
        fe = fun(xe)
        if fe < fr
            new_verts = [verts[1:end-1, :]; xe]
        else
            new_verts = [verts[1:end-1, :]; xr]
        end
    # 5. Contraction
    else
        xc = xo + rho*(verts[end, :]' - xo)
        fc = fun(xc)
        if fc < f[end]
            new_verts = [verts[1:end-1, :]; xc]
    # 6. Shrink
        else
            new_verts = zeros(verts)
            new_verts[1, :] = verts[1, :]
            for k =  1:nverts
                new_verts[k, :] = sigma*(verts[k,:] - verts[1,:])
            end
        end
    end

    return new_verts
end


function nelder_mead(fun, x; niter=50, atol=1e-8, verbose=false)
    msg = "Maximum number of iterations reached."
    f_old = fun(mean(x, 1))
    for cont = 1:niter
        if verbose
            println("n: $(cont), x: $(mean(x, 1))")
        end
        x = nelder_mead_step(fun, x)
        df = fun(mean(x, 1)) - f_old
        f_old = fun(mean(x, 1))
        if norm(df) < atol
            msg = "Extremum found with desired accuracy."
            break
        end
    end
    return mean(x, 1), f_old, msg
end


function rosen(x)
    return (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2
end


x = [1 0;
    1 1;
    2 0]
println(nelder_mead(rosen, x, verbose=false))

with result

([0.999947 0.999877], 2.9076931147093985e-8, "Extremum found with desired accuracy.")

Comparison Python/Julia

Regarding number of lines we have: 38 in Python and 39 in Julia. The comparison in execution time is done with %timeit magic command in IPython and @benchmark in Julia.

For Python:

%timeit nelder_mead(rosen, x)

with result

100 loops, best of 3: 7.82 ms per loop

For Julia:

@benchmark grad_descent(rosen, rosen_grad, [2.0, 1.0])

with result

BenchmarkTools.Trial:
  memory estimate:  162.23 KiB
  allocs estimate:  4780
  --------------
  minimum time:     462.926 μs (0.00% GC)
  median time:      506.511 μs (0.00% GC)
  mean time:        552.411 μs (3.86% GC)
  maximum time:     5.179 ms (80.31% GC)
  --------------
  samples:          9008
  evals/sample:     1

In this case, we can say that the Python code is roughly 15 times slower than the Julia one.

Comments

Comments powered by Disqus