An R programmer looks at Julia
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
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 Shahregarding 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 Juliacode 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 dpqr 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 videocan 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 nonvectorized 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 lisplike list), vectors with attributes (called structures in Rbut 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 64bit 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 S4like dispatch, not S3like 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/octavetypically implement the highlevel operations in “glue” code that, in turn, accesses library code. In Julia one can shortcircuit this process and call the library code from the highlevel 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 Juliais 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 modelfitting 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<br />
or
X * beta # in Julia<br />
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<br /> name::String # name of the link<br /> linkFun::Function # link function mu > eta<br /> linkInv::Function # inverse link eta > mu<br /> muEta::Function # derivative eta > d mu/d eta<br />end<br />
with specific instances
logitLink =<br /> Link("logit",<br /> mu > log(mu ./ (1.  mu)),<br /> eta > 1. ./ (1. + exp(eta)),<br /> eta > (e = exp(abs(eta)); f = 1. + e; e ./ (f .* f)))<br /><br />logLink =<br /> Link("log",<br /> mu > log(mu),<br /> eta > exp(eta),<br /> eta > exp(eta))<br /><br /><br />identityLink =<br /> Link("identity",<br /> mu > mu,<br /> eta > eta,<br /> eta > ones(eltype(eta), size(eta)))<br /><br />inverseLink =<br /> Link("inverse",<br /> mu > 1. ./ mu,<br /> eta > 1. ./ eta,<br /> eta > 1. ./ (eta .* eta))<br />
for the common links. Here I have used the shortcut syntax
mu > log(mu ./ (1.  mu))<br />
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 componentwiseforms of multiplicative operators (“./” and “.*”). These forms apply to scalars as well as vector. I use explicit floatingpoint constants, such as 1., to give the compiler a hint that the result should be a floatingpoint scalar or vector.
A distribution characterizes the variance as a function of the mean, and the “deviance residuals” used to evaluate the modelfitting criterion. For some distributions, including the Gaussian distribution the modelfitting criterion (residual sum of squares) is different from the deviance itself but minimizing this simpler criterion provides the maximumlikelihood (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<br /> name::String # the name of the distribution<br /> canonical::Link # the canonical link for the distribution<br /> variance::Function # variance function mu > var<br /> devResid::Function # vector of squared deviance residuals<br /> deviance::Function # the scalar deviance<br /> mustart::Function # derive a starting estimate for mu<br /> validmu::Function # check validity of the mu vector<br /> valideta::Function # check validity of the eta vector<br />end<br />
with specific distributions defined as
## utilities used in some distributions<br />logN0(x::Number) = x == 0 ? x : log(x)<br />logN0{T<:Number}(x::AbstractArray{T}) = reshape([ logN0(x[i])  i=1:numel(x) ], size(x))<br />y_log_y(y, mu) = y .* logN0(y ./ mu) # provides correct limit at y == 0<br /><br />BernoulliDist =<br /> Dist("Bernoulli",<br /> logitLink,<br /> mu > max(eps(Float64), mu .* (1.  mu)),<br /> (y, mu, wt)> 2 * wt .* (y_log_y(y, mu) + y_log_y(1.  y, 1.  mu)),<br /> (y, mu, wt)> 2. * sum(y .* log(mu) + (1.  y) .* log(1.  mu)),<br /> (y, wt)> (wt .* y + 0.5) ./ (wt + 1.),<br /> mu > all((0 < mu) & (mu < 1)),<br /> eta > true)<br />gammaDist =<br /> Dist("gamma",<br /> inverseLink,<br /> mu > mu .* mu,<br /> (y, mu, wt)> 2 * wt .* (logN0(y ./ mu)  (y  mu) ./ mu),<br /> (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))),<br /> (y, wt)> all(y > 0) ? y : error("nonpositive response values not allowed for gammaDist"),<br /> mu > all(mu > 0.),<br /> eta > all(eta > 0.))<br /><br />GaussianDist =<br /> Dist("Gaussian",<br /> identityLink,<br /> mu > ones(typeof(mu), size(mu)),<br /> (y, mu, wt)> (r = y  mu; wt .* r .* r),<br /> (y, mu, wt)> (n = length(mu); r = y  mu; n * (log(2*pi*sum(wt .* r .* r)/n) + 1) + 2  sum(log(wt))),<br /> (y, wt)> y,<br /> mu > true,<br /> eta > true)<br /><br />PoissonDist =<br /> Dist("Poisson",<br /> logLink,<br /> mu > mu,<br /> (y, mu, wt)> 2 * wt .* (y .* logN0(y ./ mu)  (y  mu)),<br /> (y, mu, wt)> 2 * sum(dpois(y, mu, true) * wt),<br /> (y, mu)> y + 0.1,<br /> mu > all(mu > 0),<br /> eta > true)<br />
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<br /> dist::Dist <br /> link::Link<br /> eta::Vector{Float64} # linear predictor<br /> mu::Vector{Float64} # mean response<br /> offset::Vector{Float64} # offset added to linear predictor (usually 0)<br /> wts::Vector{Float64} # prior weights<br /> y::Vector{Float64} # response<br />end<br />
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<br />function GlmResp(dist::Dist, link::Link, y::Vector{Float64})<br /> n = length(y)<br /> wt = ones(Float64, (n,))<br /> mu = dist.mustart(y, wt)<br /> GlmResp(dist, link, link.linkFun(mu), mu, zeros(Float64, (n,)), wt, y)<br />end<br /><br />## another outer constructor using the canonical link for the distribution<br />GlmResp(dist::Dist, y::Vector{Float64}) = GlmResp(dist, dist.canonical, y)<br />
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)<br />devResid( r::GlmResp) = r.dist.devResid(r.y, r.mu, r.wts)<br />drsum( r::GlmResp) = sum(devResid(r))<br />muEta( r::GlmResp) = r.link.muEta(r.eta)<br />sqrtWrkWt(r::GlmResp) = muEta(r) .* sqrt(r.wts ./ variance(r))<br />variance( r::GlmResp) = r.dist.variance(r.mu)<br />wrkResid( r::GlmResp) = (r.y  r.mu) ./ r.link.muEta(r.eta)<br />wrkResp( r::GlmResp) = (r.eta  r.offset) + wrkResid(r)<br />
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))<br />
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<br /> X::Matrix{Float64} # model matrix<br /> beta0::Vector{Float64} # base coefficient vector<br /> delb::Vector{Float64} # increment<br />end<br />
and an outer constructor
## outer constructor<br />predD(X::Matrix{Float64}) = (zz = zeros((size(X, 2),)); predD(X, zz, zz))<br />
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)<br /> w = sqrtWrkWt(r)<br /> wX = diagmm(w, p.X)<br /> (wX' * (w .* wrkResp(r))), (wX' * wX)<br />end<br />
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))<br />
and finally we get to the glmFit function
function glmFit(p::predD, r::GlmResp, maxIter::Uint, minStepFac::Float64, convTol::Float64)<br /> if (maxIter < 1) error("maxIter must be positive") end<br /> if (minStepFac < 0  1 < minStepFac) error("minStepFac must be in (0, 1)") end<br /> cvg = false<br /><br /> devold = typemax(Float64) # Float64 version of Inf<br /> for i=1:maxIter<br /> dev = increment(r, p)<br /> if (dev < devold)<br /> p.beta0 = p.beta0 + p.delb<br /> else<br /> error("code needed to handle the stepfactor case")<br /> end<br /> if abs((devold  dev)/dev) < convTol<br /> cvg = true<br /> break<br /> end<br /> devold = dev<br /> end<br /> if !cvg<br /> error("failure to converge in $maxIter iterations")<br /> end<br />end<br /><br />glmFit(p::predD, r::GlmResp) = glmFit(p, r, uint(30), 0.001, 1.e6)<br />
Defining two methods allows for default values for some of the arguments. Here we have a kind of allornone 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")<br /><br />## generate a Bernoulli response<br /><br />srand(123454321)<br />n = 10000<br />X = hcat(fill(1, (n,)), 3. * randn((n, 2)))<br />beta = randn((3,))<br />eta = X * beta<br />mu = logitLink.linkInv(eta)<br />y = float64(rand(n) < mu)<br />println("Coefficient vector for simulation: $(beta)")<br /><br />## establish the glmResp object<br /><br />rr = GlmResp(BernoulliDist, y)<br />pp = predD(X)<br /><br />glmFit(pp, rr)<br />println("Converged to coefficient vector: $(pp.beta0 + pp.delb)")<br />println("Deviance at estimates: $(deviance(rr))")<br /><br />## reinitialize objects for timing<br />rr = GlmResp(BernoulliDist, y)<br />pp = predD(X)<br />println("Elapsed time for fitting binary response with X $(size(X,1)) by $(size(X,2)): $(@elapsed glmFit(pp, rr)) seconds")<br />
produces
Coefficient vector for simulation: [0.54105, 0.308874, 1.02009]<br />Converged to coefficient vector: [0.551101, 0.301489, 1.02574]<br />Deviance at estimates: 6858.466827205941<br />Elapsed time for fitting binary response with X 10000 by 3: 0.09874677658081055 seconds<br />
Rbloggers.com offers daily email updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.