Algorithmic Dimensions

There are many dimensions on which we might compare a machine learning or data mining algorithm. A few of the first that come to mind are:

1) Sample complexity, convergence

How much predictive power is the algorithm able to extract from a given number of examples?

All else being equal, if algorithm A behaves the same as algorithm B with twice as many examples, we would generally prefer algorithm A. This can vary in importance depending on how scarce or expensive it is to acquire data.

2) Speed

How quickly does the algorithm run?

Obviously, a tool is more useful if it is faster. However, if that faster algorithm comes at the expense of sample complexity, one would need to measure the expense of running a program longer against the expense of gathering more data.

3) Guarantees

Does the algorithm just have good performance (along whatever dimension) in practice, or is it proven? Is it always good, or only in certain situations? How predictable is the performance?

When SVMs showed up, it looked like good practical algorithms came from good theory. These days, it seems clear that powerful intuition is also an excellent source of practical algorithms (deep learning).

Perhaps theory will some day catch up, and the algorithm with the best bounds will also perform the best. However, this is not always the case today. We often face the choice between an algorithms with good theoretical guarantees, and algorithms that work well, and are appealing otherwise. Some examples of theory versus practice are 1) Upper-Confidence Bounds versus Thomson Sampling with bandit algorithms 2) Running a convex optimization algorithm with a theoretically derived Lipschitz constant versus a smaller one that still seems to work and 3) Doing model selection via VC-dimension generalization bounds versus using K-fold cross-validation.

4) Memory usage

How much space does the algorithm need?

I work with a lot of compute-heavy applications where this is almost a wall: we don’t care about memory usage until we run out of it, after which we care a great deal. With simpler algorithms and larger datasets, this is often more nuanced, with concerns about different cache sizes.

5) Handholding

Can it be used out of the box, or does it require an expert to “tune” it for best performance?

A classic example of this is stochastic gradient descent. In principle, for a wide range of inputs, convergence is guaranteed by iteratively setting where is a noisy unbiased estimate of the gradient at and is some sequence of step-sizes that obeys the Robbins-Monro conditions [1]. However, in practice, the difference between, say, , and can be enormous, and finding these constants is a bit of a dark art.

[1] ( and .)

6) Implementability

How simple is the algorithm? How easy is it to implement?

This is quite complex and situational. These days, I’d consider an algorithm that consists of a few moderate-dimensional matrix multiplications or singular value decompositions “simple”. However, that’s due to a huge amount of effort designing reliable algorithms for these problems, and the ubiquity of easy to use libraries.

7) Amenability to parallelization

Does the algorithm lend itself to parallelization? (And what type of parallelization?)

And if it does, under what model? Map-reduce, GPU, MPI, and openMP all have different properties.

8) “Anytime-ness”

Can the algorithm be implemented in any anytime manner?

That is, does the algorithm continuously return a current “best-guess” of the answer that is refined over time in a sensible manner? This can help diagnose problems before running the full algorithm, enormously useful in debugging large systems.

Note this is distinct from being an online algorithm, which I’m not mentioning here, since it’s a mix of speed and memory properties.

9) Transparency, interpretability

Can the final result be understood? Does it give insight into how predictions are being made?

Galit Shmueli argues that “explanatory power” and “predictive power” are different dimensions. However, there are several ways in which interpretability is important even when prediction is the final goal. Firstly, the insight might convince a decision maker that the machine learning system is reliable. Second, this can be vital in practice for generalization. The world is rarely independent and identically distributed. If a domain expert can understand the predictive mechanism, they may be able to assess if this will still hold in the future, or captures something true only in the training period. Third, understanding what the predictor is doing also often yields ways to improve it. For example, the adjacent visualization of the outputs of a decision tree in two-dimensions suggests the need for non axis-aligned splits.

10) Generality

What class of problems can the algorithm address?

All else being equal, we prefer an algorithm that can, say, optimize all convex losses over one that can only fit the logistic loss. However, this often comes at a cost— a more general-purpose algorithm cannot exploit the extra structure present in a specialized problem (or, at least, has more difficulty doing so). It’s instructive how many different methods are used by LIBLINEAR depending on the particular loss and regularization constant.

11) Extendability

Does the algorithm have lots of generalizations and variants that are likely to be interesting?

More of an issue when reviewing a paper than deciding what is the final best algorithm to use once all the dust has settled.

12) Insight

Does the algorithm itself (as opposed to its results) convey insight into the problem?

Gradient descent for maximum likelihood learning of an exponential family is a good example, as it reveals the moment-matching conditions. Insight of this type, however, doesn’t suggest that one should actually *run* the algorithm.

13) Model robustness

How does the performance of the algorithm hold up to violations of its modeling assumptions?

Suppose we are fitting a simple line to some 2-D data. Obviously, all else being equal, we would prefer an algorithm that still does something reasonable when the actual dependence when the expected value of is not linear in . The example I always harp on here is the pseudolikelihood. The original paper pointed out that this will have somewhat worse sample complexity that the full likelihood, and (much!) better time complexity. Many papers seem to attribute the bad performance of the pseudolikelihood in practice to this sample complexity, when the true cause is that the likelihood does something reasonable (minimizes KL divergence) when there is model mis-specification, but the pseudolikelihood does not.

—

I often ponder is how much improvement in one dimension is enough to “matter”. Personally, I would generally consider even a small constant factor (say 5-10%) improvement in sample complexity quite important. Meanwhile, it would be rare to get exited about even, say, a factor of two improvement in running time.

What does this reflect? Firstly, generalization is the fundamental goal of data analysis, and we are likely to be willing to compromise most things if we can really predict better. Second, we instinctively distrust running times. Theory offers few tools for understanding constant factors, as these are highly architecture and implementation dependent. In principle, if one could be completely convincing that a factor of two improvement was truly there, I think this probably would be significant. (This might be true, say, if all algorithms are bottlenecked by a few matrix multiplications, and a new algorithm provably reduces the number needed.) However, this is rare.

In some places, I think constant factor skepticism can lead up astray. In reality, a factor of 30 improvement in speed is probably better than changing a complexity to . (Calculate such that .) This is particularly true when the lower-complexity algorithm has higher constant factors. As an example, I’ve always found the algorithm for projection onto the ball to be faster than the algorithm in practice.

—

Given that there are so very many dimensions, can a new algorithm really improve on all simultaneously? It seems rare to even improve on a single dimension without a corresponding cost somewhere else. There will often be a range of techniques that are Pareto optimal, depending on one’s priorities. Understanding this range is what makes having an expert around so important.

Ideally, as a community, we would be able to provide a “consumer” of machine learning an easy way to find the algorithm for them. Or, at least, we might be able to point them in the direction of possible algorithms. One admirable attempt along this line is the following table from The Elements of Statistical Learning:

(Incidentally, notice how many of the desiderata do not overlap with mine.)

Some other situations show a useful contrast. For example, take this decision tree for choosing an algorithm for unconstrained optimization, due to Dianne O’Leary:

Essentially, this amounts to the principle that one should use the least general algorithm available, so that it can exploit as much structure of the problem as possible. Though one can quibble with the details (pity the subgradient methods) it at least comes close to giving the “right” algorithm in each situation.

This doesn’t seem possible with machine learning, since there doesn’t exist a single hierarchy Rather, ML problems are a tangle of model specification, computational and architecture requirements, implementation constraints, user risk-tolerances and so on. It won’t be easy to automate away the experts. (Even ignoring the possible misalignment of incentives in that the field has the domain of automating itself.)

(This post incorporates comments from Aditya Menon and Cheng Soon Ong.)

Favorite things NIPS

I always enjoy reading conference reports, so I thought I’d mention a few papers that caught my eye.  (I welcome any corrections to my summaries of any of these.)

1. Recent Progress in the Structure of Large-Treewidth Graphs and Some Applications, by Chandra Chekuri. Not a paper, but rather an exceptionally clear tutorial on material that should probably be mandatory for anyone interested in computational complexity and graphical models. The talk isn’t online yet, but slides are available.

2. Learning Distributed Representations for Structured Output Prediction, by Vivek Srikumar and Christopher Manning. The idea is that when doing structured prediction, labels are often “similar” in some sense, and so should be encouraged to have similar weights. This is done by defining weights not over a discrete label space, but over a finite-dimensional distributed representation, with one vector for each label. Learning alternates between updating the representation of each label and regular structured learning.

3. Clamping Variables and Approximate Inference, by Adrian Weller and Tony Jebara. The paper provides a first-principles based proof that the Bethe approximation gives a lower-bound on the true partition function for attractive binary pairwise models. (A result originally proved by Ruozzi in 2012) The basic idea is that “clamping” (i.e. brute-force summing over all values of a variable, and using the Bethe approximation on the resulting model) can only increase the Bethe approximation. But, if you clamp everything, you have the exact partition function, so Bethe must be a lower-bound.  The proofs provide a lot of insight into why the results are true.

4. Making Pairwise Binary Graphical Models Attractive, by Nicholas Ruozzi and Tony Jebara. This builds on similar ideas to the above paper. There is an idea of a cover of a graph, where vertices have multiple copies, but neighborhoods look locally the same. The paper gives an appropriate cover construction such that the cover is attractive (up to a “flipping” of the variables). Thus, the Bethe partition function of the cover lower-bounds the true partition function of the cover. This is shown to be a lower-bound on the tree-reweighted partition function. As far as I can tell, Bethe on the cover may upper or lower-bound the true partition function, so Bethe-on-cover isn’t necessarily better than tree-reweighted.  Still, since tree-reweighted is typically a loose upper-bound, it is likely to be better in practice.

5. Distributed Bayesian Posterior Sampling via Moment Sharing, by Minjie Xu, Balaji Lakshminarayanan, Yee Whye Teh, Jun Zhu, and Bo Zhang. You want to do Bayesian MCMC inference using a cluster. This paper partitions the data over the nodes of the cluster, does local MCMC approximations, and uses an exponential-family approximation and expectation propagation to send information about the posterior between the nodes, thus approximating full/normal Bayesian inference. A couple days ago Andrew Gelman mentioned a similar method, considering other local approximations.

6. Mode Estimation for High Dimensional Discrete Tree Graphical Models, by Chao Chen, Han Liu, Dimitris Metaxas, and Tianqi Zhao. The goal is to find the M highest modes in a tree-structured graphical model. They find a set of necessary local conditions to be a mode, and then propose a junction-tree type algorithm operating with local mode configurations as values to find global modes.

7. Efficient Inference of Continuous Markov Random
Fields with Polynomial Potentials, by Shenlong Wang, Alexander Schwing, and Raquel Urtasun. The goal is to find a local minimum of a polynomial MRF by use of the CCCP method where one iteratively splits the objective into convex and concave parts and linearizes the concave part to get an upper-bound. I’d have expected this to be easy, but it turns out to be quite difficult to split up a polynomial into convex and concave parts. This paper gives an SDP problem to find a function to add to the objective to make it convex, yielding the desired split.

8. Hardness of parameter estimation in graphical models, by Guy Bresler, David Gamarnik, and Devavrat Shah. In principle, this is the ultimate non-surprising result: given the mean parameters (data marginals), it is computationally intractable to find the natural parameters (MRF weights) that produce these marginals. Formally given an oracle to do maximum likelihood learning in the normal way (find parameters to match a given mean vector), you could use that oracle to approximate the partition function, which is known to be hard. The really surprising thing here is that this wasn’t proved until 2014!

9. A* Sampling, by Chris Maddison, Daniel Tarlow, and Tom Minka. A method for exact sampling for low-dimensional continuous distributions, based on fundamentally different principles from existing methods. Intuitively, the algorithm simultaneously generates Gumbel-type noise for every possible (real-valued) configuration and optimizes the resulting function by branch-and-bound.

10. The workshops, as usual, were great aside from the usual annoyance of it only being possible to be in one place at one time.

11. Montreal.  Charming, friendly.

Truncated Bi-Level Optimization

In 2012, I wrote a paper that I probably should have called “truncated bi-level optimization”.  I vaguely remembered telling the reviewers I would release some code, so I’m finally getting around to it.

The idea of bilevel optimization is quite simple.  Imagine that you would like to minimize some function .  However, itself is defined through some optimization.  More formally, suppose we would like to solve

Or, equivalently,

where is defined as .  This seems a little bit obscure at first, but actually comes up in several different ways in machine learning and related fields.

Hyper-parameter learning

The first example would be in learning hyperparameters, such as regularization constants.  Inevitably in machine learning, one fits parameters parameters to optimize some tradeoff between the quality of a fit to training data and a regularization function being small.  Traditionally, the regularization constant is selected by optimizing on a training dataset with a variety of values, and then picking the one that performs best on a held-out dataset.  However, if there are a large number of regularization parameters, a high-dimensional grid-search will not be practical.  In the notation above, suppose that is a vector of regularization constants, and that are training parameters.  Let,   be the regularized empirical risk on a training dataset, and let be how well the parameters perform on some held-out validation dataset.

Energy-based models

Another example (and the one suggesting the notation) is an energy-based model.  Suppose that we have some “energy” function which measures how well an output fits to an input .  The energy is parametrized by .  For a given training input/output pair , we might have that measures how how the predicted output compares to the true output , where .

Computing the gradient exactly

Even if we just have the modest goal of following the gradient of to a local minimum, even computing the gradient is not so simple.  Clearly, even to evaluate requires solving the “inner” minimization of . It turns out that one can compute through first solving the inner minimization, and then solving a linear system.

1. Input
2. Solve .
3. Compute:
4.    (a) the loss
5.    (b) the gradient
6.    (c) the Hessian
7. Solve the linear system .
8. Compute the parameter gradient
9. Return and .

This looks a bit nasty, since we need to compute second-derivative matrices of .  In fact, as long as one has a routine to compute and , this can be avoided through Efficient Matrix-vector products.  This is essentially proposed by Do, Foo, and Ng in “Efficient multiple hyperparameter learning for log-linear models”.

Overall, this is a decent approach, but it can be quite slow, simply because one must solve an “inner” optimization in order to compute each gradient of the “outer” optimization.  Often, the inner-optimization needs to be solved to very high accuracy in order to estimate a gradient accurately enough to reduce — higher accuracy than is needed when one is simply using the predicted value itself.

Truncated optimization

To get around this expense, a fairly obvious idea is to re-define the problem.  The slowness of exactly computing the gradient stems from needing to exactly solve the inner optimization.  Hence, perhaps we re-define the problem such that an inexact solve of the inner problem nevertheless yields an “exact” gradient?

Re-define the problem as solving

,

where denotes an approximate solve of the inner optimization.  In order for this to work, must be defined in such a way that is a continuous function of .  With standard optimization methods such as gradient descent or BFGS, this can be achieved by assuming there are a fixed number of iterations applied, with a fixed step-size.  Since each iteration of these algorithms is continuous, this clearly defines as a continuous function.  Thus, in principle, it could be optimized efficiently through automatic differentiation of the code that optimizes .  That’s fine in principle, but often inconvenient in practice.

It turns out, however, that one can derive “backpropagating” versions of algorithms like gradient descent, that take as input only a procedure to compute along with it’s first derivatives.  These algorithms can then produce the gradient of in the same time as automatic differentiation.

Back Gradient-Descent

If the inner-optimization is gradient descent for steps with a step-size of , the algorithm to compute the loss is simple:

1. Input
2. For
3.     (a)
4. Return

How to compute the gradient of this quantity?  The following algorithm does the trick.

1. 
2. 
3. For
4.     (a)
5.     (b)
6. Return .

Similar algorithms can be derived for the heavy-ball algorithm (with a little more complexity) and limited memory BFGS (with a lot more complexity).

Code

So, finally, here is the code, and I’ll give a simple example of how to use it. There are just four simple files:

1. back_gd.m
2. back_hb.m
3. back_lbfgs.m
4. demo.m

I think the meanings of this are pretty straightforward, so I’ll just quickly step through the demo here.  I’ll start off by grabbing taking one of Matlab’s built-in datasets (on cities) so that we are trying to predict a measure of crime from measures of climate, housing, health, transportation, arts, recreation, and economy, as well as a constant.  There are 329 data, total, which I split into a training set of size 40, a validation set of size 160, and a test set of size 129.

load cities ratings
X = ratings;
for i=1:size(X,2)
    X(:,i) = X(:,i) - mean(X(:,i));
    X(:,i) = X(:,i) / std( X(:,i));
end

% predict crime from climate, housing, health, trans, edu, arts, rec, econ,
Y = X(:,4);
X = [X(:,[1:3 5:9]) 1+0*X(:,1)];

p = randperm(length(Y));
X = X(p,:);
Y = Y(p);

whotrain = 1:50;
whoval   = 51:200;
whotest  = 201:329;

Xtrain = X(whotrain,:);
Xval   = X(whoval  ,:);
Xtest  = X(whotest ,:);
Ytrain = Y(whotrain);
Yval   = Y(whoval  );
Ytest  = Y(whotest );

Next, I’ll set up some simple constants that will be used later on, and define the optimization parameters for minFunc, that I will be using for the outer optimization. In particular, here I choose the inner optimization to use 20 iterations.

opt_iters = 20;
ndims = size(Xtrain,2); w0    = zeros(ndims,1); ndata = size(Xtrain,1); ndata_val = size(Xval,1);
options = []; options.Method  = 'gd'; options.LS      = 1; options.MaxIter = 100;

Now, I’ll define the training risk function ( in the notation above). The computes the risk with a regularization constant of , as well as derivatives. I’ll also define the validation risk ( in the notation above).

function [R dRdw dRdloga] = training_risk(w,loga)
        a         = exp(loga);
        R         = sum( (Xtrain*w - Ytrain).^2 )/ndata + a*sum(w.^2);
        dRdw      = 2*Xtrain'*(Xtrain*w-Ytrain)  /ndata + 2*a*w;
        dRda      = sum(w.^2);
        dRdloga = dRda*a;
    end

    function [R g] = validation_risk(w)
        R = sum( (Xval*w - Yval).^2 ) / ndata_val;
        g = 2*Xval'*(Xval*w-Yval)     / ndata_val;
    end

Now, before going any further, let’s do a traditional sweep through regularization constants to see what that looks like.

LAMBDA = -5:.25:2;
    VAL_RISK = 0*LAMBDA;
    for i=1:length(LAMBDA)
        VAL_RISK(i) = back_lbfgs(@training_risk,@validation_risk,w0,LAMBDA(i),opt_iters);
    end

Plotting, we get the following:

This is a reasonable looking curve. Instead, let’s ask the algorithm to find the constant by gradient descent.

eval = @(loga) back_lbfgs(@training_risk,@validation_risk,w0,loga,opt_iters);        
loga = 0;
[loga fval] = minFunc(eval,loga,options);

Running the optimization, we see:

Iteration   FunEvals     Step Length    Function Val        Opt Cond
         1          2     1.00000e+00     8.74176e-01     3.71997e-03
         2          3     1.00000e+00     8.73910e-01     1.86453e-04
         3          4     1.00000e+00     8.73909e-01     1.06619e-05
         4          5     1.00000e+00     8.73909e-01     3.60499e-08

This leads to a regularizer of the form:

.

We can plot this on the graph, and see it matches the result of cross-validation.

If we actually compute the test-set error, this is 0.708217.

OK, let’s be a little bit more adventurous, and use a third-order regularizer. This is done like so:

    function [R dRdw dRdloga] = training_risk2(w,loga)
        a         = exp(loga(1));
        b         = exp(loga(2));
        R         = sum( (Xtrain*w - Ytrain).^2 ) / ndata + a*sum(abs(w).^2)     + b*sum(abs(w).^3);
        dRdw      = 2*Xtrain'*(Xtrain*w-Ytrain)   / ndata + a*abs(w).^1.*sign(w) + b*abs(w).^2.*sign(w);
        dRda      = sum(abs(w).^2);
        dRdb      = sum(abs(w).^3);
        dRdloga = [dRda*a; dRdb*b];
    end

Running the optimization, we see:

Iteration   FunEvals     Step Length    Function Val        Opt Cond
         1          2     1.00000e+00     8.74445e-01     1.17262e-02
         2          3     1.00000e+00     8.73685e-01     3.21956e-03
         3          4     1.00000e+00     8.73608e-01     1.41744e-03
         4          5     1.00000e+00     8.73598e-01     8.20040e-04
         5          6     1.00000e+00     8.73567e-01     1.39830e-03
         6          7     1.00000e+00     8.73513e-01     2.52994e-03
         7          8     1.00000e+00     8.73471e-01     1.77157e-03
...
        23         28     1.00000e+00     8.70741e-01     8.42628e-06

With a final regularizer of the form

and a hardly-improved test error of 0.679155.

Finally, let’s fit a fourth-order polynomial.

    function [R dRdw dRdloga] = training_risk3(w,loga)
        a         = exp(loga);
        R         = sum( (Xtrain*w - Ytrain).^2 ) / ndata;
        dRdw      = 2*Xtrain'*(Xtrain*w-Ytrain)   / ndata;
        for ii=1:length(a)
            b=ii+1;
            R    = R    + a(ii)*sum(abs(w).^b);
            dRdw = dRdw + a(ii)*abs(w).^(b-1).*sign(w);
            dRda(ii,1) = sum(abs(w).^b);
        end
        dRdloga = dRda.*a;
    end

    eval = @(loga) back_lbfgs(@training_risk3,@validation_risk,w0,loga,opt_iters);
    loga = [0; 0; 0; 0];
    loga = minFunc(eval,loga,options);

Running the optimization, we see

Iteration   FunEvals     Step Length    Function Val        Opt Cond
         1          2     1.00000e+00     8.73982e-01     1.02060e-02
         2          3     1.00000e+00     8.73524e-01     3.30445e-03
         3          4     1.00000e+00     8.73447e-01     1.74779e-03
         4          5     1.00000e+00     8.73435e-01     1.43655e-03
         5          6     1.00000e+00     8.73345e-01     4.98340e-03
         6          7     1.00000e+00     8.73295e-01     1.85535e-03
         7          8     1.00000e+00     8.73231e-01     1.81136e-03
...
        38         41     1.00000e+00     8.69758e-01     3.44848e-06

Yielding the (strange) regularizer

and a final test-error of 0.67187.

Reducing Sigmoid computations by (at least) 88.0797077977882%

A classic implementation issue in machine learning is reducing the cost of computing the sigmoid function

.

Specifically, it is common to profile your code and discover that 90% of the time is spent computing the in that function.  This comes up often in neural networks, as well as in various probabilistic architectures, such as sampling from Ising models or Boltzmann machines.  There are quite a few classic approximations to the function, using simple polynomials, etc. that can be used in neural networks.

Today, however, I was faced with a sampling problem involving the repeated use of the sigmoid function, and I noticed a simple trick that could reduce the number of sigmoids by about 88% without introducing any approximation.  The particular details of the situation aren’t interesting, but I repeatedly needed to do something like the following:

1. Input
2. Compute a random number
3. If
4.   Output
5. Else
6.     Output

Now, let’s assume for simplicity that is positive.  (Otherwise, sample using and then switch the sign of the output.)  There are two observations to make:

1. If is large, then you are likely to output
2. Otherwise, there are easy upper and lower bounds on the probability of outputting

This leads to the following algorithm:

1. Input
2. Compute a random number
3. If
4.     If or
5.         Output
6.     Else
7.         Output
8. Else
9.     If
10.         Output
11.     Else if or
12.         Output
13.     Else
14.         Output

The idea is as follows:

1. If , then we can lower-bound the probability of outputting +1 by a pre-computed value of , and short-circuit the computation in many cases.
2. If , then we can upper bound the sigmoid function by .
3. If , then we can also lower bound by respectively.  (This constant was found numerically).

The three cases are illustrated in the following figure, where the input is on the x-axis, and the random number is on the y-axis.

Since, for all at least a fraction of the numbers will be short-circuited, sigmoid calls will be reduced appropriately.  If is often near zero, you will do even better.

Obviously, you can take this farther by adding more cases, which may or may not be helpful, depending on the architecture, and the cost of branching vs. the cost of computing an .

Windows Binaries for CRF toolbox

I’m happy to report that, thanks to Alexei Skurikhin, Windows binaries are now available for all the functions in my CRF toolbox.  Hopefully this makes for an easier out of the box experience from Windows folks who, based on my email, have struggled to compile things in the past.

Etymology of Latin words/abbreviations

I often find the academic convention of using Latin words and abbreviations a bit silly.  (Though I’m as guilty as anyone.) Today, I decided to look into a couple common ones, to see if there is any good reason not to use the ordinary English equivalent.

• “e.g.” is an abbreviation of the Latin “exemplī grātiā”, meaning “for example”.  Here, “exemplī” is not hard to guess and “grātiā” means “for the sake of”.
• “i.e.” is an abbreviation of the Latin “id est”, meaning “that is”.
• “via”.  This is not an abbreviation, but was the Latin word for a road.  The plural is “viae”, which you should probably work into your next paper submission.  (I think the usage would be “Linear convergence and entropic sparsity control viae quasi-stochasticity and entropy regulators”?)
• “et al.”.  Here “et” is Latin for “and” while “al.” is an abreviation for “alii”, meaning “others”.  I seem to recall a major statistics journal maintaining the noble tradition of actually writing “and others”, but can’t recall which one.

I see no particular reason we haven’t adopted English words for all these conventions, but in the context of the many local optima society is in… this hardly seems worth changing.

Julia, Matlab, and C

Julia is a new language in the same arena as Matlab or R. I’ve had failed attempts to quit the Matlab addiction in the past, making me generally quite conservative about new platforms. However, I’ve recently been particularly annoyed by Matlab’s slow speed, evil license manager errors, restrictions on parallel processes, C++ .mex file pain, etc., and so I decided to check it out. It seems inevitable that Matlab will eventually displaced by something. The question is: is that something Julia?

“We want a language that’s open source, with a liberal license. We want the speed of C with the dynamism of Ruby. We want a language that’s homoiconic, with true macros like Lisp, but with obvious, familiar mathematical notation like Matlab. We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell. Something that is dirt simple to learn, yet keeps the most serious hackers happy. We want it interactive and we want it compiled.”

Essentially, the goal seems to be a faster, freer Matlab that treats users like adults (macros!) and doesn’t require writing any .mex files in C++ or Fortan. Sounds too good to be true? I decided to try it out. My comparisons are to Matlab (R2012a) and C (gcc 4.2 with -O2 and -fno-builtin to prevent compile-time computations) all on a recent MacBook Pro (2.3 GHz Intel Core i7 w/ 8GB RAM).

Installation was trivial: I just grabbed a pre-compiled binary and started it up.

I deliberately used naive algorithms, since I am just testing raw speed. It should be a fair comparison, as long as the algorithm is constant. Please let me know about any bugs, though.

Click here to skip to the results.

First benchmark: Fibonnaci

% Matlab
function f=fib(n)
  if n <= 2
    f=1.0;
  else
    f=fib(n-1)+fib(n-2);
  end
end

// C
double fib(int n){
  if(n<=2)
    return(1.0);
  else
    return(fib(n-2)+fib(n-1));
}

% julia
function fib(n)
  if n <= 2
    1.0
  else
    fib(n-1)+fib(n-2);
  end
end

Clarity is basically a tie. Running them for n=30 we get:

time in matlab (fib): 14.344231
time in c (fib):       0.005887
time in julia (fib):   0.237832

Second benchmark: Matrix Multiplication

This is a test of the naive O(N^3) matrix multiplication algorithm.

% matlab
function C=mmult(A,B,C)
  [M,N] = size(A);
  for i=1:M
    for j=1:M
      for k=1:M
        C(i,j) = C(i,j) + A(i,k)*B(k,j);
      end
    end
  end
end

// C
#define M 500
void mmult(double A[M][M],double B[M][M], double C[M][M]){
  //double C[M][M];
  int i,j,k;
  for(i=0; i<M; i++)
    for(j=0; j<M; j++){
      C[i][j] = 0;
      for(k=0; k<M; k++)
	C[i][j] += A[i][k]*B[k][j];
    }
}

# julia
function mmult(A,B)
  (M,N) = size(A);
  C = zeros(M,M);
  for i=1:M
    for j=1:M
      for k=1:M
        C[i,j] += A[i,k]*B[k,j];
      end
    end
  end
  C;
end

Here, I think that Matlab and Julia and a bit clearer, and Julia wins though the wonders of having “+=”. The timing results on 500×500 matrices are:

time in matlab (matmult): 1.229571
time in c (matmult):      0.157658 
time in julia (matmult):  0.5029549

Third Benchmark: numerical quadrature

Here, we attempt to calculate the integral by numerical quadrature, using a simple midpoint rule with computations at points.

% matlab
function val=numquad(lb,ub,npoints)
  val = 0.0;
  for x=lb:(ub-lb)/npoints:ub
    val = val + sin(x)/npoints;
  end
end

// C
double numquad(double lb,double ub,int npoints){
  double val = 0.0;
  int i;
  for(i=0; i<=npoints; i++){
    double x = lb + (ub-lb)*i/npoints;
    val += sin(x)/npoints;
  }
  return(val);
}

# julia
function numquad(lb,ub,npoints)
  val = 0.0
  for x=lb:(ub-lb)/npoints:ub
    val += sin(x)/npoints
  end
  val
end

The timings are:

time in matlab (numquad): 0.446151
time in c (numquad):      0.167112 
time in julia (numquad):  0.256597

Fourth Benchmark: Belief Propagation

Finally, I decided to try a little algorithm similar to what I actually tend to implement for my research. Roughly speaking, Belief Propagation is a repeated sequence of matrix multiplications, followed by normalization.

% matlab
function x=beliefprop(A,x,N)
  for i=1:N
    x = A*x;
    x = x/sum(x);
  end
end

// C
void beliefprop(double A[25][25], double x[25], int N){
  int i, n, j;
  double x2[25];
  for(n=0; n<N; n++){
    for(i=0; i<25; i++){
      x2[i]=0;
      for(j=0; j<25; j++)
	x2[i] += A[i][j]*x[j];
    }   
    for(i=0; i<25; i++)
      x[i]=x2[i];
    double mysum = 0;
    for(i=0; i<25; i++)
      mysum += x[i];
    for(i=0; i<25; i++)
      x[i] /= mysum;
  }
  return;
}

% julia
function beliefprop(A,x,N)
  for i=1:N
    x = A*x;
    x /= sum(x);
  end
  x
end

Here, I think we can agree that Matlab and Julia are clearer. (Please don’t make fun of me for hardcoding the 25 dimensions in C.) Using a matrix package for C would probably improve clarity, but perhaps also slow things down. The results are:

time in matlab (beliefprop): 0.627478
time in c (beliefprop):      0.074355 
time in julia (beliefprop):  0.376427

Fifth Benchmark: BP in log-space

In practice, Belief Propagation is often implemented in log-space (to help avoid numerical under/over-flow.). To simulate an algorithm like this, I tried changing to propagation to take an exponent before multiplication, and a logarithm before storage.

% matlab
function x=beliefprop2(A,x,N)
  for i=1:N
    x = log(A*exp(x));
    x = x - log(sum(exp(x)));
  end
end

// C
void beliefprop2(double A[25][25], double x[25], int N){
  int i, n, j;
  double x2[25];
  for(n=0; n<N; n++){
    for(i=0; i<25; i++){
      x2[i]=0;
      for(j=0; j<25; j++)
	x2[i] += A[i][j]*exp(x[j]);
    }   
    for(i=0; i<25; i++)
      x[i]=log(x2[i]);
    double mysum = 0;
    for(i=0; i<25; i++)
      mysum += exp(x[i]);
    double mynorm = log(mysum);
    for(i=0; i<25; i++)
      x[i] -= mynorm;
  }
  return;
}

# julia
function beliefprop2(A,x,N)
  for i=1:N
    x = log(A*exp(x));
    x -= log(sum(exp(x)));
  end
  x
end

Life is too short to write C code like that when not necessary. But how about the speed, you ask?

time in matlab (beliefprop2): 0.662761
time in c (beliefprop2):      0.657620 
time in julia (beliefprop2):  0.530220

Sixth Benchmark: Markov Chain Monte Carlo

Here, I implement a simple Metropolis algorithm. For no particular reason, I use the two-dimensional distribution:

% matlab
function mcmc(x,N)
f = @(x) exp(sin(x(1)*5) - x(1)^2 - x(2)^2);
p = f(x);
for n=1:N
  x2 = x + .01*randn(size(x));
  p2 = f(x2);
  if rand < p2/p
    x = x2;
    p = p2;
  end
end
end

// C
double f(double *x){
  return exp(sin(x[0]*5) - x[0]*x[0] - x[1]*x[1]);
}
#define pi 3.141592653589793
void mcmc(double *x,int N){
  double p = f(x);
  int n;
  double x2[2];
  for(n=0; n<N; n++){
    // run Box_Muller to get 2 normal random variables
    double U1 = ((double)rand())/RAND_MAX;
    double U2 = ((double)rand())/RAND_MAX;
    double R1 = sqrt(-2*log(U1))*cos(2*pi*U2);
    double R2 = sqrt(-2*log(U1))*sin(2*pi*U2);
    x2[0] = x[0] + .01*R1;
    x2[1] = x[1] + .01*R2;
    double p2 = f(x2);
    if(((double)rand())/RAND_MAX< p2/p){
      x[0] = x2[0];
      x[1] = x2[1];
      p = p2;
    }
  }
}

% julia
function mcmc(x,N)
f(x) = exp(sin(x[1]*5) - x[1]^2 - x[2]^2);
p = f(x);
for n=1:N
  x2 = x + .01*randn(size(x));
  p2 = f(x2);
  if rand() < p2/p
    x = x2;
    p = p2;
  end
end
end

Again, I think that C is far less clear than Matlab or Julia. The timings are:

time in matlab (mcmc): 7.747716
time in c (mcmc):      0.150776 
time in julia (mcmc):  0.479628

Table

All times are in seconds. (Lower is better.)

           Matlab     C     Julia
fib        14.344   0.005   0.237
matmult     1.229   0.157   0.502
numquad     0.446   0.167   0.256
bp          0.627   0.074   0.376
bp2         0.662   0.657   0.530
mcmc        7.747   0.150   0.479

Conclusions

I’m sure all these programs can be sped up. In particular, I’d bet that an expert could optimize the C code to beat Julia on bp2 and mcmc. These are a test of “how fast can Justin Domke make these programs”, not the intrinsic capabilities of the languages. That said, Julia allows for optional type declarations. I did experiment with these but found absolutely no speed improvement. (Which is a good or a bad thing, depending on how you look at life.)

Another surprise to me was how often Matlab’s JIT managed a speed within a reasonable factor of C. (Except when it didn’t…)

The main thing that at Matlab programmer will miss in Julia is undoubtedly plotting. The Julia designers seem to understand the importance of this (“non-negotiable”). If Julia equalled Matlab’s plotting facilities, Matlab would be in real trouble!

Overall, I think that the killer features of freedom, kinda-sorta-C-like speed, and ease of use make Julia more likely as a Matlab-killer than other projects such as R, Sage, Octave, Scipy, etc. (Not to say that those projects have not succeeded in other ways!) Though Julia’s designers also seem to be targeting current R users, my guess is that they will have more success with Matlab folks in the short term, since most Matlab functionality (other than plotting) already exists, while reproducing R’s statistical libraries will be quite difficult. I also think that Julia would be very attractive to current users of languages like Lush. Just to never write another .mex file, I’ll very seriously consider Julia for new projects. Other benefits such as macros, better parallelism support are just bonuses. As Julia continues to develop, it will become yet more attractive.

There was an interesting discussion on Lambda the Ultimate about Julia back when it was announced

Posted in Uncategorized | Tagged c++, efficiency, julia, matlab | 32 Comments