Showing posts with label Optimization. Show all posts
Showing posts with label Optimization. Show all posts

Tuesday, June 30, 2015

Thursday, June 25, 2015

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.

Monday, June 22, 2015

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.

Monday, May 4, 2015

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.

Tuesday, January 6, 2015

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

MethodTime (ms)Compared to PythonCompared to Numpy
Pure Python183x1x0.03
Numpy5.97x31x1
Naive Cython7.76x24x0.8
Optimised Cython2.18x84x2.7
Cython calling C2.22x82x2.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:

Python9.51s
Naive numpy64.7 ms
Numba6.72ms
Cython6.57ms
Parakeet12.3 ms
Cython6.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)

Thursday, January 1, 2015

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.



Saturday, November 1, 2014

No Free Lunch

Until today, I was absoltely sure that "No Free Lunch Theorem" is rather proverb to say "There is no easy way to solve problem in very general case, better results can be obtained only in particular cases".

I'd never expect there is really a mathematical theorem which proves that (for optimization problems):