pylspack: Parallel algorithms and data structures for sketching, column subset selection, regression and leverage scores

1 Introduction

• Sketching uses random linear maps to reduce problem size and derive fast exact or approximate solutions with high probability.The introduction identifies random sampling and random projections as typical approaches.Source: §1 · p. 1

1.1 Contributions and structure of this paper

• Parallel algorithms and data structures for CountSketch and Gaussian projections support column subset selection, least-squares regression, and leverage-score computation.The work also analyzes tall-and-skinny Gramian computation and implements the methods in pylspack.Source: §1.1 · p. 2
• pylspack is a multithreaded Python library with a C++ core and OpenMP parallelization, compatible with SciPy and NumPy matrix structures.Experiments report that it can significantly outperform existing open-source and vendor libraries, especially for tall-and-skinny matrices.Source: §1.1 · pp. 2–3

1.2 Related work

• Existing parallel-sketching work focuses on a limited set of transforms, while OSNAP and CountSketch have been overlooked.The related-work discussion also describes software limitations in Python support and shared-memory focus.Source: §1.2 · p. 3

2 Sketching algorithms and data structures

• CountSketch-focused methods are presented as adaptable to any OSNAP, since CountSketch is an OSNAP with one nonzero per column.The section first proposes data structures for CountSketch transforms.Source: §2 · p. 4
• Gaussian and CountSketch transforms are combined in parallel to obtain smaller sketching sizes.The section describes this combination as a common RandNLA technique.Source: §2 · p. 4
• A memory-efficient parallel Gaussian sketching algorithm is included for completeness and as part of the software library.The section describes it as similar to existing implementations.Source: §2 · p. 4

2.1 CountSketch

• The CountSketch data structure targets low memory, parallel construction and application, and computation of GSA in O((mdr + nnz(A))/p) operations.The output SA is stored as a dense row-major array.Source: §2.1 · p. 4

2.1.1 Vector representation

• A vector of n signed integers represents CountSketch by recording each column’s nonzero row and sign, but this representation does not lend itself to parallel algorithms.The section says additional partitioning by rows can remedy the parallelism issue.Source: §2.1.1 · p. 4, Figure 1

2.1.2 Partitioning

• Block CountSketch partitions S into a two-dimensional array of submatrices, enabling row-major computation of SA that eliminates write conflicts.The partition dimensions are br × bc blocks, each of size nr × nc.Source: §2.1.2 · p. 5, Figure 2 · p. 5

2.2.1 Construction

• Construction is parallelized across columns, while multiplication SA is parallelized across rows of the partitioned matrix.Columns are independent, though dependencies remain within each column; complexity depends on the data structure used for the blocks.Source: §2.2 · p. 5
• For tall-and-skinny inputs with n much larger than d^2, the product SA is treated and stored as dense, with total work O(nnz(A)).The text notes that each row of A is used exactly once; probabilistic guarantees for parallel workload are analyzed separately.Source: §2.2.1 · pp. 5–6

2.2.3 Computing the product 𝐺𝑆𝐴(CountGauss)

• The CountGauss product GSA is computed in batches to reduce temporary-storage requirements without changing complexity.Parallelization occurs within each individual step of the algorithm.Source: §2.2.3 · p. 6

2.3 Choosing a data structure for the submatrices

• Unordered COO supports constant-time insertion by appending entries and computes SA in O(nnz(A)) by iterating over nonzeros of S.For each nonzero, the corresponding row of A is scaled and used to update the result.Source: §2.3 · p. 6
• The CountGauss algorithm processes CountSketch rows in batches, forms B = S_iA, and updates C with the Gaussian product G B.It allocates auxiliary matrices G and B for these batched computations.Source: §2.3 · p. 7
• Block Column CountSketch stores one signed column index per nonzero and r×b_n pointers to submatrices, with n additional sign bits.Each block has one row; insertion and multiplication complexity match COO.Source: §2.3 · p. 7

2.4 Probabilistic bounds for processor workload

• For a randomly selected subset of CountSketch rows, workload is modeled as a weighted sum of independent Bernoulli variables, with weights given by row sparsities.The expected workload is tied to the fraction of rows assigned to the subset.Source: §2.4 · p. 7
• Hoeffding’s inequality is applied to centered, bounded workload contributions to bound fluctuations for a fixed processor row set.The analysis then extends the bound across processors with a union bound.Source: §2.4 · pp. 8–9
• The workload bounds hold simultaneously for all p processors with probability at least 1 − 2p/n.The simultaneous guarantee follows by a union bound over processor row sets.Source: §2.4 · p. 9
• Block-COO-CountSketch and Block-Column-CountSketch satisfy the stated data-structure conditions deterministically except for condition (2b), which holds in expectation and is well centered with high probability.Source: §2.4 · p. 9

2.5 Gaussian sketching

• For dense A, Gaussian sketching can use BLAS routines; for sparse CSR A, the algorithm exploits sparsity to compute GA.The resulting matrix is C = GA for an m×n Gaussian matrix G.Source: §2.5 · pp. 9–10
• Assigning each process a block of rows of G lets it compute its entire product G_pA without write conflicts or synchronization, with equal operation counts per process.Source: §2.5 · p. 9
• Fast Gaussian random-number generation is crucial because otherwise generating the variables can overwhelm total performance.The text specifically notes this concern for the parallel algorithm.Source: §2.5 · p. 9

3 Gram matrix and Euclidean row norms

• When A is sufficiently tall, its Gram matrix is likely dense, and a simple outer-product algorithm uses Θ(R) scalar multiplications, where R is the required count for standard inner- or outer-product algorithms.General sparse-matrix multiplication methods often target distinct sparse inputs and sparse outputs, unlike this tall-and-skinny Gram case.Source: §3 · pp. 10–11

3.1 Parallel Gram 1: A memory efficient algorithm

• The memory-efficient Gram algorithm needs no intermediate storage or processor communication, but its time is ˜Θ(n + nnz²(A)/p) under balanced row and column nonzero distributions.Sorted CSR column indices permit binary searches to filter columns for each processor’s assigned output rows.Source: §3.1 · p. 11

3.2 Parallel Gram 2: A faster algorithm

• For tall-and-skinny inputs with n ≫ d², row-block partitioning gives total time ˜O(nnz²(A)/p + d²), substantially less than the memory-efficient variant.Each processor forms a local dense Gram matrix; partial results are then summed, requiring p·d² additional storage.Source: §3.2 · pp. 11–12

3.3 Computing Euclidean row norms

• For sparse CSR inputs, squared row norms of AB can be computed in O(nnz²(A)) operations, avoiding explicit formation of the product.The outer loop has no write conflicts and can be parallelized; workers dynamically take another row after finishing one.Source: §3.3 · p. 12

4.1 Column subset selection

• For column subset selection, CountGauss sketches A with m ≈ 2d and r ≈ d², then pivoted QR selects the top-k columns.The stated parameter choice is intended to provide decent speed and approximation guarantees.Source: §4.1 · p. 13

4.2 Least squares regression

• Least-squares options include a Gram-pseudoinverse direct solver, approximate sketch-and-solve, and sketch-based preconditioning followed by an iterative solve.The direct Gram approach forms AᵀA; sketch-and-solve solves the sketched problem, while preconditioning uses a parallel iterative method.Source: §4.2 · pp. 13–14

4.3 Leverage scores

• For leverage scores, the algorithms use a Gram pseudoinverse, column subset selection, or approximate methods suited to dense matrices.The exact approach uses the Gram matrix pseudoinverse as an orthogonalizer; the subset approach selects columns first.Source: §4.3 · p. 14

5 Software: the pylspack package

• The publicly available package provides Python bindings for C++ kernels parallelized with OpenMP, compatible with NumPy and SciPy data structures.Its API includes sketching, Gram-matrix, squared-row-norm, column-selection, and leverage-score functions; least-squares solvers and preconditioning can be derived from existing calls.Source: §5 · pp. 14–16

5.1 Pseudo-random number generation

• Each thread constructs its own random engine with a seed based on thread id, using STL random-number methods without external dependencies.The authors identify specialized parallel generators such as Random123 as future work.Source: §5.1 · p. 16

5.2 Unit Testing

• The test suite checks kernels across input parameters, sizes, sparsity patterns, and thread configurations, including floating-point errors and algorithm-specific guarantees.Automated testing combines two compilers, three Python versions, and single- and multi-threaded execution.Source: §5.2 · pp. 16–17

6 Experiments

• Experiments compare methods on sparse inputs and, where applicable, dense inputs, across two matrix sizes, on a single-socket AMD EPYC machine with 16 physical cores and up to 32 threads.The sizes were 2,097,152 × 512 (sparse density 5%) and 131,072 × 8,192; the latter was tested in dense and sparse forms, with 2.5% density in the sparse form.Source: §6 · p. 17, §6.1 · p. 17

6.1 Test matrices and parameters

• The test setup uses target dimension k=512 for the tall matrix and k=100 for the shorter-and-wider matrix; sketching runs use 15 repetitions in three rounds.The paper notes that the shorter-and-wider matrix does not meet the stated CountSketch subspace-embedding size requirement.Source: §6.1 · p. 17

6.2 Gaussian sketching

• For Gaussian projections, csrjlt is clearly faster than GRP and JLT on both matrices and scales better as thread count increases.On the tall matrix, JLT is slightly faster with one thread; csrjlt slows at four threads on the short matrix, which the authors attribute to processor architecture and memory-bound sparse multiplication.Source: §6.2 · p. 18

6.3 Sparse sketching

• Across sparse and dense inputs and both matrix sizes, pylspack’s sparse-sketching methods are faster and scale better; at 32 threads, they exceed competitors by more than one order of magnitude.The comparison includes sklearn, whose methods are not multithreaded, and Skylark, which targets distributed rather than shared-memory parallelism.Source: §6.3 · p. 18

6.4 Multi-level sketching

• csrcgs/rmcgs are many times faster and scale better than the alternatives for computing GSA, and their batched organization reduces intermediate storage requirements.For the tall matrix, the reported intermediate storage is about 4.2 MB for csrcgs/rmcgs versus roughly 200 MB for SA.Source: §6.4 · pp. 19–20

6.5 Gram matrix computation

• For tall-and-skinny Gram computation, csrrk is over one order of magnitude faster than the other methods at 32 threads; the authors recommend it when the Gramian is expected to be dense.For the short matrix, GM-MKL-sparse is about twice as fast as csrrk, but storing a dense-output Gramian sparsely entails significant redundancy.Source: §6.5 · p. 21

• The algorithms scale well on shared-memory machines and can outperform existing open-source and vendor libraries in performance and memory requirements.Source: §7 · p. 21
• The Gramian runtime comparison covers tall sparse 2M×512 inputs at ∼5% density and rectangular 131K×8K inputs at ∼2.5% density.Source: Figure 8 · p. 22
• Large fluctuations in expected processor workload during CountSketch sketching occur with only small probability, motivating further load-balancing analysis and more aggressive algorithms.Source: §7 · p. 22
• Leverage scores can be computed from the Gramian and its pseudoinverse followed by row norms, but the paper states there is no general reduction in the opposite direction.Source: §7 · p. 22