The following function(s) realize the trapezoidal rule for integrating a function along a discrete set of points representing the curve.
We implement an unvectorized and a vectorized form and attempt to speed it up.

In [1]:

# Discretized sine function
x = linspace(0, pi, 100);
y = sin(x);

Here is a vectorized form, i.e. how you would code it in MATLAB or R (or Python?).

In [2]:

function trapz1(x, y)
   if length(y) != length(x)
       error("Vectors must be of same length")
   end
   sum( (x[2:end] - x[1:end-1]).*(y[2:end] + y[1:end-1]) ) / 2
end

Out[2]:

trapz1 (generic function with 1 method)

In [3]:

println(trapz1(x, y)); gc()
@time [trapz1(x, y) for i in 1:1000];

1.9998321638939924
elapsed time: 0.020384185 seconds (6921872 bytes allocated)

And this is an unvectorized, explicit coding of trapezoidal integration.

In [4]:

function trapz2(x, y)
    local n = length(x)
    if length(y) != n
        error("Vectors 'x', 'y' must be of same length")
    end
    r = 0
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    r / 2
end

Out[4]:

trapz2 (generic function with 1 method)

In [5]:

println(trapz2(x, y)); gc()
@time [trapz2(x, y) for i in 1:1000];

1.9998321638939929
elapsed time: 0.009617445 seconds (3215904 bytes allocated)

Why is the following redefinition of trapz2 different (and faster)?

In [6]:

function trapz3(x, y)
    local n = length(x)
    if length(y) != n
        error("Vectors 'x', 'y' must be of same length")
    end
    r = 0.0
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    r / 2
end

Out[6]:

trapz3 (generic function with 1 method)

In [7]:

println(trapz3(x, y)); gc()
@time [trapz3(x, y) for i in 1:1000];

1.9998321638939929
elapsed time: 0.001451867 seconds (47904 bytes allocated)

Results and comparison with R and Python

          Timings      Result               µs/loop
trapz1  0.020384185  1.9998321638939924       20.4
trapz2  0.009617445  1.9998321638939929        9.6
trapz3  0.001451867  1.9998321638939929        1.5
trapz   0.000740450  1.9998321638939929        0.75

R:        unvect. 300 µs,  vect. 130 µs,  inline  --
Python:   unvect. 119 µs,  vect.  39 µs
+numba:           < 1 µs,         54 µs

A final version of trapz, with error handling, type stability, and omitting boundary checking, will be shown later.

Automated Differentiation (AD)

For explicitely given functions (as formulas), Julia provides analytical differentiation. Automatic differentiation (AD) is a technique to numerically evaluate the derivative of a function implemented as a program in a programming language, with the aim to compute derivatives to machine accuracy. This is important for solving (partial) diffeential equations, nonlinear optimization, sensitivity analysis, or design engineering.

The following example uses the lambertW function, the inverse function of x -> x * exp(x). This function is best be calculated through the Newton-Raphson-Haley method, an iteratve approach. We will apply automated differentiation -- defined in package DualNumbers -- to this iterative procedure to determine exactly the derivative of the Lambert_W function.

In [8]:

using DualNumbers

In [9]:

function lambertW(x::Number)
    if x <  -1.0/exp(1.0); return NaN, NaN;  end
    if x == -1.0/exp(1.0); return -Inf, Inf; end

    w0 = 1.0
    w1 = w0 - (w0 * exp(w0) - x)/((w0 + 1.0) * exp(w0) - 
        (w0 + 2.0) * (w0 * exp(w0) - x)/(2.0 * w0 + 2.0))

    while w0 != w1
        w0 = w1
        w1 = w0 - (w0 * exp(w0) - x)/((w0 + 1.0) * exp(w0) - 
            (w0 + 2.0) * (w0 * exp(w0) - x)/(2.0 * w0 + 2.0))
    end

    return w1
end

Out[9]:

lambertW (generic function with 1 method)

The derivative can be calculated directly because it is the inverse of $x \to x e^x$
(remember your analysis course).

In [10]:

function d_lambertW(x::Real)
    w1 = lambertW(x)
    return 1.0 / (1 + w1) / exp(w1)
end;

Now compute the exact derivative at 1.5 and compare with automated differentiation.

In [11]:

x0 = 1.5;
lambertW(1.5), d_lambertW(1.5)

Out[11]:

(0.7258613577662263,0.2803861212064392)

In [12]:

lambertW(Dual(x0, 1.0))  # epsilon(...)

Out[12]:

0.7258613577662263 + 0.2803861212064392du

Compare this with numerical derivatives calculates through the "central difference" formula:

In [13]:

# using Calculus
# Calculus.finite_difference(lambertW, 1.5)

Modelling optimization problems

The following optimization packages are available for Julia:

Optim w/ Nelder-Mead, (L-)BFGS, CG, etc.
NLopt w/, e.g., augmented Lagrangian
Convex programm.: ECOS, CVX
COIN-OR projects: Cbc, Clp, Ipopt
Global Optimiz.: BlackBoxOptim w/ differential evolution (DE), CMAES
commercial prog.: CPLEX, Gurobi, Mosek[, SCIP]
Modeling: JuMP, MathProgBase, CVX

JuMP is an algebraic modeling language for Mathematical Programming in Julia. It supports linear (LP) and quadratic (QP) optimization problems with linear constraints (but not quadratic constraints) and mixed-integer linear problems (MILP).
It interfaces to a number of open-source and commercial solvers (Clp, Cbc, GLPK, Gurobi, MOSEK, and CPLEX) via a generic solver-independent interface.
JuMP is extremely fast in generating a model, even compared to commercial solvers (like AMPL or GAMS).

As an example, we will maximize the following expression $5x_1 + 3.5x_2 + 4.5x_3$ with constraints \[3x_1 + 5x_2 + 4x_3 \le 540\] \[6x_1 + x_2 + 3x_3 \le 480\]

In [14]:

using JuMP, Cbc
m = Model( solver=CbcSolver() );

In [15]:

@defVar( m, x[1:3] >= 0)
@setObjective( m, Max, 5x[1] + 3.5x[2] + 4.5x[3] )

Out[15]:

5 x[1] + 3.5 x[2] + 4.5 x[3]

In [16]:

# @addConstraint( m, 3x[1] +  5x[2] +  4x[3] <= 540 )
# @addConstraint( m, 6x[1] +   x[2] +  3x[3] <= 480 )

A = [3. 5. 4.; 6. 1. 3.];
b = [540., 480.];

for i = 1:2
    @addConstraint( m, sum{A[i, j]*x[j], j=1:3} <= b[i] )
end

In [17]:

Out[17]:

$$ \begin{alignat*}{1} \max \quad &5 x[1] + 3.5 x[2] + 4.5 x[3]\\ \text{Subject to} \quad& 3 x[1] + 5 x[2] + 4 x[3] <= 540 \\ & 6 x[1] + x[2] + 3 x[3] <= 480 \\ & x_{i} \geq 0 \quad \forall i \in \{ 1..3 \} \\ \end{alignat*} $$

In [18]:

status = solve(m);
println("Objective value: ", getObjectiveValue(m))
println(getValue(x))

WARNING: Dual solutions not available

Objective value: 640.0
x
[1] = 20.0
[2] = 0.0
[3] = 120.0

Knapsack problem:
As a second problem, we solve a knapsack problem: Given the odd square numbers $1^2, ..., 43^2$, are there numbers that sum up to exactly 2014?

In [19]:

w = [1:2:43].^2;
n = length(w)     # n = 22

Out[19]:

In [20]:

m = Model( solver=CbcSolver() );
@defVar( m, x[1:n], Bin );

@setObjective( m, Max, sum{w[i]*x[i], i=1:n});

In [21]:

@addConstraint(m, sum{w[i]*x[i], i=1:n} <= 2014);

In [22]:

status = solve(m);
println("Objective value: ", getObjectiveValue(m))
println(getValue(x))

Objective value: 2014.0
x
[ 1] = 1.0
[ 2] = 0.0
[ 3] = 0.0
[ 4] = 0.9999999999999999
[ 5] = 0.9999999999999999
[ 6] = 0.0
[ 7] = 0.0
[ 8] = 0.0
[ 9] = 0.0
[10] = 0.0
[11] = 0.0
[12] = 1.0
[13] = 1.0
[14] = 0.9999999999999999
[15] = 0.0
[16] = 0.0
[17] = 0.0
[18] = 0.0
[19] = 0.0
[20] = 0.0
[21] = 0.0
[22] = 0.0

Get and check the result:

In [23]:

idx = [iround(getValue(x[i])) for i in 1:n]
println(idx)

{1,0,0,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0}

In [24]:

# sum(w[[1,4,5,12,13,14]]) == 2014  # = 1 + 7^2 + 9^2 + 23^2 + 25^2 + 27^2
s = w[find(idx)];
println(s)
println(sum(s))

[1,49,81,529,625,729]
2014

Applications

How to speed-up functions

Automated Differentiation (AD)

Modelling optimization problems