# Numerical methods challenge: Day 12

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.

## Hermite interpolation: Inverting Vandermonde matrix

Today we are comparing Lagrange interpolation with Hermite interpolation. For this example we are using the confluent Vandermonde matrix $V$. The code is based on an old code of mine that is present in this repo.

As in the case of Lagrange interpolation we solve the system

\begin{equation*} V\mathbf{c} = I \end{equation*}

where $\mathbf{c}$ is the vector of coefficients and $I$ is the identity matrix. This method is not stable and should not be used for the computation of higher order interpolants, even for optimally chosed sampling. It will start failing around 20 points. A better approach is to use the barycentric form of the interpolation.

In the example below we use Chebyshev nodes. The nodes are given by

\begin{equation*} x_k = \cos\left(\frac{2k-1}{2n}\pi\right), \quad k = 1, \ldots, n \end{equation*}

where $n$ is the degree of the polynomial.

Following are the codes.

### Python

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt

def vander_mat(x):
n = len(x)
van = np.zeros((n, n))
power = np.array(range(n))
for row in range(n):
van[row, :] = x[row]**power
return van

def conf_vander_mat(x):
n = len(x)
conf_van = np.zeros((2*n, 2*n))
power = np.array(range(2*n))
for row in range(n):
conf_van[row, :] = x[row]**power
conf_van[row + n, :] = power*x[row]**(power - 1)
return conf_van

def inter_coef(x, inter_type="lagrange"):
if inter_type == "lagrange":
vand_mat = vander_mat(x)
elif inter_type == "hermite":
vand_mat = conf_vander_mat(x)
coef = np.linalg.solve(vand_mat, np.eye(vand_mat.shape))
return coef

def compute_interp(x, f, x_eval, df=None):
n = len(x)
if df is None:
coef = inter_coef(x, inter_type="lagrange")
else:
coef = inter_coef(x, inter_type="hermite")
f_eval = np.zeros_like(x_eval)
nmat = coef.shape
for row in range(nmat):
for col in range(nmat):
if col < n or nmat == n:
f_eval += coef[row, col]*x_eval**row*f[col]
else:
f_eval += coef[row, col]*x_eval**row*df[col - n]
return f_eval

n = 7
x = -np.cos(np.linspace(0, np.pi, n))
f = lambda x: 1/(1 + 25*x**2)
df = lambda x: -50*x/(1 + 25*x**2)**2
x_eval = np.linspace(-1, 1, 500)
interp_f = compute_interp(x, f(x), x_eval, df=df(x))
plt.plot(x_eval, f(x_eval))
plt.plot(x_eval, interp_f)
plt.plot(x, f(x), ".")
plt.ylim(0, 1.2)
plt.show()


### Julia

using PyPlot

function vander_mat(x)
n = length(x)
van = zeros(n, n)
power = 0:n-1
for row = 1:n
van[row, :] = x[row].^power
end
return van
end

function conf_vander_mat(x)
n = length(x)
conf_van = zeros(2*n, 2*n)
power = 0:2*n-1
for row = 1:n
conf_van[row, :] = x[row].^power
conf_van[row + n, :] = power.*x[row].^(power - 1)
end
return conf_van
end

function inter_coef(x; inter_type="lagrange")
if inter_type == "lagrange"
vand_mat = vander_mat(x)
elseif inter_type == "hermite"
vand_mat = conf_vander_mat(x)
end
coef = vand_mat \ eye(size(vand_mat))
return coef
end

function compute_interp(x, f, x_eval; df=nothing)
n = length(x)
if df == nothing
coef = inter_coef(x, inter_type="lagrange")
else
coef = inter_coef(x, inter_type="hermite")
end
f_eval = zeros(x_eval)
nmat = size(coef)
for row = 1:nmat
for col = 1:nmat
if col <= n || nmat == n
f_eval += coef[row, col]*x_eval.^(row - 1)*f[col]
else
f_eval += coef[row, col]*x_eval.^(row - 1)*df[col - n]
end
end
end
return f_eval
end

n = 7
x = -cos.(linspace(0, pi, n))
f = 1./(1 + 25*x.^2)
df = -50*x./(1 + 25*x.^2).^2
x_eval = linspace(-1, 1, 500)
interp_f = compute_interp(x, f, x_eval, df=df)
plot(x_eval, 1./(1 + 25*x_eval.^2))
plot(x_eval, interp_f)
plot(x, f, ".")
ylim(0, 1.2)
show()


In both cases the result is the plot below.

And, if we try with a high $n$, say $n=43$, we can see the problems.

### Comparison Python/Julia

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

For Python:

%%timeit -n 100
n = 7
x = -np.cos(np.linspace(0, np.pi, n))
f = lambda x: 1/(1 + 25*x**2)
df = lambda x: -50*x/(1 + 25*x**2)**2
x_eval = np.linspace(-1, 1, 500)
interp_f = compute_interp(x, f(x), x_eval, df=df(x))


with result

100 loops, best of 3: 18.1 ms per loop


For Julia:

function bench()
n = 7
x = -cos.(linspace(0, pi, n))
f(x) = 1./(1 + 25*x.^2)
df(x) = -50*x./(1 + 25*x.^2).^2
x_eval = linspace(-1, 1, 500)
interp_f = compute_interp(x, f(x), x_eval, df=df(x))
end
@benchmark bench()


with result

BenchmarkTools.Trial:
memory estimate:  3.13 MiB
allocs estimate:  836
--------------
minimum time:     10.318 ms (0.00% GC)
median time:      10.449 ms (0.00% GC)
mean time:        11.362 ms (1.74% GC)
maximum time:     26.646 ms (0.00% GC)
--------------
samples:          100
evals/sample:     1


In this case, we can say that the Python code is roughly as fast as the Julia one.

### Comparison Hermite/Lagrange interpolation

We want to compare Hermite interpolation with Lagrange interpolation for the same number of degrees of freedom. We use the same function for test purposes

\begin{equation*} f(x) = \frac{1}{1 + 25x^2} \end{equation*}

This is the Python code that does that part

n_dof = np.array(range(1, 20))
error_herm = np.zeros(19)
error_lag = np.zeros(19)
for cont, n in enumerate(n_dof):
f = lambda x: 1/(1 + 25*x**2)
df = lambda x: -50*x/(1 + 25*x**2)**2
x = -np.cos(np.linspace(0, np.pi, n))
x2 = -np.cos(np.linspace(0, np.pi, 2*n))
x_eval = np.linspace(-1, 1, 500)
herm = compute_interp(x, f(x), x_eval, df=df(x))
lag = compute_interp(x2, f(x2), x_eval)
fun = f(x_eval)
error_herm[cont] = np.linalg.norm(fun - herm)/np.linalg.norm(fun)
error_lag[cont] = np.linalg.norm(fun - lag)/np.linalg.norm(fun)

plt.plot(2*n_dof, error_lag)
plt.plot(2*n_dof, error_herm)
plt.xlabel("Number of degrees of freedom")
plt.ylabel("Relative error")
plt.legend(["Lagrange", "Hermite"])
plt.show()


And this is the comparison of the relative errors

In general, the Lagrange approximation is closer for this function.

# Numerical methods challenge: Day 11

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.

## Lagrange interpolation: Inverting Vandermonde matrix

Today we have Lagrange interpolation, yet one more time. I will have a different approach to compute the interpolation; I will form the Vandermonde matrix $V$ and solve the system

\begin{equation*} V\mathbf{c} = I \end{equation*}

where $\mathbf{c}$ is the vector of coefficients and $I$ is the identity matrix. This method, and the previous one are not stable and should not be used for the computation of higher order interpolants, even for optimally chosed sampling. It will start failing around 40 points. A better approach is to use the barycentric form of the interpolation.

In the example below we use Chebyshev nodes. The nodes are given by

\begin{equation*} x_k = \cos\left(\frac{2k-1}{2n}\pi\right), \quad k = 1, \ldots, n \end{equation*}

where $n$ is the degree of the polynomial.

Following are the codes.

### Python

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt

def vander_mat(x):
n = len(x)
van = np.zeros((n, n))
power = np.array(range(n))
for row in range(n):
van[row, :] = x[row]**power
return van

def inter_coef(x):
vand_mat = vander_mat(x)
coef = np.linalg.solve(vand_mat, np.eye(len(x)))
return coef

def compute_interp(x, f, x_eval):
n = len(x)
coef = inter_coef(x)
f_eval = np.zeros_like(x_eval)
for row in range(n):
for col in range(n):
f_eval += coef[row, col]*x_eval**row*f[col]
return f_eval

n = 11
x = -np.cos(np.linspace(0, np.pi, n))
f = 1/(1 + 25*x**2)
x_eval = np.linspace(-1, 1, 500)
interp_f = compute_interp(x, f, x_eval)
plt.figure()
plt.plot(x_eval, 1/(1 + 25*x_eval**2))
plt.plot(x_eval, interp_f)
plt.plot(x, f, ".")
plt.ylim(0, 1.2)
plt.show()


### Julia

using PyPlot

function vander_mat(x)
n = length(x)
van = zeros(n, n)
power = 0:n-1
for row = 1:n
van[row, :] = x[row].^power
end
return van
end

function inter_coef(x)
vand_mat = vander_mat(x)
coef = vand_mat \ eye(length(x))
return coef
end

function compute_interp(x, f, x_eval)
n = length(x)
coef = inter_coef(x)
f_eval = zeros(x_eval)
for row = 1:n
for col = 1:n
f_eval += coef[row, col]*x_eval.^(row - 1)*f[col]
end
end
return f_eval
end

n = 11
x = - cos.(linspace(0, pi, n))
f = 1./(1 + 25*x.^2)
x_eval = linspace(-1, 1, 500)
interp_f = compute_interp(x, f, x_eval)
plot(x_eval, 1./(1 + 25*x_eval.^2))
plot(x_eval, interp_f)
plot(x, f, ".")
ylim(0, 1.2)
show()


In both cases the result is the plot below.

And, if we try with a high $n$, say $n=45$, we can see the problems.

### Comparison Python/Julia

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

For Python:

%%timeit -n 100
n = 11
x = -np.cos(np.linspace(0, np.pi, n))
f = 1/(1 + 25*x**2)
x_eval = np.linspace(-1, 1, 500)
interp_f = compute_interp(x, f, x_eval)


with result

100 loops, best of 3: 7.86 ms per loop


For Julia:

function bench()
x = - cos.(linspace(0, pi, n))
f = 1./(1 + 25*x.^2)
x_eval = linspace(-1, 1, 500)
interp_f = compute_interp(x, f, x_eval)
return nothing
end
@benchmark bench()


with result

BenchmarkTools.Trial:
memory estimate:  32.23 MiB
allocs estimate:  8277
--------------
minimum time:     114.282 ms (1.50% GC)
median time:      122.061 ms (1.46% GC)
mean time:        129.733 ms (1.90% GC)
maximum time:     163.716 ms (1.98% GC)
--------------
samples:          39
evals/sample:     1


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

# Numerical methods challenge: Day 10

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.

## Lagrange interpolation: Runge phenomenon

Today we have Lagrange interpolation, again. Technically, I am not posting about a different method, but just using the same algorithm for interpolation. The difference is that I will change the sampling, that is, I will use non-uniform sampling.

The problem with uniform interpolation is known as Runge phenomenon and is illustrated in the following image.

One way to mitigate the problem is to use non-uniform sampling, such as Chebyshev nodes or Lobatto nodes. The former set minimizes the Runge phenomenon, while the latter maximizes the Vandermonde determinant.

In the example below we use Lobatto sampling. Lobatto nodes are the zeros of

\begin{equation*} (1 - x^2) P'_N(x) \end{equation*}

where $P_N$ refers to the Nth Legendre polynomial. The use of these nodes is useful in numerical integration and spectral methods. Finding the zeroes of these polynomials might be cumbersome in general. Nevertheless, we use an approach originally implemented in MATLAB by Greg von Winckel that use Chebyshev nodes as first guess and then update this guess using Newton-Raphson method.

Following are the codes.

### Python

from __future__ import division
from numpy import (zeros_like, pi, cos, zeros, amax, abs,
array, linspace, prod)
import matplotlib.pyplot as plt

def lagrange(x_int, y_int, x_new):
y_new = zeros_like(x_new)
for xi, yi in zip(x_int, y_int):
y_new += yi*prod([(x_new - xj)/(xi - xj)
for xj in x_int if xi!=xj], axis=0)
return y_new

def gauss_lobatto(N, tol=1e-15):
x = -cos(linspace(0, pi, N))
P = zeros((N, N))  # Vandermonde Matrix
x_old = 2
while amax(abs(x - x_old)) > tol:
x_old = x
P[:, 0] = 1
P[:, 1] = x
for k in range(2, N):
P[:, k] = ((2 * k - 1) * x * P[:, k - 1] -
(k - 1) * P[:, k - 2]) / k
x = x_old - (x * P[:, N - 1] - P[:, N - 2]) / (N * P[:, N - 1])
print(x)
return array(x)

runge = lambda x: 1/(1 + 25*x**2)
x = linspace(-1, 1, 100)
x_int = linspace(-1, 1, 11)
x_int2 = gauss_lobatto(11)
x_new = linspace(-1, 1, 1000)
y_new = lagrange(x_int, runge(x_int), x_new)
y_new2 = lagrange(x_int2, runge(x_int2), x_new)
plt.plot(x, runge(x), "k")
plt.plot(x_new, y_new)
plt.plot(x_new, y_new2)
plt.legend(["Runge function",
"Uniform interpolation",
"Lobatto-sampling interpolation"])
plt.xlabel("x")
plt.ylabel("y")
plt.show()


### Julia

using PyPlot

function lagrange(x_int, y_int, x_new)
y_new = zeros(x_new)
for (xi, yi) in zip(x_int, y_int)
prod = ones(x_new)
for xj in x_int
if xi != xj
prod = prod.* (x_new - xj)/(xi - xj)
end
end
y_new += yi*prod
end
return y_new
end

function gauss_lobatto(N; tol=1e-15)
x = -cos.(linspace(0, pi, N))
P = zeros(N, N)  # Vandermonde Matrix
x_old = 2
while maximum(abs.(x - x_old)) > tol
x_old = x
P[:, 1] = 1
P[:, 2] = x
for k = 3:N
P[:, k] = ((2 * k - 1) * x .* P[:, k - 1] -
(k - 1) * P[:, k - 2]) / k
end
x = x_old - (x .* P[:, N] - P[:, N - 1]) ./ (N* P[:, N])
end
return x
end

runge(x) =  1./(1 + 25*x.^2)
x = linspace(-1, 1, 100)
x_int = linspace(-1, 1, 11)
x_int2 = gauss_lobatto(11)
x_new = linspace(-1, 1, 1000)
y_new = lagrange(x_int, runge(x_int), x_new)
y_new2 = lagrange(x_int2, runge(x_int2), x_new)
plot(x, runge(x), "k")
plot(x_new, y_new)
plot(x_new, y_new2)
legend(["Runge function",
"Uniform interpolation",
"Lobatto-sampling interpolation"])
xlabel("x")
ylabel("y")


In both cases the result is plot shown above.

# Numerical methods challenge: Day 9

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.

## Lagrange interpolation

Today we have Lagrange interpolation. This is defined by

\begin{equation*} L(x) = \sum_{j=0}^{k} y_j \prod_{m\neq j}\frac{x - x_m}{x_j - x_m} \end{equation*}

with $x$ the point/points where we want to interpolate.

We will test the method with a sigmoid

\begin{equation*} f(x) = \frac{1}{1 + e^{-x}} \end{equation*}

Following are the codes.

### Python

from __future__ import division
from numpy import zeros_like, linspace, exp, prod
import matplotlib.pyplot as plt

def lagrange(x_int, y_int, x_new):
y_new = zeros_like(x_new)
for xi, yi in zip(x_int, y_int):
y_new += yi*prod([(x_new - xj)/(xi - xj)
for xj in x_int if xi!=xj], axis=0)
return y_new

x_int = linspace(-5, 5, 11)
y_int = 1/(1 + exp(-x_int))
x_new = linspace(-5, 5, 1000)
y_new = lagrange(x_int, y_int, x_new)
plt.plot(x_int, y_int, "ok")
plt.plot(x_new, y_new)
plt.show()


### Julia

using PyPlot

function lagrange(x_int, y_int, x_new)
y_new = zeros(x_new)
for (xi, yi) in zip(x_int, y_int)
prod = ones(x_new)
for xj in x_int
if xi != xj
prod = prod.* (x_new - xj)/(xi - xj)
end
end
y_new += yi*prod
end
return y_new
end

x_int = linspace(-5, 5, 11)
y_int = 1./(1 + exp.(-x_int))
x_new = linspace(-5, 5, 1000)
y_new = lagrange(x_int, y_int, x_new)
plot(x_int, y_int, "ok")
plot(x_new, y_new)


In both cases the result is the following plot.

### Comparison Python/Julia

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

For Python:

%timeit lagrange(x_int, y_int, x_new)


with result

1000 loops, best of 3: 1.55 ms per loop


For Julia:

@benchmark newton_opt(rosen, rosen_grad, rosen_hess, [2.0, 1.0])


with result

BenchmarkTools.Trial:
memory estimate:  1.97 MiB
allocs estimate:  254
--------------
minimum time:     737.665 μs (0.00% GC)
median time:      811.633 μs (0.00% GC)
mean time:        916.450 μs (10.77% GC)
maximum time:     3.119 ms (64.40% GC)
--------------
samples:          5433
evals/sample:     1


In this case, we can say that the Python code is roughly 2 times slower than the Julia one, where probably I am not using the best approach for Julia.

# Numerical methods challenge: Day 8

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.

## Newton's method for optimization

Today we have Newton's method for optimization. The main difference of this method with gradient descent is the use of higher derivatives, namely, the Hessian matrix of the objective function. The use of higher derivatives provides information of the curvature besides the slope information used in gradient descent. The following image compare the path for gradient descent (green) and Newton's method (red).

Mathematically, the update is written as

\begin{equation*} \mathbf{x}_k = \mathbf{x}_{k-1} - H(\mathbf{x}_k)^{-1} \nabla f(\mathbf{x_k}) \end{equation*}

where $H(\mathbf{x}_k)$ is the Hessian matrix at step k.

We will test the method with the Rosenbrock's function

\begin{equation*} f(x_1, x_2) = (1-x_1)^2 + 100(x_2-{x_1}^2)^2 \end{equation*}

Following are the codes.

### Python

from __future__ import division, print_function
from numpy import array
from numpy.linalg import norm, solve

def newton_opt(fun, grad, hess, x, niter=50, gtol=1e-8, verbose=False):
msg = "Maximum number of iterations reached."
for cont in range(niter):
if verbose:
print("n: {}, x: {}, g: {}".format(cont, x, g))
x = x - solve(hess(x), g)
if norm(g) < gtol:
msg = "Extremum found with desired accuracy."
break
return x, fun(x), msg

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

return array([
-2*(1 - x) - 400*x*(x - x**2),
200*(x - x**2)])

def rosen_hess(x):
return array([[-400*(x-x**2) + 800*x**2 + 2, -400*x],
[-400*x, 200]])

print(newton_opt(rosen, rosen_grad, rosen_hess, [2.0, 1.0], verbose=True))


with result

n: 0, x: [2.0, 1.0], g: [ 2402.  -600.]
n: 1, x: [ 1.99833611  3.99334443], g: [  1.99888520e+00  -5.53708323e-04]
n: 2, x: [ 1.00055248  0.0055331 ], g: [ 398.44998412 -199.11443262]
n: 3, x: [ 1.00054972  1.00109974], g: [  1.09944359e-03  -1.52451385e-09]
n: 4, x: [ 1.         0.9999997], g: [  1.20876952e-04  -6.04384753e-05]
(array([ 1.,  1.]), 0.0, 'Extremum found with desired accuracy.'


### Julia

function newton_opt(fun, grad, hess, x; niter=50, gtol=1e-8, verbose=false)
msg = "Maximum number of iterations reached."
for cont = 1:niter
if verbose
println("n: $(cont), x:$(x), g: $(g)") end x = x - hess(x)\g g = grad(x) if norm(g) < gtol msg = "Extremum found with desired accuracy." break end end return x, fun(x), msg end function rosen(x) return (1 - x)^2 + 100*(x - x^2)^2 end function rosen_grad(x) return [-2*(1 - x) - 400*x*(x - x^2); 200*(x - x^2)] end function rosen_hess(x) return [-400*(x - x^2) + 800*x^2 + 2 -400*x; -400*x 200] end println(newton_opt(rosen, rosen_grad, rosen_hess, [2.0, 1.0], verbose=true))  with result n: 1, x: [2.0, 1.0], g: [2402.0, -600.0] n: 2, x: [1.99834, 3.99334], g: [1.99889, -0.000553708] n: 3, x: [1.00055, 0.0055331], g: [398.45, -199.114] n: 4, x: [1.00055, 1.0011], g: [0.00109944, -1.52451e-9] n: 5, x: [1.0, 1.0], g: [0.000120877, -6.04385e-5] ([1.0, 1.0], 0.0, "Extremum found with desired accuracy.")  ### Comparison Python/Julia Regarding number of lines we have: 34 in Python and 37 in Julia. The comparison in execution time is done with %timeit magic command in IPython and @benchmark in Julia. For Python: %timeit newton_opt(rosen, rosen_grad, rosen_hess, [2.0, 1.0])  with result 1000 loops, best of 3: 247 µs per loop  For Julia: @benchmark newton_opt(rosen, rosen_grad, rosen_hess, [2.0, 1.0])  with result BenchmarkTools.Trial: memory estimate: 5.48 KiB allocs estimate: 120 -------------- minimum time: 5.784 μs (0.00% GC) median time: 6.030 μs (0.00% GC) mean time: 6.953 μs (10.00% GC) maximum time: 572.279 μs (95.96% GC) -------------- samples: 10000 evals/sample: 6  In this case, we can say that the Python code is roughly 40 times slower than the Julia one. ### Comparison with gradient descent We see an improvement in the number of iterations compared with gradient descent, that is, we moved from 40 iterations to 4 iterations, even if we demand the method to have higher accuracy, $10^{-12}$, for example. The appearance of this faster convergence does not come for free, of course. When using Newton's method we have two major drawbacks: • We need to compute the Hessian of the function, that might prove really difficult, even if we have the analytical expression for our function. • We need to solve a system of equation in each iteration. # 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. ## 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$. 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 _ References ---------- ..  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<=fr and fr<f[-2]: new_verts = np.vstack((verts[:-1, :], xr)) # 4. Expansion elif fr<f: 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)**2 + 100*(x - x**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<=fr && fr<f[end - 1] new_verts = [verts[1:end-1, :]; xr] # 4. Expansion elseif fr<f 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)^2 + 100*(x - x^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. # Numerical methods challenge: Day 6 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. ## Gradient descent Today we have gradient descent method. The main idea in this method is to use the steepest direction to move from a point to a new one in the neighbourhood of a differentiable function. This is illustrated in the following image, where we can see that the steps are given perpendicular to the isocontours of the function. If we think the function as the topography of certain region, it would probably be safe to think the isocontour of the functios as the path that cows would like to walk, while the gradient direction the path for a mountain goat. Mathematically, the update is written as \begin{equation*} \mathbf{x}_k = \mathbf{x}_{k-1} - \gamma_k \nabla f(\mathbf{x_k}) \end{equation*} where $\gamma_k$ is the size of the kth step. Intuitively, we can see that the size of the step should decrease while we approach the extremum of the function. This is normally achieving by a line search. For gradient descent method, we could use the Barzilai-Borwein method: \begin{equation*} \gamma_{n} = \frac{(\mathbf x_{n} - \mathbf x_{n-1})^T[\nabla f(\mathbf x_{n}) - \nabla f(\mathbf x_{n-1})]}{||\nabla f(\mathbf x_{n}) - \nabla f(\mathbf x_{n-1})||^2} \end{equation*} We will test the method with the Rosenbrock's function \begin{equation*} f(x_1, x_2) = (1-x_1)^2 + 100(x_2-{x_1}^2)^2 \end{equation*} Gradient descent present some problems with this function, that is not convex. Nevertheless, if we use an initial point that is close enough to the minimum $[1, 1]$, we can achieve convergence. Following are the codes. ### Python from __future__ import division, print_function from numpy import array from numpy.linalg import norm def grad_descent(fun, grad, x, niter=50, gtol=1e-8, verbose=False): msg = "Maximum number of iterations reached." g_old = grad(x) gamma = 0.1 for cont in range(niter): if verbose: print("n: {}, x: {}, g: {}".format(cont, x, g_old)) dx = -gamma*g_old x = x + dx g = grad(x) dg = g - g_old g_old = g gamma = dx.dot(dg)/dg.dot(dg) if norm(g) < gtol: msg = "Extremum found with desired accuracy." break return x, fun(x), msg def rosen(x): return (1 - x)**2 + 100*(x - x**2)**2 def rosen_grad(x): return array([ -2*(1 - x) - 400*x*(x - x**2), 200*(x - x**2)]) print(grad_descent(rosen, rosen_grad, [2.0, 1.0], verbose=True))  with result n: 0, x: [2.0, 1.0], g: [ 2402. -600.] n: 1, x: [-238.2 61. ], g: [ -5.40030319e+09 -1.13356480e+07] n: 2, x: [ 1.87289769 61.50393131], g: [-43446.62297136 11599.23711122] n: 3, x: [ 1.87482914 61.50341566], g: [-43485.61077531 11597.68626767] n: 4, x: [ -0.2532018 62.07096513], g: [ 6277.59251909 12401.37079452] n: 5, x: [ 1.40217916e-02 6.25988648e+01], g: [ -353.0701476 12519.73363327] n: 6, x: [ 2.98797952e-04 6.30854769e+01], g: [ -9.53932693e+00 1.26170954e+04] n: 7, x: [ 3.49096082e-03 5.88633946e+01], g: [ -84.18892291 11772.67648837] n: 8, x: [ 0.42114221 0.46054405], g: [-48.8618897 56.6366569] n: 9, x: [ 0.66471507 0.17821457], g: [ 69.42537577 -52.72631034] n: 10, x: [ 0.50504193 0.29948111], g: [-9.96223909 8.88275049] n: 11, x: [ 0.52491812 0.28175867], g: [-2.25608389 1.24392746] n: 12, x: [ 0.53044731 0.27871006], g: [-0.37379949 -0.53285773] n: 13, x: [ 0.53133016 0.27996858], g: [-0.43934324 -0.46863181] n: 14, x: [ 0.53252827 0.28124656], g: [-0.4365402 -0.46795943] n: 15, x: [ 0.75411231 0.51877873], g: [ 14.56231057 -9.98132891] n: 16, x: [ 0.7050077 0.55243611], g: [-16.21302574 11.08005 ] n: 17, x: [ 0.73088975 0.53474821], g: [-0.69854407 0.10967699] n: 18, x: [ 0.73204207 0.53456728], g: [-0.14989113 -0.26366293] n: 19, x: [ 0.73228024 0.53498623], g: [-0.16984866 -0.24962496] n: 20, x: [ 0.73260201 0.53545913], g: [-0.16949786 -0.24931553] n: 21, x: [ 0.93339279 0.83080364], g: [ 14.95730865 -8.08369381] n: 22, x: [ 0.89610321 0.85095683], g: [-17.39715957 9.59117536] n: 23, x: [ 0.91610476 0.83992984], g: [-0.41766894 0.13638094] n: 24, x: [ 0.91659561 0.83976956], g: [-0.02823623 -0.07559088] n: 25, x: [ 0.91662795 0.83985612], g: [-0.03817128 -0.0701336 ] n: 26, x: [ 0.91667285 0.83993863], g: [-0.03814167 -0.07009733] n: 27, x: [ 0.99186712 0.97813179], g: [ 2.23273615 -1.13372137] n: 28, x: [ 0.98342663 0.98241763], g: [-6.0476644 3.05793919] n: 29, x: [ 0.98959509 0.97929862], g: [ -2.08791162e-02 3.50119013e-05] n: 30, x: [ 0.98961644 0.97929858], g: [-0.00409146 -0.00842531] n: 31, x: [ 0.9896206 0.97930713], g: [-0.00421464 -0.00835884] n: 32, x: [ 0.98963284 0.97933141], g: [-0.0042095 -0.00834897] n: 33, x: [ 0.99991222 0.99971914], g: [ 0.04194186 -0.02106056] n: 34, x: [ 0.99597256 1.00169739], g: [-3.88678974 1.9472097 ] n: 35, x: [ 0.99987194 0.99974387], g: [ -2.42382231e-04 -6.86507784e-06] n: 36, x: [ 0.99987219 0.99974388], g: [ -5.01566886e-05 -1.02746923e-04] n: 37, x: [ 0.99987224 0.99974398], g: [ -5.10336616e-05 -1.02258276e-04] n: 38, x: [ 0.99987255 0.99974461], g: [ -5.09073070e-05 -1.02006794e-04] n: 39, x: [ 0.99999998 0.99999996], g: [ 5.05854337e-06 -2.54465444e-06] n: 40, x: [ 0.99998735 1.00000631], g: [-0.01266881 0.00632184] (array([ 1., 1.]), 7.2961338114681859e-21, 'Extremum found with desired accuracy.')  ### Julia function grad_descent(fun, grad, x; niter=50, gtol=1e-8, verbose=false) msg = "Maximum number of iterations reached." g_old = grad(x) gamma = 0.1 for cont = 1:niter if verbose println("n:$(cont), x: $(x), g:$(g_old)")
end
dx = - gamma*g_old
x = x + dx
dg = g - g_old
g_old = g
gamma = dx' * dg / (dg' * dg)
if norm(g) < gtol
msg = "Extremum found with desired accuracy."
break
end
end
return x, fun(x), msg
end

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

return [-2*(1 - x) - 400*x*(x - x^2);
200*(x - x^2)]
end



with result

n: 1, x: [2.0, 1.0], g: [2402.0, -600.0]
n: 2, x: [-238.2, 61.0], g: [-5.4003e9, -1.13356e7]
n: 3, x: [1.8729, 61.5039], g: [-43446.6, 11599.2]
n: 4, x: [1.87483, 61.5034], g: [-43485.6, 11597.7]
n: 5, x: [-0.253202, 62.071], g: [6277.59, 12401.4]
n: 6, x: [0.0140218, 62.5989], g: [-353.07, 12519.7]
n: 7, x: [0.000298798, 63.0855], g: [-9.53933, 12617.1]
n: 8, x: [0.00349096, 58.8634], g: [-84.1889, 11772.7]
n: 9, x: [0.421142, 0.460544], g: [-48.8619, 56.6367]
n: 10, x: [0.664715, 0.178215], g: [69.4254, -52.7263]
n: 11, x: [0.505042, 0.299481], g: [-9.96224, 8.88275]
n: 12, x: [0.524918, 0.281759], g: [-2.25608, 1.24393]
n: 13, x: [0.530447, 0.27871], g: [-0.373799, -0.532858]
n: 14, x: [0.53133, 0.279969], g: [-0.439343, -0.468632]
n: 15, x: [0.532528, 0.281247], g: [-0.43654, -0.467959]
n: 16, x: [0.754112, 0.518779], g: [14.5623, -9.98133]
n: 17, x: [0.705008, 0.552436], g: [-16.213, 11.08]
n: 18, x: [0.73089, 0.534748], g: [-0.698544, 0.109677]
n: 19, x: [0.732042, 0.534567], g: [-0.149891, -0.263663]
n: 20, x: [0.73228, 0.534986], g: [-0.169849, -0.249625]
n: 21, x: [0.732602, 0.535459], g: [-0.169498, -0.249316]
n: 22, x: [0.933393, 0.830804], g: [14.9573, -8.08369]
n: 23, x: [0.896103, 0.850957], g: [-17.3972, 9.59118]
n: 24, x: [0.916105, 0.83993], g: [-0.417669, 0.136381]
n: 25, x: [0.916596, 0.83977], g: [-0.0282362, -0.0755909]
n: 26, x: [0.916628, 0.839856], g: [-0.0381713, -0.0701336]
n: 27, x: [0.916673, 0.839939], g: [-0.0381417, -0.0700973]
n: 28, x: [0.991867, 0.978132], g: [2.23274, -1.13372]
n: 29, x: [0.983427, 0.982418], g: [-6.04766, 3.05794]
n: 30, x: [0.989595, 0.979299], g: [-0.0208791, 3.50119e-5]
n: 31, x: [0.989616, 0.979299], g: [-0.00409146, -0.00842531]
n: 32, x: [0.989621, 0.979307], g: [-0.00421464, -0.00835884]
n: 33, x: [0.989633, 0.979331], g: [-0.0042095, -0.00834897]
n: 34, x: [0.999912, 0.999719], g: [0.0419419, -0.0210606]
n: 35, x: [0.995973, 1.0017], g: [-3.88679, 1.94721]
n: 36, x: [0.999872, 0.999744], g: [-0.000242382, -6.86508e-6]
n: 37, x: [0.999872, 0.999744], g: [-5.01567e-5, -0.000102747]
n: 38, x: [0.999872, 0.999744], g: [-5.10337e-5, -0.000102258]
n: 39, x: [0.999873, 0.999745], g: [-5.09073e-5, -0.000102007]
n: 40, x: [1.0, 1.0], g: [5.05854e-6, -2.54465e-6]
n: 41, x: [0.999987, 1.00001], g: [-0.0126688, 0.00632184]
([1.0, 1.0], 7.296133811468186e-21, "Root 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 grad_descent(rosen, rosen_grad, [2.0, 1.0])


with result

1000 loops, best of 3: 386 µs per loop


For Julia:

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


with result

BenchmarkTools.Trial:
memory estimate:  16.91 KiB
allocs estimate:  251
--------------
minimum time:     6.479 μs (0.00% GC)
median time:      7.393 μs (0.00% GC)
mean time:        13.437 μs (18.45% GC)
maximum time:     2.029 ms (95.94% GC)
--------------
samples:          10000
evals/sample:     5


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

# Numerical methods challenge: Day 5

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.

## Broyden's method

Today we have Broyden's method. The main idea in this method is to compute the Jacobian matrix just at the first iteration and change it for each iteration doing rank-1 updates.

\begin{equation*} \mathbf{x}_k = \mathbf{x}_{k-1} - J_n \mathbf{F}(\mathbf{x_k}) \end{equation*}

where we need estimate the Jacobian matrix at step $n$ by

\begin{equation*} \mathbf J_n = \mathbf J_{n - 1} + \frac{\Delta \mathbf f_n - \mathbf J_{n - 1} \Delta \mathbf x_n}{\|\Delta \mathbf x_n\|^2} \Delta \mathbf x_n^{\mathrm T} \end{equation*}

that correspond to rank-1 updates of the Jacobian matrix.

We will test the method with the function $\mathbf{F}(x, y) = (x + 2y - 2, x^2 + 4x - 4)$ with solution $\mathbf{x} = (0, 1)$.

Following are the codes.

### Python

from __future__ import division, print_function
from numpy import array, outer, dot
from numpy.linalg import solve, norm, det

def broyden(fun, jaco, x, niter=50, ftol=1e-12, verbose=False):
msg = "Maximum number of iterations reached."
J = jaco(x)
for cont in range(niter):
if det(J) < ftol:
x = None
msg = "Derivative near to zero."
break
if verbose:
print("n: {}, x: {}".format(cont, x))
f_old = fun(x)
dx = -solve(J, f_old)
x = x + dx
f = fun(x)
df = f - f_old
J = J + outer(df - dot(J, dx), dx)/dot(dx, dx)
if norm(f) < ftol:
msg = "Root found with desired accuracy."
break
return x, msg

def fun(x):
return array([x + 2*x - 2, x**2 + 4*x**2 - 4])

def jaco(x):
return array([
[1, 2],
[2*x, 8*x]])

print(broyden(fun, jaco, [1.0, 2.0]))


### Julia

function broyden(fun, jaco, x, niter=50, ftol=1e-12, verbose=false)
msg = "Maximum number of iterations reached."
J = jaco(x)
for cont = 1:niter
if det(J) < ftol
x = nothing
msg = "Derivative near to zero."
break
end
if verbose
println("n: $(cont), x:$(x)")
end
f_old = fun(x)
dx = -J\f_old
x = x + dx
f = fun(x)
df = f - f_old
J = J + (df - J*dx) * dx'/ (dx' * dx)
if norm(f) < ftol
msg = "Root found with desired accuracy."
break
end
end
return x, msg
end

function fun(x)
return [x + 2*x - 2, x^2 + 4*x^2 - 4]
end

function jaco(x)
return [1 2;
2*x 8*x]
end

println(broyden(fun, jaco, [1.0, 2.0]))


### 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 broyden(fun, jaco, [1.0, 2.0])


with result

1000 loops, best of 3: 703 µs per loop


For Julia:

@benchmark broyden(fun, jaco, [1.0, 2.0])


with result

BenchmarkTools.Trial:
memory estimate:  14.41 KiB
allocs estimate:  220
--------------
minimum time:     12.099 μs (0.00% GC)
median time:      12.867 μs (0.00% GC)
mean time:        15.378 μs (10.78% GC)
maximum time:     3.511 ms (97.53% GC)
--------------
samples:          10000
evals/sample:     1


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

### Comparison Newton/Broyden

Following, we are comparing Newton's and Broyden's method. We are using the function $\mathbf{x}^T \mathbf{x} + \mathbf{x}$ for this test.

#### Python

The code for the function and Jacobian is

from numpy import diag
fun = lambda x: x**2 + x
jaco = lambda x: diag(2*x + 1)


and the results are:

 n Newton (μs) Broyden (μs) 2 500 664 10 541 717 100 3450 4800

#### Julia

The code for the function and Jacobian is

fun(x) = x' * x + x
jaco(x) = diagm(2*x + 1)


and the results are:

 n Newton (μs) Broyden (μs) 2 1.76 1.65 10 56.42 5.12 100 1782 367

In this case, we are comparing the mean values of the results from @benchmark.

# Numerical methods challenge: Day 4

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.

## Newton's method: vector case

Today we have Newton's method one more time. In this case, the function is vector-valued. This implies a slight modification over the original post. The new approximation is computed from the old one using

\begin{equation*} \mathbf{x}_k = \mathbf{x}_{k-1} - J(\mathbf{x}_k)^{-1} \mathbf{F}(\mathbf{x_k}) \end{equation*}

where we need to use the Jacobian matrix $J$.

We will test the method with the function $\mathbf{F}(x, y) = (x + 2y - 2, x^2 + 4x - 4)$ with solution $\mathbf{x} = (0, 1)$.

Following are the codes.

### Python

from __future__ import division, print_function
from numpy import array
from numpy.linalg import solve, norm, det

def newton(fun, jaco, x, niter=50, ftol=1e-12, verbose=False):
msg = "Maximum number of iterations reached."
for cont in range(niter):
J = jaco(x)
f = fun(x)
if det(J) < ftol:
x = None
msg = "Derivative near to zero."
break
if verbose:
print("n: {}, x: {}".format(cont, x))
x = x - solve(J, f)
if norm(f) < ftol:
msg = "Root found with desired accuracy."
break
return x, msg

def fun(x):
return array([x + 2*x - 2, x**2 + 4*x**2 - 4])

def jaco(x):
return array([
[1, 2],
[2*x, 8*x]])

print(newton(fun, jaco, [1.0, 10.0]))


### Julia

function newton(fun, jaco, x, niter=50, ftol=1e-12, verbose=false)
msg = "Maximum number of iterations reached."
for cont = 1:niter
J = jaco(x)
f = fun(x)
if det(J) < ftol
x = nothing
msg = "Derivative near to zero."
break
end
if verbose
println("n: $(cont), x:$(x)")
end
x = x - J\f
if norm(f) < ftol
msg = "Root found with desired accuracy."
break
end
end
return x, msg
end

function fun(x)
return [x + 2*x - 2, x^2 + 4*x^2 - 4]
end

function jaco(x)
return [1.0 2.0;
2*x 8*x]
end

println(newton(fun, jaco, [1.0, 10.0]))


### Comparison

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

For Python:

%timeit newton(fun, jaco, [1.0, 10.0])


with result

1000 loops, best of 3: 284 µs per loop


For Julia:

@benchmark newton(fun, jaco, [1.0, 10.0])


with result

BenchmarkTools.Trial:
memory estimate:  10.44 KiB
allocs estimate:  192
--------------
minimum time:     6.818 μs (0.00% GC)
median time:      7.167 μs (0.00% GC)
mean time:        9.607 μs (16.53% GC)
maximum time:     2.953 ms (97.40% GC)
--------------
samples:          10000
evals/sample:     4


In this case, we can say that the Python code is roughly 40 times slower than the Julia one. This is an improvement compared to the previous examples, where the ratio was around 100. The reason for this "improvement" might be in the inversion of the Jacobian, that calls a numpy routine, doing the weight-lifting for us.

# Numerical methods challenge: Day 3

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.

## Newton's method

Today's method is Newton's method. This method is used to solve the equation $f(x) = 0$ for $x$ real, and $f$ is a differentiable funciton. It starts with an initial guess $x_0$ and it succesively refine it by finding the intercept of the tangent line to the function with zero. The new approximation is computed from the old one using

\begin{equation*} x_k = x_{k-1} - \frac{f(x)}{f'(x)} \end{equation*}

Convergence of this method is generally faster than bisection method. Nevertheless, the convergence is not guaranteed. Another drawback is the need of the derivative of the function.

We will use the function $f(x) = \cos(x) - x$ to test the codes, and the initial guess is 1.0.

Following are the codes.

### Python

from __future__ import division, print_function
from numpy import abs, cos, sin

def newton(fun, grad, x, niter=50, ftol=1e-12, verbose=False):
msg = "Maximum number of iterations reached."
for cont in range(niter):
x = None
msg = "Derivative near to zero."
break
if verbose:
print("n: {}, x: {}".format(cont, x))
if abs(fun(x)) < ftol:
msg = "Root found with desired accuracy."
break
return x, msg

def fun(x):
return cos(x) - x

return -sin(x) - 1.0



### Julia

function newton(fun, grad, x, niter=50, ftol=1e-12, verbose=false)
msg = "Maximum number of iterations reached."
for cont = 1:niter
x = nothing
msg = "Derivative near to zero."
break
end
if verbose
println("n: $(cont), x:$(x)")
end
if abs(fun(x)) < ftol
msg = "Root found with desired accuracy."
break
end
end
return x, msg
end

function fun(x)
return cos(x) - x
end

return -sin(x) - 1.0
end



### Comparison

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

For Python:

%timeit newton(fun, grad, 1.0)


with result

10000 loops, best of 3: 27.3 µs per loop


For Julia:

@benchmark newton(fun, grad, 1.0)


with result

BenchmarkTools.Trial:
memory estimate:  48 bytes
allocs estimate:  2
--------------
minimum time:     327.925 ns (0.00% GC)
median time:      337.956 ns (0.00% GC)
mean time:        351.064 ns (0.80% GC)
maximum time:     8.118 μs (92.60% GC)
--------------
samples:          10000
evals/sample:     226