Machine learning in COMET experiment (part II)

In previous post I've written a short explanation of COMET - Japanese experiment in particle physics. This second post will be devoted to machine learning approach we developed for tracking (currently this is uses the information which is probably unavailable in online, but I still hope we will be able to get some online information and thus apply part of the approach as high-level trigger).

Local wire features

Anyway, the data we are going to use in classification is:

• energy deposit for each wire
• time of deposit for each wire

Let's look at their distributions:

Energy deposited on each 'wire' by signal and background tracks

Time elapsed starting from the moment of triggering till the time when particles are detected.

How do we obtain these values?
COMET uses straw chambers (which we call wires for simplicity), these are long tubes, which fill detector. They are filled with gas, which is ionized by moving charged particles (mostly, electrons).

After ionization, electrons and ions are moving to opposite directions, so we probably can estimate the moment of ionization by drift time (which is of course much greater compared to time of particle flight in detector). Drift time also gives approximate information about the distance between center of wire and track.

Straw tubes used in many detector systems

Finally, each wire has one more characteristic, namely, the distance to target:

Signal hits are more frequent on the inner layers by construction of experiment.

The time is counted from the moment of when trigger worked (it's another subdetector system). The event is 'recorded' starting from that moment.

The basic algorithm that was proposed (baseline algorithm) is using the cut on energy deposition. As you can see, there is really significant difference: energy deposited by background hits is higher. The reason is that background particles are moving slower, so they ionize more particles.

ROC AUC when we use the only feature - energy, is around 0.95, which seems to be very high. Nevertheless, it's not enough, since we have around 4400 wires, of which around 1000 gets activated (this number is called occupancy) within each measurement (event), while the signal track usually contains around 80 points.

So the event is represented using 4400 pairs (energy, time), of which most are zeros.

So the noise which passes through such basic filter is still very significant and there is large room for improvements.

First, let's combine all the information we have about the single wire (distance from center, time and energy deposition), let's call them wire features:

ROC curves (using physical notation here), we see that usage of wire features made us able to twice decrease the background efficiency. 

Neighbors features

Let me remind you, how the whole picture of hits in COMET detector looks like (we use here projection on the plane, orthogonal to beam line):

Blue dots designate background hits, most of which are tracks of different charged particles, say protons; but some of them are simply noise.
Signal hits, which we are looking for,  are red points forming an arc of circle. It's approximate radius is known ahead.

Simple but useful observation one can see in this plot: signal hits nearly always come together.
This implies that we can try using features collected from neighbours of point.

This drives to dramatic improvement of  classification: we almost get rid of random noise tracks, still there are some coupled misidentified tracks. AUC is about 0.995, fpr (background efficiency) decreased by 4-5 times. In principle it suffices to use only information from left and right wires from the same layer.

Ok, I believe that was simple. Is there still something we missed and that could be improved?
For sure, this is the whole shape of the track. It's time to try using this information.

Hough transform and circle shape

Hough transform was initially developed to detect lines and circles. Being quite trivial, it is one of the most effective algorithms in high energy physics.

After using GBDT trained on wire features and features of its neighbours, we are getting quite clean picture with very few false positives. All we want is to detect and remove possible isolated 'islands' with misidentified background hits.

This is done by very approximate reconstruction of track centers. Since we know approximate radius of track center, we can use the Hough transform with fixed radius. It looks like:

Visualization of Hough transform for circles with fixed radius. We are trying to reconstruct the center of track, going through red points. Assuming that we know the radius of fitted track, all possible centers are laying on the circle with center in red point.

We discretize the space of possible track centers and for each point we reconstruct how likely it is the center of track. It is done using sparse matrices and some normalization + regularization (because otherwise tracks with few or many points will have very low/high probabilities).

Once we computed Hough transform, we leave only those centers, where hough transform is high, applying some nonlinear transformation there and applying inverse hough + some filtering. This way we obtain for each wire the probability that it belongs to some signal track.

Finally, we collect all the information we get for each wire:

• local features (energy deposition, timing, layer_id)
• features collected from neighbors
• result of inverse hough transform

And train GBDT on these features to obtain final classifier. It's ROC AUC is 0.9993 (100 times less probability of misordering)

Final classifier is red and it is very close to ideal one. The ROC AUC is about 0.9993

When we are comparing ROC curves at the threshold of interest (with very high signal sensitivity), things are bit worse, but still very impressing: 

At stable benchmark, the background yield decreased by factor of 34. Original ROC curve is not seen on the plot, since it is much lower.

Visualization of all steps

Initial picture of hits in CyDET. Red are signal tracks, blue are background ones.

After we apply initial GBDT (which uses wire features + neighbors), we have some bck hits (to the right, for instance). See the 'island' of misidentified background to the right.
By greed dots we denote the possible centers of tracks. The bigger the point, the greater value of Hough transform image.

Now we apply some non-linear transformations to leave only centers with very high probability. Then applying inverse hough transform and apply second GBDT, which incorporates also information from inverse Hough transform.

Conclusion

I've described how simple machine learning techniques coupled with well-known algorithms can produce a very good result, superior to many complex approaches.
These results are to be checked on better simulation data, when it will be ready. Also we have some ideas about possibilities of classification for online triggers, which use same ideas.

Links

1. Repository with algorithms and results of tracking.
   Most of plots are taken from it, thanks to Ewen Gillies.
2. Reproducible experiment platform used in experiments, gradient boosting was used from scikit-learn
3. Review of straw chambers. Straw chambers are used within many different experiments due to their high resolution and cheapness.

Multimixture fitting

I was wondering how one can modify Expectation-Maximization procedure for fitting mixtures (well, gaussian mixtures, because it's the only distribution that can be fitted easily) to support really many overlapping summands in mixture.

Randomization probably can be a solution to this problem.

Let me first remind how EM works. There are two steps that are computed iteratively

1. (Expectation) where we compute probability that each particular event belongs to each distribution
2. (Maximization) where given the probabilities we maximize parameters of each distribution.

What if we sample events according to distribution from expectation step? At each stage we will attribute each event to one (in simplest case) component of mixture, or maybe several of them. This kind of randomization should prevent us from 'shrinking' of distribution.

Ok, this again needs time for experiments.

LibFM in python

LibFM is library for factorization machines using an approach proposed by Steffen Rendle. However it seems that there is no python wrapper for this famous library.

However, if you're looking for python version of it, take a look at these projects:

Update: I've tested some LibFM implementations on several datasets. I included libFM and options proposed in this post.

• https://github.com/coreylynch/pyFM
  Lovely minimalistic implementation of Factorization Machines using cython (previous version used numpy).
  This library has an interface similar to scikit-learn.
• http://ibayer.github.io/fastFM/index.html
  fastFM is another option. Library supports both classification and regression.
  Contains three different solvers (SGD, ALS, MCMC). ALS = alternative least squares.
  In the paper written by author of algorithm, he argues that other algorithms are comparable to SGD.
  
  This library is completely following scikit-learn interface (even deriving from BaseEstimator and appropriate mixin classes).
• I was also looking for code in theano, but the only code I found was very dirty and minimalistic (so I'm not hoping it is usable)
  https://github.com/instagibbs/FactorizationMachine

Learning to rank (software, datasets)

Since for some time I've been working on ranking. I was going to adopt pruning techniques to ranking problem, which could be rather helpful, but the problem is I haven't seen any significant improvement with changing the algorithm.

Ok, anyway, let's see what we have in this area.

Datasets for ranking (LETOR datasets)

1. MSLR-WEB10k and MSLR-WEB30k
   You'll need much patience to download it, since Microsoft's server seeds with the speed of 1 Mbit or even slower.
   The only difference between these two datasets is the number of queries (10000 and 30000 respectively). They contain 136 columns,  mostly filled with different term frequencies and so on.
2. Apart from these datasets, LETOR3.0 and LETOR4.0 are available, which were published in 2008 and 2009. Those datasets are smaller. From LETOR4.0 MQ-2007 and MQ-2008 are interesting (46 features there). MQ stays for million queries.
3. Yahoo! LETOR dataset, from challenge organized in 2010. There are currently two versions: 1.0(400Mb) and 2.0 (600Mb). Here is more info about two sets within this data
   

   	Set 1			Set 2
   	Train	Val	Test	Train	Val	Test
   # queries	19,944	2,994	6,983	1,266	1,266	3,798
   # urls	473,134	71,083	165,660	34,815	34,881	103,174
   # features	519			596
4. There is also Yandex imat'2009 (Интернет-Математика 2009) dataset, which is rather small. (~100000 query-pairs in test and in the same train, 245 features). 

And seems that these are most datasets.

Algorithms

There are plenty of algorithms and their modifications created specially for LETOR (with papers).

Implementations

There are many algorithms developed, but checking most of them is real problem, because there is no available implementation one can try. But constantly new algorithms appear and their developers claim that new algorithm provides best results on all (or almost all) datasets. 

This of course hardly believable, specially provided that most researchers don't publish code of their algorithm. In theory,  one shall publish not only the code of algorithms, but the whole code of experiment.

However, there are some algorithms that are available (apart from regression, of course).

1. LEMUR.Ranklib project incorporates many algorithms in C++
   http://sourceforge.net/projects/lemur/
   the best option unless you need implementation of something specific. Currently contains
   
   MART (=GBRT), RankNet, RankBoost, AdaRank, Coordinate Ascent, LambdaMART and ListNet
2. LEROT: written in python online learning to rank framework.
   Also there is less detailed, butlonger list of datasets: https://bitbucket.org/ilps/lerot#rst-header-data
3. IPython demo on learning to rank
4. Implementation of LambdaRank (in python specially for kaggle ranking competition)
5. xapian-letor is part of xapian project, this library was developed at GSoC 2014. Though I haven't found anythong on ranking in documentation, some implementations can be found in C++ code:
   https://github.com/xapian/xapian/tree/master/xapian-letor
   https://github.com/v-hasu/xapian/tree/master/xapian-letor
6. Metric learning to rank (mlr) for matlab
7. SVM-Rank implementation in C++
8. ListMLE, ListNET rankers (seems these were used in xapian)
9. SVM-MAP implementation in C++

Some comparison (randomly sampled pictures from net):

Taken from http://www2009.org/pdf/T7A-LEARNING%20TO%20RANK%20TUTORIAL.pdf

Comparison from http://www.ke.tu-darmstadt.de/events/PL-12/papers/07-busa-fekete.pdf,
though paper was about comparison of nDCG implementations.

Machine learning used in tracking of COMET (japanese particle physics experiment), part I

Recently I worked together with intern from Imperial College over the tracking system of COMET for two months, and here I'm going to briefly sum our (impressing) results. To begin with, let's explain what this experiment about.

From physics courses we know that there are conservation laws, most famous are conservation of  energy and momenta, but there are other. For instance, conservation of different leptonic family numbers.

For instance, electronic number is number of electrons + number of electronic neutrinos. For some time these numbers (electronic number, muonic number, tauonic number) were considered to be preserved, but in the late 1960's it was proved, that neutrino oscillations change the leptonic number of system (neutrino can become electronic from muonic, for instance). This observation is called lepton flavour violation. However, in the Standard Model overall leptonic number (which is sum of three named family numbers) is conserved.

Since that time physicists are searching for other processes which change the leptonic number and include charged leptons (electrons, muons, tauons), which could be a key for new physics, however no success in this direction and the results are currently mostly upper limits we can prove for some processes. Such possible processes have the name of CLFV (charged lepton flavour violation).

Details of COMET experiment

COMET is small experiment (which is only in plans at this moment) built to detect one specific decay (muon to electron conversion on nucleus), which is suppressed in the standard model:
$$ \mu^{-} + N \to e^{-} + N.$$

The frequency of this process is extremely small in SM: it's probability is about $10^{-52}$, which in particular means that we are sure this will not happen during the experiment.

COMET is sending lots of muons (with fixed energy) on aluminum, and makes the muons to interact with it. The process we are looking for is
$$ \mu^{-} + Al \to e^{-} + Al,$$
while this isn't the only way we can obtain electron from muon
$$ \mu^{-} \to e^{-} + \nu_{\mu} + \overline{\nu}_e.$$

Pay attention,  that second process is 'normal' in the sense it doesn't violate the leptonic conservation laws.

The only way to determine whether we had decay-in-orbit (DIO) or not is to measure the energy of resulting electron (since neutrinos are very elusive particles and we are unable to detect it). Fortunately, the distributions of energies of electron is quite different for these processes and background electrons usually have lower energy.

Distribution of energies for signal process and main background (Decay-in-orbit) in COMET

Some other illustrations of COMET experiment

The scheme of planned experiment COMET. Detector (with aluminum target) will be in the end of this pipeline.
This is what planned during phase II, but for the first time experiment will be simpler

COMET during phase I will be much simpler.

Detector of COMET (CyDET). The detector is filled with sensitive wires (white).
When muon hits the aluminum target (in the center), the electron produced travels over helix trajectories of larger radius in magnetic field and hits wires. The radius of trajectory depends on the transverse momentum of electron.

Data collected from COMET

The main information (apart one from trigger) is the data taken from wires, which are registering when charged particle flies nearby (by adsorbing the particles after ionization) and also it measures the time when this happened.

Another characteristic, which helps in distinguishing electrons, is the radius of trajectory, which is computed (in constant magnetic field) by formula
$$r = \dfrac{p_T}{eB}, $$
so for electron we know the distribution of gyroradius from distribution of energies. Gyroradius is how we actually will reconstruct energy of electron.

Typical signal event in COMET, particle is helixing in magnetic field, leaving energy in drift cambers.

In orthogonal projection of CyDET, the collected information from wires will look like this:

Typical picture of signal event in experiment. Red dots are corresponding to signal hits and form a circle of needed radius.
Also one can see here many background hits, which come from other particles flying through detector.

But the shape of circle with some radius is not the only useful information. During signal hits usually less energy is disposed (which may sound a bit contradictory to what I wrote earlier, but the faster particle flies, the less energy it disposes).

The energy deposited is very strong feature, though leaves significant amount of bck hits, which can spoil the picture, when we have high occupancy.

Apart from this, we can use the 'hit time' (measured from the moment when trigger detected particle).

The time after trigger tends to be smaller, while background hits are distributed uniformly across the time.
This is actually to be checked with better monte-carlo simulations (which are quite poor at the moment).

The main goal of our part is classification of signal vs background hits.

In the next post I'll write about how we applied machine learning to this problem and achieved promising results.

PS. There is second part of this post, devoted to machine learning.

Links

1. Official COMET site.
2. Detailed presentation about CLFV, COMET (and e2mu).
3. Раритеты микромира, an article in russian about rare processes (CLFV is one of examples)  
4. COMET phase-I proposal
5. Repository with algorithms and results of tracking.

Automatic reweighting with gradient boosting

In high energy physics there is frequently used trick: reweighting of distribution. The most frequent reason to reweight data is to introduce some kind of constistency between real data and simulated one. The typical picture of how it looks like you can find below:

The aim is to find such weights for MC, that makes original distribution look like target distribution. I've found a way to use here gradient boosting over regression trees, though it is quite different from usual GBRT in loss, update rule and splitting criterion.

Update rule of reweighting GB
Here is the definition of weights:
$$w = \begin{cases}w, \text{event from target distribution} \\ e^{pred} w, \text{event from original distribution} \end{cases}, $$
so as you see, we are looking for multiplier $e^{pred}$, which will reweight the original distribution.

Pred is raw prediction of gradient boosting.

Splitting criterion of reweighting GB
The most obvious way to select splitting is to maximize binned $\chi^2$ statistics, so we are looking for the tree, which splits the space on the most 'informative' cells, where the difference is significant.

$$ \chi^2 = \sum_{\text{bins}} \dfrac{ (w_{target} - w_{original} )^2}{w_{target} + w_{original} } $$

Here I selected a bit more symmetric version, though this was not necessary.

Computing optimal value in the leaf
since we are going to remove difference in distributions, the optimal value is obviously:
$$ \text{leaf_value} = \log{\dfrac{w_{target}}{ w_{original} }} $$

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).

Back to the future (vector representation of future, an approach to analysis of time series)

Today I generated an interesting idea on analysis of time series. Those are frequently used in pattern recognition, mathematical finance and forecastings (weather prediction, sales prediction, number of site visitors etc.).
Apart from different window-based approaches, there are two basic models that by their nature represent the sequential structure in data:

• hidden Markov chain (HMC)
• recurrent neural network (RNN)

Both of the above are actually Markov models and have some internal state (in first case it's probabilities of each of internal states of model, in the case of RNN the state is equal do dumped states of all elements with delay). The notable difference between these models is first is generative, while the second is discriminative.

Some pros and cons of both models I see:

• HMC training is inexpensive, but requires predicted values to be binned (so it doesn't work with continuous variables, even with variables with many possible values)
• RNN requires quite a long training, but has hard times to find stable mapping (since it should be continuous), so requires some regularizations built in its architecture. Training is usually based on prediction of one or several next points in sequence, which is not always adequate.

Apart of setting several next observation in sequence as a target, there is plenty of various targets, like running mean over $n$ next observations. However, I'd prefer machine learning to find out these targets for me. 

Possible approach I came to:

1. train HMC in inverse time, so HMC predicts past, not future
2. it's predictions (probabilities of hidden states) can be treated as vector representation of future, because it is computed on information from future
3. RNN is trained to predict hidden states of HMC
4. Later, some other model is used based on the output of RNN to predict the value of interest. 

So the trick is that we use HMC predictions as some reliable informative target about the future, like this is done is deep learning, for instance.

Theano-based libraries for machine learning

I wrote several times about mathematical vector engine theano and its benefits.

If you're going to dive deep into neural networks, I recommend learning it and using pure theano. However, there are numerous neural network libraries based on theano, let's list some of them:

• Theanets, http://theanets.readthedocs.org/ (pay attention - it is different from theanet, which I haven't found useful)
  theanets is a good option to start,  quite efficient and simple, provides also posssibilities for recurrent neural networks. However, notice, that rprop implementation uses mini-batches, which makes it unstable.
• Keras, http://keras.io/
  so far seems to be very adequate theano library, contains several minibatch-based optimizers and several loss functions, mostly for regression. Authors compare it to Torch.
• Pylearn2, http://deeplearning.net/software/pylearn2/
  this library was written by LISA-lab, authors of theano, however library though being very advanced, itself is terribly complex. Usually it is easier (at lest for me) to write things from the scratch then writing YAML configuration
• others which I consider to be not as mature and partially forgotten:
  lasagne, blocks, crino, deepANN (last one in deprecated), .
• finally, want to note my nano-library (500 lines of code!), which allows using 5 trainers + 6 losses on feedforward networks. Supports weights, and extremely flexible, because it's main paradygm: write an expression. Yes, use theano and just write activation function, blackbox optimization methods will do everything for you. This allows writing amazingly complex activations since you're not more restricted to 'layers' model and fitting of arbitrary functions:
  
  https://github.com/iamfullofspam/hep_ml/blob/master/hep_ml/nnet.py
  
  One more notable thing: it uses scikit-learn interface, so you can use at as a part of, say,  pipeline. Or run AdaBoost over neural networks (which is very fast,  by the way).

Tensor train

One of probable approaches to build graphical models with categorical variables is tensor decomposition:
http://bayesgroup.ru/wp-content/uploads/2014/05/icml2014_NROV-1.pdf
Notably both tensor decomposition (tensor train format) and method to use it for graphical models were developed at my faculty, though by different people at different chairs.
One more interesting question is interpretation of those hidden variables emerging in the middle.
At this moment I'm thinking over possibility to build this into GB-train, since trained model for each event provides sequence of binary visible variables. In principle, this may be written in quite simple way. Provided, that $x_i$ is boolean variable corresponding to $i$th cut in the tree (in the train, to be more precise).

For instance, one can write partition function as $$ Z = A_1[x_1, y] A_2[x_2, y] \dots A_n[x_n, y] $$ or as $$ Z = A_1[x_1] B_1[y] A_2[x_2] B_2[y] \dots A_n[x_n] B_n[y] $$

In both cases it's quite simple to estimate posterior probability since we have only limited set of options (targets) to be checked. But which one should be preferred and why?

Anomaly detection at CERN experiments

Yesterday I dot a brilliant idea on how to implement automatic anomaly detection at CERN experiments. Today this work is done manually - many students / PhD students are looking at different distributions online. This first is quite inreliable, second - it's quite expensive, since you need many people to work all the time (nobody is paid for this - but you anyway spend money on travels).

So, the basic idea is quite simple: one can bin each variable and look at distributions within each of bins. Knowing, that number of events observed inside each bin is Poisson-distributed, one can detect anomalies.

However, this detects only deviations of single variable. How to compute deviations of many variables?
Inside LHCb experiment, for instance, we have topological trigger, which uses gradient-boosted regression trees to filter out events. Trees are actually splitting data into bins, so one can use this and eatimate for one tree the probability of observing anomaly. Here we can apply Wilks theorem, but only for every particular tree, since bins of different trees are correlated.

Links on deep learning

Didn't know where to put it, so just for the memory will post it to blog.

• https://charlesmartin14.wordpress.com/2015/03/25/why-does-deep-learning-work/
  There is much fuzz nowadays about why deep learning works at all (there is no any deep theory under today), and I love reading these hypothetical explanations (though I'm absolutely sure all of them are wrong. A good explanation of success should give you new ideas about what new things will work).
  
  In this couple of articles author is arguing that the action of RBM can be derived as an action of renormalization group. BTW, this is not a first physical analogy in neural network. Apart from RBMs which use Gibbs-like distribution, there were explanations of Hopfiled neural networks via spin glasses and derivation of update rule from mean-field theory.
  
  
• http://colah.github.io/posts/2014-07-NLP-RNNs-Representations/
  This impressive post was written a year ago and shows recent fresh ideas about deep representations of objects. Namely, I was surprised with how common space of representation for different objects may help in translations.
  
  
• http://ai.stanford.edu/~ang/papers/icml09-ConvolutionalDeepBeliefNetworks.pdf
  Finally, a link to a paper where convolutional RBMs were introduced. Using softmax to combine with poll-layer is a good idea.

PS. Found a link of recommended reading for new LISA-lab students. http://www.datakit.cn/blog/2014/09/22/Reading_lists_for_new_LISA_students.html

Vector representations of categories

After thinking a while over the last point in my list of topics in machine learning I was going to think over...

I finally came to the thought that I'd better use some graphical model, since there is a great majority of possible algorithms without any warranty of their workability.

So, first thing - probably one can use LDA (Latent Dirichlet allocation, wiki) to find good representation for categories. In this case we should use bag-of-words model, and `categories` will stand for words.
Yes, the model of LDA doesn't seem to be appropriate, though I'm sure it will give reasonable representations.

Another graphical model I came with, is much simpler (and probably was developed a long time ago, but I don't know).
It's generative model looks looks this:

1. First, a `topic` of event is generated. That's an n-dimensional vector $t_{event} \in \mathbb{R}^n$ drawn from gaussian distribution.
2. For each category, say, `site_id`, each value of site_id has it's own topic $t_{cat.value}$ of same dimension $n$. The greater dot product $(t_{cat.value}, t_{event})$, the higher probability that this value of category will be chosen. I shall stress here, that I'm not thinking over the arbitrary number of categories, I'm currently interested in the case, where the number and types of categories are fixed.
3. To be more precise, the probability of each value within category to be drawn for the event with $t_{event}$ vector is propostional to $$ p(\text{cat_value} | \text{event}) \sim p(\text{cat.value}) \times e^{(t_{cat.value}, t_{event})} $$ so, to compute final probability one shall use softmax function.

That said, currently the problem I expect to meet is extremely low speed of training when the number of values for some category is very large (say, there are cases when number of categories one shall proceed is greater then 10000)

PS. After writing this I understood that usage of decision trees to generate ids of leaves - new 'categories' + usage of LibFM over these new features was a very interesting and good idea.

Plans for ML experiments

Well, since I defended my PhD thesis, I now have more spare time to write something to my blog.

First, I'd better fix the thoughts about things I want to experiment with in the nearest time

• probabilistic gradient boosting, firstly those are modifications of AdaBoost, which model not predictions but distributions for particular events
  
• separate experiments for Oblivious DT probabilistic pruning.
  
• pruning + boosting over bagging + probabilistic GB + optimal selection of values for leaves for RankBoost with proper loss function. I expect this mixture of algorithms and techniques to provide very efficient and reliable method of ranking. At this moment, pruning works fine for classification purposes.
  
• rotation-invariant Conformal FT-like neural networks. Seems that I resolved main issues with final formula, but there are still some problems with pretraining, since I don't use binary hidden variables. PCA is the strongest candidate on pretraining at the moment.
  
• finally, after Avazu CTR contest I came to strong opinion that good results in this area may be achieved only after finding good vector representations for categories (like word2vec for words). This maybe a bit tricky, since I have only ideas about heuristics that may appear useful, not some optimization approach.

Good list of NN training methods

Dived deeper into the methods of training NNs.
Good yet incomplete list of what people did in this area is given in this article dated 2006.

But there is no my favourite Rprop and it's modifications (IRprop+, IRprop-).

BTW, recently dived deeper into experiments with neural networks, and I decided to improve IRprop by keeping track about not simply moving along each axis, but as well along directions
$w_i + w_j$, $w_i - w_j$, being sure that this should speed up training progress.

Well, it increased the speed for the first time, but very fast it stops decreasing loss function and when it is close to the minimal value, it starts serious oscillations and becomes simply unstable.

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)

Why using HDF5?

Update. Everything below is inessential, since I've found the stackoverflow answer about hdf5.

The only thing  don't agree with is blaze, I've tried it and it is obviously raw and needs much time to become not even stable but at least really useful.


My current workflow is completely based on IPython, and I'm working much with pandas (which I personally consider as a good example of poor library design).

Nevertheless,  I moved recently to HDF, though installing pyTables (which is needed to use hdf with pandas) isn't as straightforward as I expected.

And now I convert all the data to hdf.

• First, this usually results in less (about 2-3 times) space needed to store data (but that depends on the dataset, for some ot them there is no difference between csv and hdf).
• Second, your data now stored in binary format, so all the types are strictly defined - no parsing is needed, no guessing of types
• Thus, read/write operations are orders of times faster
• And no approximations are made for float numbers.

Hope, there are enough arguments to move to hdf.

Neural Networks

Just to describe one of my experiments with neural networks.

Neural networs initially were developed as simulation of real neurons, first training rules (i.e. Hebb's rule) were 'reproducing' the behaviour we observe in nature.

But I don't expect this aproach to be very fruitful today. I prefer thinking of neural network as of just one of ways to define function (which is usually called activation function).

For instance, one-layer perceptron's activation function may be written down as
$$f(x) = \sigma( a^i \, x_i )$$
following the Einstein rule, I omit the summation over $i$. $a_i$ are weights.

Activation function for two-layer perceptron ($a^i_j$ and $b^j$ are weights):
$$f(x) =  \sigma( b^j \, \sigma( a^i_j \, x_i )) $$

If one operates the vector variables, and $Ab$ is matrix-by-vector dot product, $\sigma x$ denotes elementwise sigmoid function, then activation function can be written down in pretty simple way:
$$f(x) = \sigma b \sigma A x $$

This is how one can define two-layer perceptron in theano, for instance. Three- or four- layer perceptron isn't more complicated really.

But defining function is only the part of the story - what about training of network? 

I'm sure that the most efficient algorithms won't come from neurobiology, but from pure mathematics. And that is how it is done in today's guides to neural networks: you define activation function, define some figure of merit (logloss for instance), and then use your favourite way of optimization.

I hope that soon the activation functions will be inspired by mathematics, though I didn't succeed much n this direction.

One of activation functions I tried is the following:

First layer:

$$y = \sigma A x $$

Second (pairwise) layer:

$$f(x) = \sigma (b^{ij} y_i y_j ) $$

The difference here that we can use now not only activation of neurons, but introduce some pairwise interaction between them. Unfortunately,  I didn't feel much difference between this modification and simple two-layer network.

Thank to theano,  this is very simple to play with different activation functions :)

So many names

It is often for some simple term to have may different names. But terms used in Machine Learning are absolutely special.

ROC curve! Receiver operating characteristic, sounds like something terribly new.
But this is the same as P-P plot in stats.
(HEP people draw it mirrored frequently)

What about tpr? True positive rate in machine learning. But HEP people use 'signal efficiency', or just 's'.
And survival function in statistic means really the same, though rarely applied to comparison of distributions.

Log-likelihood in stats was renamed to logloss in machine learning. Those who train neural networks, usually call it cross-entropy loss function, but scikit-learn developers call this binomial deviance. And I am still not sure that I know all the names.

For the first time such things just drive you crazy. Later there is no difference for you.
But this turns into real problem much later, because the same things are discussed in different terms, and you do not know which one to use when searching on the internet.