 # Polynomial expansion of the binary classification function

This paper describes a novel method to approximate the polynomial coefficients of regression functions, with particular interest on multi-dimensional classification. The derivation is simple, and offers a fast, robust classification technique that is resistant to over-fitting.

## Code Repositories

### multinomial

A library to store covariants of a multivariate polynomial, with fast evaluation and random access

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Procedure to calculate the classification polynomial

The goal of binary classification is to find the distinction between a signal and background probability distribution. The optimal separation contours are described by Neyman and Pearson (1933), and it is well known these contours can be found by binomial regression (see Bishop, 2006)

. In neural networks it is typically a regression between target values

, and the optimal response function is related to the purity of signal:

 F(x)=s(x)−b(x)s(x)+b(x)=2P(s|x)−1.

By reordering, performing a Fourier transformation and using the Taylor expansion of F(x), it becomes

 ∞∑k=01k!Fk∫Rxk(s(x)+b(x))eiωxdx=∫R(s(x)−b(x))eiωxdx. (1)

The Taylor series of characteristic functions can be expressed with the moments of the corresponding distribution, which is used in the following definition of the and functions and their Fourier transforms and :

 ^g(ω)=∞∑k=0i−k(⟨xk⟩s+⟨xk⟩b)^gkωkk!
 ∂j∂ωjg(ω)=∞∑k=0i−kωkk!⋅^gk+j

 ^h(ω)=∞∑k=0i−k(⟨xk⟩s−⟨xk⟩b)^hkωkk!.

Substituting and back into eq. (1), and exploiting that the Fourier transform of can be expressed with the th derivative of , one gets an equation that is true for every , hence the equation should hold for the coefficients of for every as follows:

 ∞∑k=0i−kωkk!∞∑j=01j!Fj^gk+j=∞∑k=0i−kωkk!^hk
 ^hk=∞∑j=01j!Fj^gk+j. (2)

These equations can be solved for either by a deconvolution, or by finding the solution for the matrix equation

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝^h0^h1⋮^hk⋮⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝^g0^g1⋯^gk⋯^g1^g2⋯^gk+1⋯⋮⋮⋱⋮^gk^gk+1⋯^g2k⋯⋮⋮⋱⋮⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝F0F1⋮Fk⋮⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (3)

where the coefficients were suppressed into the unknowns, which also simplifies the later evaluation of the function.

A possible approximation is using the upper left part of the matrix, and solve the finite system of equations. An example can be seen on fig. 1

, where a 20 degree polynomial was used as a classifier on a Gaussian mixture sample with

events, while the testing was done on an independent events. The resulting separation power is very similar to the theoretical optimum. Figure (c)c clearly shows, that the purity , evaluated in bins of has a monotonic dependence on itself. Lower order approximations of may produce a non-linear, but still monotonically ascending curves, which feature is a requirement for a good classifier, as one can safely say that the events right to a certain value are more signal like than the events to the left.

## 2 Optimisations for multi-dimensional input

In case the dimension of the input is greater than one, the th moment of the distributions become

th order tensors, and similarly

, and . The equation that connects these three tensor series is similar to eq. (2), except that is has free tensor indices:

 ^hkμ1...μk=∞∑j=0^gk+jμ1...μkν1...νjFjν1...νj. (4)

Although this is a tensor equation, the indices of and can be serialised, while can be turned into a rectangular matrix with the serialised indices of as columns, and as rows. These system of equations can be rewritten as a block matrix equation, similar to eq. (3). However, a -dimensional, th order symmetric tensor has only free parameters from the possible . The difference can be many orders of magnitude even for small and , therefore to speed up computations and use less memory, it is beneficiary to compactify the tensors in question, in a way described by Ballard et al. (2011).

For symmetric tensors with degree

, the component belonging to an index vector

is the same for any permutation of . These set of indices can be uniquely identified with the monomial of the index , which is a -dimensional vector, where is the multiplicity of the value in the index vector . The multiplicity of a given monomial is the multinomial coefficient:

 Multiplicity of {m1,...,md}=(nm1,...,md)=n!m1!⋯md!.

The tensor multiplications in question, between the tensors and , can be simplified by running over only the free parameters, indexed with the monomials:

 gk+jμ1...μkν1...νjFjν1...νj=∑m∈m(ν1,⋯,νj)^gμ1...μkmFjm(nm1,...,md).

The multiplicity factor can be factored into , just as the terms before, with the benefit that the same terms should be used when is collapsed with the tensor of an input vector , in order to evaluate . The remaining indices of still hold a -fold symmetry, just as , which can be simplified with the same procedure. In the eq. (4), there is no summation over the free indices of , hence there is no need for the multiplicity factors either. This makes it possible to create a large vector of the serialised tensors, and a symmetric block matrix , containing the rectangular matrix version of for every and .

The difficulty in creating this block matrix is, that despite most of the tensors are used multiple times, they have to be partitioned into matrices in different ways. The efficient way of storing in memory was described above, but to create the matrix versions one needs to access the tensor elements according to the simplified, monomial indices and of order and for the tensor that matches with the structure of indices of and . In case the elements in any of the tensors above are stored serially in lexical index order, then for any index it is true that . For the indices on the diagonal, where , it is possible to calculate the number of elements with higher lexically ordered index, because those indices map the free elements of a dimensional symmetric tensor of order . The same way, the position of the generic index can be found by first finding the position111In this paper it is assumed that the indexing starts from , and the first position is for the diagonal index , then calculating the position for the index, where , but for every . It is done by simply subtracting from the position the number of free elements in a dimensional symmetric tensor of order . Repeating this until the last element of the index, the formula to calculate its position reveals as

 pos(μ)=(k+d−1k)No. free elements of a sym.% tensor of order k−k∑i=1(suborder(k−i+1)+subdimension(d−μi)−1k−i+1).

To match the elements of the serialised tensor with the partitioned matrix, one only has to combine the lexically ordered indices of into one combined index with elements, which can be lexically ordered again with the help of its monomial.

## 3 Tests and conclusions

The following example was made on a three dimensional sample, consisting of twelve non-overlapping, non-symmetric Gaussian peaks as signal and a flat background. The classifier is a 20 degree multi-polynomial, which was found by solving a matrix equation with elements in a few seconds with the solver of the Lapack library (see Anderson et al., 1999). Calculating the elements of this matrix from thousand events takes about 20 seconds on a single core computer.

Figure (a)a shows the histograms of the response, and although the values slightly overshoot , the response vs. purity on fig. (b)b is still a strictly monotonically ascending curve, assuring that

approximates well the optimal classification contours. For higher dimensional inputs, it is usually enough to approximate the classifier function with a low degree polynomial to have a good estimate on the classification contours, or on the separating power of a new variable. Nevertheless, the method seems to be stable against overfitting, since as it is fed with well determined moments of the distributions, and not with the individual events itself; it is not expected to be sensitive to the high frequency noise associated with sampling. The method is also capable of fitting a non-binary target. In this case the

tensors are expectation values of target multiplied with the moments of the input parameters, , while .

As a remark, it must be noted, that since certain distributions have diverging moments, it is beneficiary to transform the distributions into compact phase spaces prior to training, in order to have evaluable results.

## Acknowlegments

I would like to acknowledge the support of my colleges, particularly A. Elizabeth Nuncio Quiroz, Ian C. Brock and Eckhard von Törne.

## References

• Anderson et al. (1999) E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999. ISBN 0-89871-447-8 (paperback).
• Ballard et al. (2011) G. Ballard, T. Kolda, and T. Plantenga.

Efficiently computing tensor eigenvalues on a gpu.

In Parallel and Distributed Processing Workshops and Phd Forum (IPDPSW), 2011 IEEE International Symposium on, pages 1340 –1348, may 2011.
• Bishop (2006) Christopher M. Bishop. Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006. ISBN 0387310738. ISBN 0-387-31073-8.
• Brun and Rademakers (1996) Rene Brun and Fons Rademakers. Root - an object oriented data analysis framework. In Proceedings AIHENP’96 Workshop, Lausanne, volume 389 (1997), pages 81–86. Nucl. Inst. & Meth. in Phys. Res. A, Sep. 1996. See also http://root.cern.ch/.
• Neyman and Pearson (1933) J. Neyman and E. S. Pearson. On the problem of the most efficient tests of statistical hypotheses. Philosophical Transactions of the Royal Society of London. Series A, 231(694-706):289–337, 1933.