Expensive lessons in Python performance tuning

Over the past year I've done a lot of number crunching in Python with the wonderful ipython/numpy/pandas/sklearn stack and using picloud to parallelize my computations. Picloud is an amazing service which lets you run your code in Amazon's behemoth data centers while hiding most of the configuration headache. Their computational model is essentially a parallel map: you submit a list of inputs and specify what function to map over them and picloud automagically launches that function on multiple machines. This is less powerful than full-blown MapReduce (no parallel reductions) but, luckily, a lot of the computations machine learning people want to perform are trivial to parallelize. These include:

• parameter searches (parallelize over distinct parameters)
• k-fold cross validation (parallelize over the folds)
• training random forests (trees don't need to know about each other)
• random starts for locally converging hill climbers like k-means, EM, etc..

As long as you have multiple things to do in parallel and the work for each item greatly exceeds the communication time, then picloud can be a miracle. However, they do charge a premium over the standard AWS instance prices and even Amazon's baseline prices can add up quickly if you're using powerful instance types. When you're running on someone else's hardware the inefficiencies of your code translate very directly into added costs. One month, I managed to rack up an extremely high bill from an embarrassingly small amount of actual work because (1) my Python code created millions of bloated objects in memory, forcing me to use the more expensive high-memory instance type, and (2) I had overlooked some very simple optimizations that later cut my runtime by an order of magnitude.  

So, I thought I might save someone else time and money by sharing some tips which made my Python code more efficient without sacrificing flexibility or maintainability. 

1. Profile everything! Profile mercilessly and relentlessly. Let no assumption about efficiency pass without being put under the cold and uncaring gaze of a profiler. This applies to both runtime (use ipython's prun command) and memory consumption (which you can inspect with the over-engineered guppy). Don't "optimize" a damn thing until you know for sure it's a top bottleneck. 
2. Be suspicious of loops (but have faith in Python's core operations). If your bottleneck is a loop performing numerical computations you can almost always speed it up by using some equivalent vectorized NumPy operations. On the other hand, once the obviously vectorizable stuff is off-loaded to library code, then I'm often surprised just how efficient Python's operations can be. I can build lists of startling length, index into dictionaries all day long, construct millions of simple objects, and it all adds up to a small fraction of my total runtime. 
3. Simplest is often fastest. I spent a whole day concocting elaborate buffering and pre-fetching schemes to traverse the lines of a file as quickly as possible without reading it all into memory. I neglected, however, to compare with the trivial solution of simply writing "for line in file:", which of course turned out to be faster and more efficient. Similarly, I spent a while tweaking a CSV parsers before finding that using "string.split" with the standard parsing routines int and float was just as fast. 
4. If you're creating millions of objects, make them namedtuples. In one feature extraction pipeline, I take a multi-million line CSV and parse each line into a Python object. My ever-frustrated internal micro-optimizer screams at the obvious inefficiency. And, sadly, he's right. Python objects are extremely heavy-weight compared to the simple struct we would create in C. Luckily, Guido and Friends foresaw the "zillions of simple objects" use-case and thus created namedtuples. They have named fields like objects, but are laid out contiguously and take up a lot less space. Alternatively, if you're using an object for more than just data storage and want to use non-trivial methods, you can lay out a full-blown object contiguously by filling out the slots attribute. How efficient are they? As a totally unscientific demonstration I constructed a million little vectors with field names "x", "y", and "z". When implemented as ordinary Python objects they took up 400MB of memory, whereas their namedtuple cousins only required 114MB. Still not nearly as compact as they could be, but a dramatic improvement nonetheless! 
5. Don't search linearly through a sorted array, use binary search (duh). On the one hand, this is CS 101 kid's stuff. On the other, I didn't know about the bisect library and had code littered with np.nonzero(sorted<=x)[0]. Once I switched to using bisect_left/bisect_right I saw a huge performance improvement. edit: In the comments, Peter pointed out that NumPy implements a faster binary search called searchsorted; you (and I) should probably use that instead!.
6. Numpy's "fancy indexing" is slow. If you index into a NumPy array with a range then you get a very lightweight view over the original array; this kind of indexing doesn't copy any data and is almost free. If, however, you index with a boolean vector or a sequence of integers ("fancy indexing"), then you actually incur a significant cost. If your code does a lot of fancy indexing then you should either switch to using the take method or rewrite your indexing as an explicit loop in native code (see #8).
7. Wes McKinney is a genius. If you're implementing anything Wes McKinney has already put in his library pandas, just stop. His code is faster, more robust, and more likely to be correct than anything you're about to write. Want rolling window aggregators? Use pandas. Need to handle missing data? Use pandas. Are you writing some sort of unbelievably ugly hack which attempts to implement joins and group-by's over NumPy arrays but actually manages to spend 3 hours computing a subtly incorrect result? (I have done this). Jesus Christ, just stop and use pandas.
8. When all else fails, try weave. Yes, weave is unmaintained, ugly, and hard to debug. But, if you find yourself thinking "if only I could just inline 5 lines of C right in the middle of this Python code", then weave lets you actually do that. If, however, that 5 lines turns into 50 then maybe it's time to learn about Cython.

Also, a few last words specific to machine learning: climb the complexity ladder carefully. It's tempting to throw a fancy tool at some new data, and doubly tempting to try every possible combination of preprocessing steps and hyperparameters. Brute force search and a muscular learning algorithm let you sidestep the need for thinking, experimentation, or understanding. Of course, in my experience, this rarely works out and ends up costing quite a bit of money.

So, do the more prudent mature thing and start simple. Collect simple statistics, plot each feature, and try to understand the data a little before you launch 100,000 jobs on picloud. Even then, start with simpler learners like  SGD Logistic Regression or small Random Forests. These will let you explore your choices of preprocessing steps and target functions without eating up too many resources. Once you're more comfortable with the problem, go ahead and train a 1000 layer deep belief network (if that's your thing).

Postscript: After almost a year, this article seems out of date with regard to the various ways you can generate native code. Particularly, before you go reaching for weave, Cython, or (gasp!) even C itself, it might be worthwhile to try selective just-in-time compilation. Both Numba and my own Parakeet will generate efficient native implementations for array-oriented programs. Copperhead is a little more limited semantically (it's really a subset of Haskell masquerading as Python), but if you can squeeze your algorithm into Copperhead's model then you'll get a blazing fast GPU kernel as a reward. Alternatively, you can skip the "pretending to be Python" charade and build data-parallel expression trees directly using Theano. It can be a little rough around the edges but, like Copperhead, if you can get it working then the (many orders of magnitude faster) GPU kernel performance is well worth the trouble.

Which classifiers are fast enough for exploring medium-sized data?

When I'm constructing a classification pipeline, there are usually a great number of choices to be made about feature selection/transformation/extraction. I often find myself thinking things like:

• "Is number of burps a useful feature?" 
• "Should I take the log of lollipops per hectare?" 
• "I don't know what structure this time series actually has, so fuck it, I'll try wavelets..."
• "Andrew Ng told me learning a sparse overcomplete feature representation might work..."

etc..

While you can choose and manipulate your features by divination, tea leaf reading, or good old-fashioned gut-following, it's better to actually check whether what you're doing improves performance on a validation set. The trouble is that repeatedly training a classifier can become cumbersome once you leave the cozy world of Small Data. Even worse, if you want to use some sort of automatic feature selection algorithm (which requires training hundreds or thousands of classifiers), you might find yourself waiting for a very very long time.

So which classifiers can I realistically use for exploring medium-sized data sets (those small enough to fit in memory, but too big for something silly like constructing a kernel matrix)? I tested a variety of scikit-learn's classifiers to see how they compare.

Data 

I generated 50,000 samples from 10 distinct 500-dimensional gaussians (with random means and random diagonal covariance matrices), and randomly assigned each sample to either a test or training set. This yielded a data set with 10 classes, 500 dimensions, 24,853 training samples and 25,147 test samples.

In addition to checking the runtime and accuracy of each classifier, I also wanted to know how well they cope with noisy/irrelevant features. To this end, I made a copy of the data with 258 of the 500 features replaced by uniform noise.

Results

Algorithm	Training/Test Time	Training/Test Accuracy	Corrupt Feature Accuracy
Linear Discriminant Analysis	9.16s / 0.23s	51.64% / 46.19%	37.39% / 30.96%
Quadratic Discriminant Analysis	14.78s / 34.27s	100% / 100%	100% / 100%
Linear SVM (trained with liblinear)	449.9s / 0.49s	38.81% / 34.95%	25.26% / 21.27%
Linear SVM (trained with SGD)	3.8s / 0.13s	50.13% / 43.68%	24.4% / 21.38%
L1-penalized SVM (trained with SGD)	27.49 s / 0.13s	50.29% / 43.39%	33.81% / 26.5%
Logistic Regression	19.74s / 0.49s	51.69% / 45.56%	37.24% / 30.37%
Gaussian Naive Bayes	0.73s / 5.14s	100% / 100%	100% / 100%
Decision Tree (CART)	20.4s / 0.1s	100% / 61%	100% / 58.9%
Random Forest (10 trees)	49.45s / 1.01s	99.96% / 84.32%	99.93% / 74.2%
Random Forest (20 trees)	98.62s / 2.02s	100% / 95.2%	100% / 88.8%
Random Forest (40 trees)	197.56s / 4.03s	100% / 99.2%	100% / 96.75%
Random Forest (80 trees)	393.5s / 8.08s	100% / 99.96%	100% / 99.26%
Random Forest (160 trees)	789.2s / 16.01s	100% / 100%	100% / 99.82%
Extremely Randomized Tree	6.58s / 0.1s	100% / 35.4.6%	100% / 29%
Extra Trees (10 trees)	57.94s / 1.04s	100% / 72.6%	100% / 55.4%
Extra Trees (20 trees)	110.84s / 2.12s	100% / 90.1%	100% / 72.56%
Extra Trees (40 trees)	221.43s / 4.23s	100% / 98.16%	100% / 88.99%
Extra Trees (80 trees)	441.7s / 8.46s	100% / 99.87%	100% / 88.99%
Extra Trees (160 trees)	878.2s / 16.84s	100% / 100%	100% / 99.4%
Gradient Boosting (4 trees per class)	59.93s / 0.23s	82.9% / 80.59%	75.22% / 72.91%
Gradient Boosting (8 trees per class)	118.66s / 0.36s	93.85% / 92.59%	89.04% / 87.19%
Gradient Boosting (16 trees per class)	234.78s / 0.61s	98.75% / 97.94%	96.28% / 94.72%
Gradient Boosting (32 trees per class)	464.89s / 1.11s	99.9% / 99.72%	99.2% / 98.52%

Observations

• If you know the process by which your data was generated then consider your work done. In this case, the data was sampled from 10 Gaussian distributions and (unsurprisingly) QDA and Gaussian Naive Bayes both do very well. If I had used Gaussians with non-zero off-diagonal covariances (interactions between features), then I would expect Gaussian Naive Bayes to do worse but QDA would still capture the best possible class boundaries. 
• SVM and Logistic Regression are very sensitive to the actual margin between classes. I retrained the linear learners with better separated class means and found their test accuracy went up to ~99% (whereas tree-based algorithms only improved slightly). 
• Stochastic gradient descent for L2-regularized learners is blazing fast. Kudos to Peter Prettenhofer for implementing something which makes every other linear learner look like a slug.
• Though still fast compared to the glacial alternative of using a convex solver, SGD for L1-penalized learning is noticeably slower. 
• As advertised, a single Extremely Randomized Tree is terrible by itself. An Extra Tree ensemble, however, seems competitive with Random Forests for clean data, but shows worse degradation when trained on corrupted features. This makes sense, since a Random Forest's tree-growing algorithm will get a chance to skip useless features (and an Extremely Randomized Tree uses all features with equal probability). 
• Only Naive Bayes and SGD-SVM are cheap enough for use with greedy feature selection. 

Next Time

• Since the base learners of a Random Forest took ~5 seconds to train, can random forests for 2 or 3 trees give good enough accuracy to be useful?
• What if we use smaller bootstrap samples (train each tree on 10% or even 1% of the data), how badly will accuracy suffer?
• Do any of the "expensive" algorithms become cheaper on a differently shaped data set (more rows, fewer features)?
• The SGD learners above performed about ~1 million weight updates (the suggested heuristic from the sklearn page). Presumably I could get a 10X performance gain by using only 10^5 updates. How badly will this impact generalization/test error?
• Simulated Gaussian data is totally unrepresentative! How do these algorithms behave on real data? I'm going to try them out on some intra-day currency time series and some other data we found lying around on kaggle.