Do you know that convolution operation is implemented in deep learning systems via matrix multiplication?

Here is a detailed description of how it is done by chunking the data
http://petewarden.com/2015/04/20/why-gemm-is-at-the-heart-of-deep-learning/

Theano library, by the way, provides many possibilities for implementation of convolution, among which CPUCorrMM (this isn't the main way however)
http://deeplearning.net/software/theano/library/tensor/nnet/conv.html

A Programming Language

Since I'm very interested in numpy and vectorization, I became curious about when how this approach appeared. I find this notion very convenient and very producing, while it leaves the space for optimizations inside.

Obviously, a predecessor of numpy is MATLAB, which appeared a long time ago - in the late 70's, while the first release happened more than 30 years ago - in 1984! Many things were inherited from this language (not the multiCPU/GPU backend, unfortunately, but I hope that things will change oce).

But more interesting moment that vector-based syntax was not fortran, as I thought earlier. The latter got it's operations over vectors only in Fortran'90 (which was much later). 

The real predecessor of MATLAB was APL (A Programming Language, not named after Apple company), this language was quite complicated and short,  and mostly this was looking not like a code, but as numerous formulas. Apart from this,  APL used many specific symbols and combined operators. For instance,  +\ seq is sum of elements in sequence.

Wikipedia claims that this returns list of prime numbers in $1, 2, ..R$. 

(~R∊R∘.×R)/R←1↓ιR

Now you understand why this language isn't popular today :)

From that moment we learnt that many things like map, where, sum can be written using words, and this will be much more clear and reliable.

Thing I also find interesting is APL interpretation process was much more complicated then one of numpy or matlab, and in particular supported lazy evaluation (thus, construction of expressions was available). This also means that different optimization were possible, like usage of multiple different vector operations at once (like a + b*c), if processor supports this. In russian this is called векторное зацепление, but I haven't found appropriate translation in wiki.

Decision train classifier

During the weekend I've been working on implementation of decision train.

Decision train is a classification/regression model, which I first introduced as a way to speed up gradient boosting with bit-hacks, so it's main purpose was to demonstrate how bit operations coupled with proper preprocessing of data plus revisiting can speed up training/application of trained formula. It may sound strange,  but this incredibly fast algorithm was written in python (all the computationally expensive things are using numpy).

In general, I found it's basic version to be no worse than GBRT from scikit-learn. And my implementation is of course much faster.

Am interesting moment of my approach: I am writing hep_ml library in modular way, for instance, considering GBRT implementation, there are separate classes/functions for:

• boosting algorithm
• base estimator (at this moment only different trees supported, but I believe that one can use any clustering algorithm, like k-means, which sounds a bit strange, but after experiments with pruning this becomes obvious, that getting estimators from some predefined pool can yield very good results in classification)
• loss function. The only difference between regression, binary classification and ranking for gradient boosting is in the loss you use. 
• splitting criterion (actually, this may be called loss function too), the FOM we minimize when building a new tree, but GBRT usually uses simply MSE as such criterion)

This 'modular structure' made it possible to write loss functions once and use with different classifiers.

Maybe, I'll find a time to provide a support of FlatnessLoss inside my neural networks (since it is flexible and uses only gradient, this should be not very complicated).

Numpy tricks (again)

A small post about different tricks I use in numpy. Since frequently I work with statistical data and different statistical criterions, frequently I have to compute the order statistics (alternatively, one can ask about which position will get every number in array after array was sorted)
There is rankdata function inside scipy which solves the named problem:

result = rankdata(array)

However, one can easily implement this is numpy (isn't beautiful?):

result = numpy.argsort(numpy.argsort(array))

which is (almost) equivalent code.
See also my previous post on numpy tricks
PS. Once upon a time I will understand why all blogging platforms suggest heavy and ugly layout templates. Even adding background image turns into non-trivial task.

Benchmarks of speed (Numpy vs all)

Personally I am a big fan of numpy package, since it makes the code clean and still quite fast.
However I am much worried about the speed, so decided to collect different benchmarks

I decided to write my own benchmarks with several non-trivial operation to check things claimed at different sources. Please find it here here. In my tests numpy is not much slower than other libraries or C++.

numpy vs cython vs weave (numpy is about 2 times slower than others)
(posted in 2011)
http://technicaldiscovery.blogspot.ru/2011/06/speeding-up-python-numpy-cython-and.html

Primarily the post is about numba, the pairwise distances are computed with cython, numpy, numba.
Numba is claimed to be the fastest, around 10 times faster than numpy.
(posted in 2013)
https://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2/

Julia is claimed by its developers to be very fast language.
Well, if that is true, there would be no need in writing easily-vectorizable operation in pure python (yep, I mean they are simply cheating), currently these 'benchmarks' are at the main page
http://julialang.org/

The following post was written in 2011, the problem observed is solution of Laplace equation as usual.
Numpy is ~10 times slower than pure C++ solution (and equal to matlab), weave shows the speed comparable with pure C++ (I don't know anyone using weave now)
http://wiki.scipy.org/PerformancePython

Fresh (2014) benchmark of different python tools, simple vectorized expression A*B-4.1*A > 2.5*B is evaluated with numpy, cython, numba, numexpr, and parakeet (and two latest are the fastest - about 10 times less time than numpy, achieved by using multithreading with two cores)
http://nbviewer.ipython.org/github/rasbt/One-Python-benchmark-per-day/blob/master/ipython_nbs/day7_2_jit_numpy.ipynb

Haven't found any general-purpose theano vs numpy benchmarks, but in the article there is comparison of neural networks and theano is expected to give much better speed than numpy/torch(c++)/matlab, specially it is fast on GPU
http://conference.scipy.org/proceedings/scipy2010/pdfs/bergstra.pdf


One more detailed review of numpy vs cython vs c (held in 2014)
http://notes-on-cython.readthedocs.org/en/latest/std_dev.html
Let me copy-paste the example results (computation of std).

Method	Time (ms)	Compared to Python	Compared to Numpy
Pure Python	183	x1	x0.03
Numpy	5.97	x31	x1
Naive Cython	7.76	x24	x0.8
Optimised Cython	2.18	x84	x2.7
Cython calling C	2.22	x82	x2.7

Simple, but comprehensive comparison of python accelerators was prepared by Jake Vanderplaas and Olivier Grisel:
http://nbviewer.ipython.org/github/ogrisel/notebooks/blob/master/Numba%20Parakeet%20Cython.ipynb
Problem is computation of pairwise distances. The results this time seem to be unbiased:

Python	9.51s
Naive numpy	64.7 ms
Numba	6.72ms
Cython	6.57ms
Parakeet	12.3 ms
Cython	6.57 ms

By the way, I was recently looking for a slow place in my code, and it totally dropped of my mind that computation of residual is long operation:
http://embeddedgurus.com/stack-overflow/2011/02/efficient-c-tip-13-use-the-modulus-operator-with-caution/
Integer division / remainder computation time significaly depends on the bit depth (16 bit operation are ~2 times faster than 32 bit ones, division takes roughly 10 times more time than multiplication).

Surprising: multiplication/addition/binary operations take the same time (compared on numpy 1.9)

Optimization of vector operations

Recently worked over optimization of some (secret) classifier. The problem was mostly not in training, but in applying of trained classifier — this code was originally written in C++ and then translated to cython (which decreased the speed by factor of 2).

This was quite easy rewrite the code using numpy and vectorized approach (initally predictions were built event-by event, after rewriting the classifier was applied tree-by-tree). However this gave only speed comparable with original C++ code (so, twixe faster cython).

What really fastened the code is switching from int8 operations to int64 (the latest are natively supported in all modern processors). So 8 operations in int8 were grouped into one 64-bit. This is done pretty simple by creating views:

In[1]: import numpy
In[2]: x = numpy.random.randint(0, 100, size=64000).astype(’int8′)
In[3]: y = numpy.random.randint(0, 100, size=64000).astype(’int8′)
In[4]: %timeit x & y 
10000 loops, best of 3: 60.6 µs per loop
In[5]: %timeit x.view(’int64′) & y.view(’int64′) 
100000 loops, best of 3: 12 µs per loop
# Checked that output is correct
In[6]: numpy.all( (x & y) == (x.view(’int64′) & y.view(’int64’)).view(’int8′) )
Out[6]: True

 

In this simple example we see 5x speed up. Views of course do not copy the data, which is very essential for the speed. This trick can be applied with: summation / subtration / binary or / binary and, but you need that the size of original array was divisible by 8.

By the way, there is awesome collection of twiddling bits, which was my starting point in bit optimizations.