We study the serial and parallel computation of Γ (Gamma), a comprehensive data summarization matrix for linear Gaussian models, widely used in big data analytics. Computing Gamma can be reduced to a single matrix multiplication with the data set, where such multiplication can be evaluated as a sum of vector outer products, which enables incremental and parallel computation, essential features for scalable computation. By exploiting Gamma, iterative algorithms are changed to work in two phases: (1) Incremental-parallel data set summarization (i.e. in one scan and distributive); (2) Iteration in main memory exploiting the summarization matrix in intermediate matrix computations (i.e. reducing number of scans). Most intermediate computations on large matrices collapse to computations based on Gamma, a much smaller matrix. We present specialized database algorithms for dense and sparse matrices, respectively. Assuming a distributed memory model (i.e. shared-nothing) and a larger number of points than processing nodes, we show Gamma parallel computation has close to linear speedup. We explain how to compute Gamma with existing database systems processing mechanisms and their impact on time complexity.
Machine learning models generally require demanding iterative computations on
matrices. This problem becomes significantly more difficult when input matrices
(i.e. the data set) are large, residing on secondary storage, and when they need
to be processed in parallel. We propose the Γ (Gamma) data summarizaton
matrix [
We introduce a general augmented matrix Z, by prefixing it with a row of 1s
and appending a [d + 2]th row to X, containing vector Y (in transposed form).
Since X is d × n, Z is (d + 2) × n, where row [
Multidimensional sufficient statistics [
i=1 xi, Q = XXT . Intuitively, n counts points, L is a simple “linear” sum (a vector) and Q is a “quadratic” matrix. We now introduce Γ , our main contribution, which generalizes previous sufficient statistics. Γ = n LT 1T · Y T
L Q XY T Y · 1 Y XT Y Y T
n P xT P yi = P xi P xiixiT P xiyi
P yi P yixiT P yi2
The fundamental property of Γ is that it can be computed by a single matrix
multiplication with Z: Γ = ZZT . Notice ZZT 6= ZTZ, a property of Gramian
matrices [
Proposition: Γ can be equivalently computed as follows: Γ = ZZT = n X zi · ziT . i=1 (1) 4
Recall from Section 2 that Θ represents a model. We start by presenting the standard iterative algorithm to compute a model Θ in two phases. Phase 1: initialize Θ; Phase 2: Iterate until covergence (i.e. until Θ does not change) reading X from secondary storage at every iteration (once or several times). Our proposed equivalent optimized iterative 2-phase algorithm based on Γ follows. Phase 1: Compute Γ , updating Γ in main memory reading X from secondary storage. Initialize Θ from Γ ; Phase 2: Iterate until convergence exploiting Γ in main memory in intermediate matrix computations.
Lemma (Succinctness–asymptotically small size): Given a large data set X where d n, matrix Γ size is O(d2) O(dn) as n → ∞.
Our time and space analysis is based on the 2-phase algorithm introduced above. Phase 1, which corresponds to the summarization matrix operator, is the most important. Computing Γ with a dense matrix multiplication is O(d2n). On the other hand, computing Γ with a sparse matrix multiplication is O(k2n) for the average case assuming k entries from xi are non-zero; assuming X is hyper-sparse k2 = O(d) then the time is O(dn) on average. In Phase 2 we take advantage of Γ to accelerate computations involving X. Since we are computing matrix factorizations derived from Γ time is Ω(d3), which for a dense matrix it may approach O(d4), when the number of iterations in the factorization numerical method is proportional to d. In short, time complexity for Phase 2 for the models we consider does not depend on n. I/O cost is just reading X with the corresponding number of I/Os in time O(dn) for dense matrix storage and O(kn) for sparse matrix storage. 5
Let N be the number of processing nodes, under a distributed memory
architecture (shared-nothing [
Recall we assume d < n and n → ∞. For Phase 1 each xi ∈ X is hashed to N processors. Then our matrix operator can compute each outer product zi · ziT in parallel, resulting in O(d2n/N ) work per processor for a dense matrix and O(dn/N ) work per processor for a sparse matrix. That is, both dense and sparse matrix multiplication algorithms become optimal as N increases or n → ∞. Given its importance from a time complexity point of view, we present the following result as a theorem:
Theorem (linear speedup): Let Tj be processing time using j nodes, where 1 ≤ j ≤ N . Under our main assumption and Θ fits in main memory then our optimized algorithm gets close to optimal speedup T1/TN ≈ O(N ).
For completeness, we provide an analysis of a parallel system where N is large. That is, N = O(n).
Lemma (sequential communication bottleneck): Assume only one node can receive partial summarizations. Then time complexity with N messages to a single coordinator is O(d2n/N + N d2).
Lemma (improved communication: global summarization with a balanced binary tree): Assume partial results Γ [I] can be sent to other workers in a binary tree fashion. Then T (d, n, N ) ≥ O(d2N +log2(N )d2) (a lower communication bound). 6
Here we connect our matrix-based parallel computation with a database system with: (1) relational queries; (2) extending it with matrix functions. We assume X can be evenly partitioned by i, a reasonable assumption when d n.
We now treat X as a table in order to process with relational algebra, the formalization of SQL. In order to simplify presentation we use X as the table name and we treat Y as another dimension. There are two major alternatives: (1) Defining a table with d columns, plus primary key i X(i, X1, X2, . . . , Xd), ideal for dense matrices. (2) Defining a table with triples, using the matrix subscripts as primary key: X(i, j, v), where i means point id, j dimension subscript and v is the matrix entry value, excluding zeroesm, where v = 0 (i.e. sparse matrix). Alternative (1) requires expressing matrix product as d2 aggregations, which requires a wrapper program to generate the query. Notice Xj is a column name: we cannot assume there is subscript-based access. Alternative (2) can express the computation as a single query, making the solution expressible within relational algebra and it is natural for sparse matrices. In the following query X1 and X2 are “alias” of X because it is referenced twice: Γ = πX1.j,X2.j,sum(X1.v∗X2.v)(X1 1X1.i=X2.i X2). The time complexity to evaluate the Gamma query depends on the join algorithm. Therefore, it can be: O(d2n) (hash join, balanced buckets), O(d2nlog(n)) (sort-merge join, merge join), to O(d2n2) (nested loop join). These bounds can be tighter if the join operation can be avoided using a subscript based mechanism in main memory when Γ fits in main memory.
The previous query solution requires an expensive join (1), which can be
slow. Therefore, our goal is to eliminate the join operation by enabling
efficient subscript-based access. Our proposal is to introduce a matrix function
in an extended database model, where tuples are transfomed into vectors and
a matrix can be the result of the query. This function bridges relational
tables on one side with matrices on the other side. Let Gamma() be an
aggregation which returns a (d + 2) × (d + 2) matrix. There are two elegant solutions:
(1) Γ = Gamma(X1, X2, . . . , Xd)(X) for X with a “horizontal” schema. (2)
Γ = Gamma(X) for a “vertical” schema, where the function assume X has
the schema defined above. In this case time complexity ranges from O(d2n) to
O(d2n log(n)), depending on the aggregation algorithm. The price we are paying
is that the query is no longer expressible within relational algebra and therefore
it is not feasible to treat Γ as a new relational operator. This second mechanism
corresponds to so-called SQL UDFs [