In January of this year I first saw mention of the Julia language in the release notes for LLVM. I mentioned this to Dirk Eddelbuettel and later we got in contact with Viral Shah regarding a Debian package for Julia.

There are many aspects of *Julia* that are quite intriguing to an
R programmer. I am interested in
programming languages for "Computing with Data", in John Chambers'
term, or "Technical Computing", as the authors of *Julia* classify it.
I believe that learning a programming language is somewhat like
learning a natural language in that you need to live with it and use
it for a while before you feel comfortable with it and with the
culture surrounding it.

A common complaint for those learning *R* is finding the name of the
function to perform a particular task. In writing a bit of *Julia*
code for fitting generalized linear models, as described below, I
found myself in exactly the same position of having to search through
documentation to find how to do something that I felt should be
simple. The experience is frustrating but I don't know of a way of
avoiding it. One word of advice for *R* programmers looking at
*Julia*, the names of most functions correspond to the
*Matlab*/*octave* names, not the *R* names. One exception is the
d-p-q-r functions for distributions, as I described in an earlier
posting.

## Evaluating a new language

There is a temptation to approach a new programming language with the
idea that it should be exactly like the language I am currently using
but better, faster, easier, ... This is, of course, ridiculous. If
you want *R* syntax, semantics, function names, and packages then use
*R*. If you consider another language you should accept that many
tasks will be done differently. Some will be done better than in your
current language and some will not be done as well. And the majority
will just be done differently.

Two primary *Julia* developers, Jeff Bezanson and Stefan Karpinski,
recently gave a presentation at the lang.NEXT conference. The
slides and
video
can help you get a flavor of the language.

So, what does *Julia* provide that *R* doesn't?

A JIT. The first time that a function is called with a particular argument signature it is compiled to machine code using the LLVM compiler backend. This, naturally, provides a speed boost. It can also affect programming style. In

*R*one instinctively writes vectorized code to avoid a performance hit. In*Julia*non-vectorized code can have similar and, under some circumstances, better performance than vectorized code.A richer type system. Scalars, vectors/arrays, tuples, composite types and several others can be defined in

*Julia*. In*R*the atomic types are actually atomic vectors and the composite types are either lists (which is actually a type of vector in*R*, not a lisp-like list), vectors with attributes (called structures in*R*but not to be confused with*C*'s structs). The class systems, S3 and S4, are built on top of these structures.Currently it is difficult to work with 64-bit integers in

*R*but not in*Julia*Multiple dispatch. In

*Julia*all functions provide multiple dispatch to handle different signatures of arguments. In some sense, functions are just names used to index into method tables. In*R*terminology all functions are generic functions with S4-like dispatch, not S3-like which only dispatches on the first argument.Parallelism. A multiple process model and distributed arrays are a basic part of the

*Julia*language.A cleaner model for accessing libraries of compiled code. As Jeff and Stefan point out, languages like

*R*and*Matlab*/*octave*typically implement the high-level operations in "glue" code that, in turn, accesses library code. In*Julia*one can short-circuit this process and call the library code from the high-level language. See the earlier blog posting on accessing the Rmath library from Julia. This may not seem like a big deal unless you have written thousands of lines of glue code, as I have.

There are some parts of *R* not present in *Julia* that programmers
may miss.

Default arguments of functions and named actual arguments in function calls. The process of matching actual arguments to formal arguments in

*R*function calls is richer. All matching in*Julia*is by position within function signature. Function design with many arguments requires more thought and care in*Julia*than in*R*.The

*R*package system and CRAN. One could also throw in namespaces and function documentation standards as part of the package system. These are important parts of*R*that are not yet available. However, that does not preclude simular facilities being developed for*Julia*. The*R*package system did not spring fully grown from the brow of Zeus.Graphics. One of the big advantages of

*R*over similar languages is the sophisticated graphics available in ggplot2 or lattice. Work is underway on graphics models for*Julia*but, again, it is early days still.

## An example: Generalized Linear Models (GLMs)

Most *R* users are familiar with the basic model-fitting functions
like lm() and glm(). These functions call model.frame and
model.matrix to form the numerical representation of the model from
the formula and data specification, taking into account optional
arguments such as subset, na.action and contrasts. After that they
call the numerical workhorses, lm.fit and glm.fit. What I will
describe is a function like glm.fit (the name in *Julia* cannot
include a period and I call it glmFit) and structures like the glm
"family" in *R*.

A glm family is a collection of functions that describe the distribution of the response, given the mean, and the link between the linear predictor and the mean response. The link function takes the mean response, mu, to the "linear predictor", eta, of the form

```
X %*% beta # in R
```

or

```
X * beta # in Julia
```

where X is the model matrix and beta is the coefficient vector. In practice we more frequently evaluate the inverse link, from eta to mu, and its derivative.

We define a composite type

```
type Link
name::String # name of the link
linkFun::Function # link function mu -> eta
linkInv::Function # inverse link eta -> mu
muEta::Function # derivative eta -> d mu/d eta
end
```

with specific instances

```
logitLink =
Link("logit",
mu -> log(mu ./ (1. - mu)),
eta -> 1. ./ (1. + exp(-eta)),
eta -> (e = exp(-abs(eta)); f = 1. + e; e ./ (f .* f)))
logLink =
Link("log",
mu -> log(mu),
eta -> exp(eta),
eta -> exp(eta))
identityLink =
Link("identity",
mu -> mu,
eta -> eta,
eta -> ones(eltype(eta), size(eta)))
inverseLink =
Link("inverse",
mu -> 1. ./ mu,
eta -> 1. ./ eta,
eta -> -1. ./ (eta .* eta))
```

for the common links. Here I have used the short-cut syntax

```
mu -> log(mu ./ (1. - mu))
```

to create anonymous functions. The muEta function for the logitLink
shows the use of multiple statements within the body of the anonymous
function. Like *R*, the last expression evaluated in a function is
the value of a *Julia* function.

I also ensure that the functions in the Link object are safe for vector arguments by using the componentwise-forms of multiplicative operators ("./" and ".*"). These forms apply to scalars as well as vector. I use explicit floating-point constants, such as 1., to give the compiler a hint that the result should be a floating-point scalar or vector.

A distribution characterizes the variance as a function of the mean, and the "deviance residuals" used to evaluate the model-fitting criterion. For some distributions, including the Gaussian distribution the model-fitting criterion (residual sum of squares) is different from the deviance itself but minimizing this simpler criterion provides the maximum-likelihood (or minimum deviance) estimates of the coefficients.

When we eventually fit the model we want the value of the deviance itself, which may not be the same as the sum of squared deviance residuals. Other functions are used to check for valid mean vectors or linear predictors. Also each distribution in what is called the "exponential family", which includes many common distributions like Bernoulli, binomial, gamma, Gaussian and Poisson, has a canonical link function derived from the distribution itself.

The composite type representing the distribution is

```
type Dist
name::String # the name of the distribution
canonical::Link # the canonical link for the distribution
variance::Function # variance function mu -> var
devResid::Function # vector of squared deviance residuals
deviance::Function # the scalar deviance
mustart::Function # derive a starting estimate for mu
validmu::Function # check validity of the mu vector
valideta::Function # check validity of the eta vector
end
```

with specific distributions defined as

```
## utilities used in some distributions
logN0(x::Number) = x == 0 ? x : log(x)
logN0{T<:Number}(x::AbstractArray{T}) = reshape([ logN0(x[i]) | i=1:numel(x) ], size(x))
y_log_y(y, mu) = y .* logN0(y ./ mu) # provides correct limit at y == 0
BernoulliDist =
Dist("Bernoulli",
logitLink,
mu -> max(eps(Float64), mu .* (1. - mu)),
(y, mu, wt)-> 2 * wt .* (y_log_y(y, mu) + y_log_y(1. - y, 1. - mu)),
(y, mu, wt)-> -2. * sum(y .* log(mu) + (1. - y) .* log(1. - mu)),
(y, wt)-> (wt .* y + 0.5) ./ (wt + 1.),
mu -> all((0 < mu) & (mu < 1)),
eta -> true)
gammaDist =
Dist("gamma",
inverseLink,
mu -> mu .* mu,
(y, mu, wt)-> -2 * wt .* (logN0(y ./ mu) - (y - mu) ./ mu),
(y, mu, wt)-> (n=sum(wt); disp=sum(-2 * wt .* (logN0(y ./ mu) - (y - mu) ./ mu))/n; invdisp(1/disp); sum(wt .* dgamma(y, invdisp, mu * disp, true))),
(y, wt)-> all(y > 0) ? y : error("non-positive response values not allowed for gammaDist"),
mu -> all(mu > 0.),
eta -> all(eta > 0.))
GaussianDist =
Dist("Gaussian",
identityLink,
mu -> ones(typeof(mu), size(mu)),
(y, mu, wt)-> (r = y - mu; wt .* r .* r),
(y, mu, wt)-> (n = length(mu); r = y - mu; n * (log(2*pi*sum(wt .* r .* r)/n) + 1) + 2 - sum(log(wt))),
(y, wt)-> y,
mu -> true,
eta -> true)
PoissonDist =
Dist("Poisson",
logLink,
mu -> mu,
(y, mu, wt)-> 2 * wt .* (y .* logN0(y ./ mu) - (y - mu)),
(y, mu, wt)-> -2 * sum(dpois(y, mu, true) * wt),
(y, mu)-> y + 0.1,
mu -> all(mu > 0),
eta -> true)
```

Finally the GlmResp type consists of the distribution, the link, the response, the mean and linear predictor. Occasionally we want to add an offset to the linear predictor expression or apply prior weights to the responses and these are included.

```
type GlmResp # response in a glm model
dist::Dist
link::Link
eta::Vector{Float64} # linear predictor
mu::Vector{Float64} # mean response
offset::Vector{Float64} # offset added to linear predictor (usually 0)
wts::Vector{Float64} # prior weights
y::Vector{Float64} # response
end
```

With just this definition we would need to specify the offset, wts, etc. every time we construct such an object. By defining an outer constructor (meaning a constructor that is defined outside the type definition) we can provide defaults

```
## outer constructor - the most common way of creating the object
function GlmResp(dist::Dist, link::Link, y::Vector{Float64})
n = length(y)
wt = ones(Float64, (n,))
mu = dist.mustart(y, wt)
GlmResp(dist, link, link.linkFun(mu), mu, zeros(Float64, (n,)), wt, y)
end
## another outer constructor using the canonical link for the distribution
GlmResp(dist::Dist, y::Vector{Float64}) = GlmResp(dist, dist.canonical, y)
```

The second outer constructor allows us to use the canonical link without needing to specify it.

There are several functions that we wish to apply to the glmResp object.

```
deviance( r::GlmResp) = r.dist.deviance(r.y, r.mu, r.wts)
devResid( r::GlmResp) = r.dist.devResid(r.y, r.mu, r.wts)
drsum( r::GlmResp) = sum(devResid(r))
muEta( r::GlmResp) = r.link.muEta(r.eta)
sqrtWrkWt(r::GlmResp) = muEta(r) .* sqrt(r.wts ./ variance(r))
variance( r::GlmResp) = r.dist.variance(r.mu)
wrkResid( r::GlmResp) = (r.y - r.mu) ./ r.link.muEta(r.eta)
wrkResp( r::GlmResp) = (r.eta - r.offset) + wrkResid(r)
```

Our initial definitions of these functions are actually method
definitions because we declare that the argument r must be a
GlmResp. However, defining a method also defines the function if it
has not already been defined. In *R* we must separately specify S3
generics and methods. An S4 generic is automatically created from a
method specification but the interactions of S4 generics and
namespaces can become tricky.

Note that the period ('.') is used to access components or members of a composite type and cannot be used in a Julia identifier.

The function to update the linear predictor requires both a GlmResp and the linear predictor value to assign. Here we use an abstract type inheriting from Number to allow a more general specification

```
updateMu{T<:Number}(r::GlmResp, linPr::AbstractArray{T}) = (r.eta = linPr + r.offset; r.mu = r.link.linkInv(r.eta); drsum(r))
```

Updating the linear predictor also updates the mean response then
returns the sum of the square deviance residuals. Note that this is
not functional semantics in that the object r being passed as an
argument is updated in place. Arguments to *Julia* functions are
passed by reference. Those accustomed to the functional semantics of
*R* (meaning that you can't change the value of an argument by passing
it to a function) should beware.

Finally we create a composite type for a predictor

```
type predD # predictor with dense X
X::Matrix{Float64} # model matrix
beta0::Vector{Float64} # base coefficient vector
delb::Vector{Float64} # increment
end
```

and an outer constructor

```
## outer constructor
predD(X::Matrix{Float64}) = (zz = zeros((size(X, 2),)); predD(X, zz, zz))
```

Given the current state of the predictor we create the weighted X'X matrix and the product of the weighted model matrix and the weighted working residuals in the accumulate function

```
function accumulate(r::GlmResp, p::predD)
w = sqrtWrkWt(r)
wX = diagmm(w, p.X)
(wX' * (w .* wrkResp(r))), (wX' * wX)
end
```

This function doesn't need to be written separately but I hope that doing so will make it easier to apply these operations to distributed arrays, which, to me, seem like a natural way of applying parallel computing to many statistical computing tasks.

A function to calculate and apply the increment is

```
increment(r::GlmResp, p::predD) = ((wXz, wXtwX) = accumulate(r, p); bb = wXtwX \ wXz; p.delb = bb - p.beta0; updateMu(r, p.X * bb))
```

and finally we get to the glmFit function

```
function glmFit(p::predD, r::GlmResp, maxIter::Uint, minStepFac::Float64, convTol::Float64)
if (maxIter < 1) error("maxIter must be positive") end
if (minStepFac < 0 || 1 < minStepFac) error("minStepFac must be in (0, 1)") end
cvg = false
devold = typemax(Float64) # Float64 version of Inf
for i=1:maxIter
dev = increment(r, p)
if (dev < devold)
p.beta0 = p.beta0 + p.delb
else
error("code needed to handle the step-factor case")
end
if abs((devold - dev)/dev) < convTol
cvg = true
break
end
devold = dev
end
if !cvg
error("failure to converge in $maxIter iterations")
end
end
glmFit(p::predD, r::GlmResp) = glmFit(p, r, uint(30), 0.001, 1.e-6)
```

Defining two methods allows for default values for some of the arguments. Here we have a kind of all-or-none approach to defaults. It is possible to create other method signatures for defaults applied to only some of the arguments. In retrospect I think I was being too cute in declaring the maxIter argument as an unsigned int. It only makes sense for this to be positive but it might make life simpler to allow it to be an integer as positivity is checked anyway.

Note the use of string interpolation in the last error message. The numeric value of maxIter will be substituted in the error message.

## Checking that it works

The script

```
load("../extras/glm.jl")
## generate a Bernoulli response
srand(123454321)
n = 10000
X = hcat(fill(1, (n,)), 3. * randn((n, 2)))
beta = randn((3,))
eta = X * beta
mu = logitLink.linkInv(eta)
y = float64(rand(n) < mu)
println("Coefficient vector for simulation: $(beta)")
## establish the glmResp object
rr = GlmResp(BernoulliDist, y)
pp = predD(X)
glmFit(pp, rr)
println("Converged to coefficient vector: $(pp.beta0 + pp.delb)")
println("Deviance at estimates: $(deviance(rr))")
## reinitialize objects for timing
rr = GlmResp(BernoulliDist, y)
pp = predD(X)
println("Elapsed time for fitting binary response with X $(size(X,1)) by $(size(X,2)): $(@elapsed glmFit(pp, rr)) seconds")
```

produces

```
Coefficient vector for simulation: [-0.54105, 0.308874, -1.02009]
Converged to coefficient vector: [-0.551101, 0.301489, -1.02574]
Deviance at estimates: 6858.466827205941
Elapsed time for fitting binary response with X 10000 by 3: 0.09874677658081055 seconds
```

so in summary...? benchmarks? do you like it more?

ReplyDeleteProgramming is very interesting and creative thing if you do it with love. Your blog code helps a lot to beginners to learn programming from basic to advance level. I really love this blog because I learn a lot from here and this process is still continuing.

ReplyDeleteLove from Pro Programmer