Title: EFFICIENT CALCULATION OF POLYNOMIAL FEATURES ON SPARSE MATRICES

Abstract: We provide an algorithm for polynomial feature expansion that both operates on and produces a compressed sparse row matrix without any densification. For a vector of dimension D, density d, and degree k the algorithm has time complexity
where k is the polynomial-feature order; this is an improvement by a factor d k over the standard method.

1. INTRODUCTION:

Polynomial feature expansion has long been used in statistics to approximate nonlinear functions Gergonne (1974) ; Smith (1918) . The compressed sparse row (CSR) matrix format is a widelyused data structure to hold design matrices for statistics and machine learning applications. However, polynomial expansions are typically not performed directly on sparse CSR matrices, nor on any sparse matrix format for that matter, without intermediate densification steps. This densification not only adds extra overhead, but wastefully computes combinations of features that have a product of zero, which are then discarded during conversion into a sparse format.

We provide an algorithm that allows CSR matrices to be the input of a polynomial feature expansion without any densification. The algorithm leverages the CSR format to only compute products of features that result in nonzero values. This exploits the sparsity of the data to achieve an improved time complexity of O(d k D k ) on each vector of the matrix where k is the degree of the expansion, D is the dimensionality, and d is the density. The standard algorithm has time complexity O(D k ). Since 0 ≤ d ≤ 1, our algorithm is a significant improvement. While the algorithm we describe uses CSR matrices, it could be modified to operate on other sparse formats.

2. PRELIMINARIES:

Matrices are denoted by uppercase bold letters thus: A. The ithe row of A is written a i . All vectors are written in bold, and a, with no subscript, is a vector.

A compressed sparse row (CSR) matrix representation of an r-row matrix A consists of three vectors: c, d, and p and a single number: the number of columns of A. The vectors c and d contain the same number of elements, and hold the column indices and data values, respectively, of all nonzero elements of A. The vector p has r entries. The values in p index both c and d. The ith entry p i of p tells where the data describing nonzero columns of a i are within the other two vectors: c pi:pi+1 contain the column indices of those entries; d pi:pi+1 contain the entries themselves. Since only nonzero elements of each row are held, the overall number of columns of A must also be stored, since it cannot be derived from the other data.

Scalars, vectors, and matrices are often referenced with the superscript k. This is not to be interpreted as an exponent, but to indicate that it is the analogous aspect of that which procedes it, but in its polynomial expansion form. For example, c 2 is the vector that holds columns for nonzero values in A's quadratic feature expansion CSR representation.

For simplicity in the presentation, we work with polynomial expansions of degree 2, but continue to use the exponent k to show how the ideas apply in the general case.

We do provide an algorithm for third degree expansions, and derive the big-O time complexity of the general case.

We have also developed an algorithm for second and third degree interaction features (combinations without repetition), which can be found in the implementation.

3. MOTIVATION:

In this section, we present a strawman algorithm for computing polynomial feature expansions on dense matrices. We then modify the algorithm slightly to operate on a CSR matrix, in order to expose its infeasibility in that context. We then show how the algorithm would be feasible with an added component, which we then derive in the following section.

3.1. DENSE EXPANSION ALGORITHM:

A natural way to calculate polynomial features for a matrix A is to walk down its rows and, for each row, take products of all k-combinations of elements. To determine in which column of A k i products of elements in A i belong, a simple counter can be set to zero for each row of A and incremented efter each polynomial feature is generated. This counter gives the column of A k i into which each expansion feature belongs.

SECOND ORDER (k = 2) DENSE POLYNOMIAL EXPANSION ALGORITHM(A) 1 N = row count of A 2 D = column count of A 3 A k = empty N × D 2 matrix 4 for i = 0 to N -1 5 c p = 0 6 for j 1 = 0 to D -1 7 for j 2 = j 1 to D -1 8 A k icp = A ij1 • A ij2 9 c p = c p + 1

3.2. IMPERFECT CSR EXPANSION ALGORITHM:

Now consider how this algorithm might be modified to accept a CSR matrix. Instead of walking directly down rows of A, we will walk down sections of c and d partitioned by p, and instead of inserting polynomial features into A k , we will insert column numbers into c k and data elements into

d k . INCOMPLETE SECOND ORDER (k = 2) CSR POLYNOMIAL EXPANSION ALGORITHM(A) 1 N = row count of A 2 p k = vector of size N + 1 3 p k 0 = 0 4 nnz k = 0 5 for i = 0 to N -1 6 i start = p i 7 i stop = p i+1 8 c i = c istart:istop 9 nnz k i = |ci| 2 10 nnz k = nnz k + nnz k i 11 p k i+1 = p k i + nnz k i //

Build up the elements of p k , c k , and d k 12 p k = vector of size N + 1 13 c k = vector of size nnz k 14 d k = vector of size nnz k 15 n = 0

16 for i = 0 to N -1 17 i start = p i 18 i stop = p i+1 19 c i = c istart:istop 20 d i = d istart:istop 21 for c 1 = 0 to |c i | -1 22 for c 2 = c 1 to |c i | -1 23 d k n = d c0 • d c1 24 c k n =? 25 n = n + 1

The crux of the problem is at line 24. Given the arbitrary columns involved in a polynomial feature of A i , we need to determine the corresponding column of A k i . We cannot simply reset a counter for each row as we did in the dense algorithm, because only columns corresponding to nonzero values are stored. Any time a column that would have held a zero value is implicitly skipped, the counter would err.

To develop a general algorithm, we require a mapping from columns of A to a column of A k . If there are D columns of A and D k columns of A k , this can be accomplished by a bijective mapping of the following form:

EQUATION

such that 0 ≤ j 0 ≤ j 1 ≤ • • • ≤ j k-1 < D where (j 0 , j 1 , . . . , j k-1 ) are elements of c and p j0j1...i k-1 is an element of c k .

4. CONSTRUCTION OF MAPPING:

Within this section, i, j, and k denote column indices. For the second degree case, we seek a map from matrix indices (i, j) (with 0 ≤ i < j < D ) to numbers f (i, j) with 0 ≤ f (i, j) < D(D-1) 2 , one that follows the pattern indicated by   

EQUATION

where the entry in row i, column j, displays the value f (i, j). We let T 2 (n) = 1 2 n(n + 1) be the nth triangular number; then in Equation 2, column j (for j > 0) contains entries with T 2 (j -1) ≤ dD, is therefore

O dD + k -1 k = O (dD + k -1)! k!(dD -1)! (8) = O (dD + k -1)(dD + k -2) . . . (dD) k! (9) = O ((dD + k -1)(dD + k -2) . . . (dD)) for k dD (10) = O d k D k (11) 6.2 EMPIRICAL

To demonstrate how our algorithm scales with the density of a matrix, we compare it to the traditional polynomial expansion algorithm in the popular machine library scikit-learn Pedregosa et al. (2011) in the task of generating second degree polynomial expansions. Matrices of size 100 × 5000 were randomly generated with densities of 0.2, 0.4, 0.6, 0.8, and 1.0. Thirty matrices of each density were randomly generated, and the mean times (gray) of each algorithm were plotted. The red or blue width around the mean marks the third standard deviation from the mean. The time to densify the input to the standard algorithm was not counted.

The standard algorithm's runtime stays constant no matter the density of the matrix. This is because it does not avoid products that result in zero, but simply multiplies all second order combinations of features. Our algorithm scales quadratically with respect to the density. If the task were third degree expansions rather than second, the plot would show cubic scaling.

The fact that our algorithm is approximately 6.5 times faster than the scikit-learn algorithm on 100× 5000 matrices that are entirely dense is likely a language implementation difference. What matters is that the time of our algorithm increases quadratically with respect to the density in accordance with the big-O analysis.

Figure 1 : Our algorithm (bottom) scales with the density of a matrix, unlike the traditional polynomial feature expansion method (top). The task was a second degree expansion, which is why the time of our algorithm scales quadratically with the density.

:



Figures and Tables:

Peer review (iclr2017):
This paper proposes an algorithm for polynomial feature expansion on CSR matrices, which reduces the time complexity of the standard method by a factor d^k where d is the density of the sparse matrix. The main contribution of this work is not significant enough. The experiments are incomplete and not convincing. The background of the problem is not sufficiently introduced. There are only two references in the introduction part (overall only three papers are cited), which are from decades ago. Many more relevant papers should be cited from the recent literature. The experiment part is very weak. This paper claims that the time complexity of their algorithm is O(d^k D^k), which is an improvement over standard method O(D^k) by a factor d^k. But in the experiments, when d=1, there is still a large gap (~14s vs. ~90s) between the proposed method and the standard one. The authors explain this as *likely a language implementation*, which is not convincing. To fairly compare the two methods, of course you need to implement both in the same programming language and run experiments in the same environment. For higher degree feature expansion, there is no empirical experiments to show the advantage of the proposed method. Some minor problems are listed below. 1) In Section 2, the notation *p_i:p_i+1* is not clearly defined. 2) In Section 3.1, typo: *efter* - *after* 3) All the algorithms in this paper are not titled. The input and output is not clearly listed. 4) In Figure 1, the meaning of the colored area is not described. Is it standard deviation or some quantile of the running time? How many runs of each algorithm are used to generate the ribbons? Many details of the experimental settings are missing.