# Numerical methods challenge: Day 7

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.

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$.

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)))
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]])


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
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]


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.