Need introductory example

[jeanqasaur]

This example needs to show:

• Why scientific computing is a mess with static types.
• How things look a lot better in Julia.

I've created a file introductory_example.tex with a skeleton for this, as well as prompts for parts of the writing.

Please work with @jiahao, @JeffBezanson, @StefanKarpinski, and @andreasnoack!

[jiahao]

Here is @andreasnoack's example which motivated him to work on Julia in the first place.

Fitting models to observed data is central to the practice of statistics. The fitting procedure estimates parameters to produce an optimal fit. Resampling methods (e.g. the bootstrap) are a class of computations in statistics use to compute the sensitivity of the fit to the values of the estimated parameters. At a high level, the pseudocode looks like this:

    for each resampling iteration do
        perturb the data slightly
        fit model to data
        save parameters
    end

    compute statistics of parameters

R has a very nice way to express the fit model to data step (in fact, it is one of the central innovations of S, the language that R was based on). In R, the fits are outsourced to packages like glmnet which are wrappers around C or Fortran code that do the heavy lifting.

The limitations of the two-language design are acute when trying to speed up this outer loop (which does the resampling) in R. There is no way to speed it up without reimplementing the entire kernel of the outer loop in C/Fortran. In the process one has to give up the convenient notation of R for the whole kernel.

From this example we can conclude that a slow high level language is not practical for scientific computing since one has to give up expressiveness for performance, and the high level language devolves into a glorified convenience shim layer around C/Fortran. Languages/libraries like NumPy and R that wrap C/Fortran under the hood produce a situation with limited composability and flexibility, since all the performance gains attained by wrapping C underneath evaporate once one wants to write a different computational kernel using that underlying C kernel as an inner loop, together with more code in the high level language.

[jiahao]

Upon reflection, this introductory example doesn't show off any multiple dispatch semantics. The multiplication operator is a better introductory example. In addition to the usual scalar_scalar, matrix_matrix examples, we can also show some nontrivial products like

1. Rayleigh quotients (for normalized vectors v) v*A*v:

which could make use of a specialized product method *(::Vector, ::Matrix, ::Vector) that knows that the result is a scalar to avoid allocating any temporary vectors (which would have necessary if using a pairwise product evaluation like v⋅(A*v) or (v*A)⋅v.

1. Projections of a vector w*w*v:

which can use a specialized product method *(::Vector, ::Vector, ::Vector) to ensure that the product is always evaluated in the order w*(w⋅v) even though the usual mathematical interpretation is the action of the projection operator (w*w') on the vector v. It is never advantageous to reify the matrix with the outer product first, since that means allocating the memory for a full matrix and furthermore can be numerically more unstable for floating-point products (since floating-point multiplies and adds are not associative, the order of operations matters).

We can also show how to define a custom ProjectionMatrix type that wraps w and whose product with a Vector can be defined with:

type ProjectionMatrix{T}
   w :: Vector{T}
end

*{T<:Number}(P::ProjectionMatrix{T}, v::Vector{T}) = P.w  * P.w * v

function *{T<:Real}(w::Vector{T}, yt::Vector{T}, v::Vector{T})
    (yt ⋅ v) * w
end

The disadvantage of these examples is that they punt on not distinguishing v from v' in the product, which seems a little awkward but inevitable due to JuliaLang/julia#4774. These examples would be cleaner if v' produced a Transpose type, but it doesn't.

[JeffBezanson]

Here's the qrfact example I was given that shows computing a result type in a non-trivial way:

function qrfact{T}(A::StridedMatrix{T}; pivot=false)
    S = typeof(one(T)/norm(one(T)))
    if S != T
        qrfact!(convert(AbstractMatrix{S},A), pivot=pivot)
    else
        qrfact!(copy(A),pivot=pivot))
    end
end

[JeffBezanson]

Edited example of a numpy routine that uses the same pattern of computing what type to use:

def vander(x, N):
  x = asarray(x)
  if x.ndim != 1:
    raise ValueError("x must be a one-dimensional array or sequence.")

  v = empty((len(x), N), dtype=promote_types(x.dtype, int))

  if N > 0:
    v[:, 0] = 1
  if N > 1:
    v[:, 1:] = x[:, None]
    multiply.accumulate(v[:, 1:], out=v[:, 1:], axis=1)

  return v