Principal Component Analysis and Higher Correlations for Distributed Data

Introduction

In modern large-scale machine learning problems the input data is often distributed among many servers, while the communication as well as time and space resources per server are limited. We consider two well-studied problems: (1) Principal Component Analysis (PCA), and (2) Generalized Higher-order correlations. Both problems study correlations between vectors. For the first problem, the vectors correspond to the rows of a matrix and we are interested in second-order correlations, while in the second problem we are interested in higher-order correlations among the vectors.

PCA is a central tool in many learning algorithms. The goal of PCA is to find a low-dimensional subspace that captures as much of the variance of a dataset as possible. By projecting the rows of a matrix onto this lower-dimensional subspace, one preserves important properties of the input matrix, but can now run subsequent algorithms in the lower-dimensional space, resulting in significant computational and storage savings. In a distributed setting, by having each server first locally project his/her own data onto a low-dimensional subspace, this can also result in savings in communication. PCA is useful for a variety of downstream tasks, e.g., for clustering or shape-fitting problems () and latent semantic analysis.

The second problem we consider is the Generalized Higher Order Correlation Problem. For this problem we assume server $t$ has an $n$-dimensional vector $a_{t}$ with non-negative entries. Note that for PCA, it is useful and more general to allow the entries to be positive, negative, or zero. On the other hand, the non-negativity assumption for Generalized Higher Order Correlations is justified both by the applications we give , as well as the fact that it is impossible to achieve low communication without this assumption, as described in more detail below.

A special case of this problem is the well-studied frequency moment problem. That is, if server $t$ holds the vector $a_{t}$, with coordinates $a_{t1},a_{t2},\ldots,a_{tn}$, then the $k$-th frequency moment of $\sum_{t=1}^{s}a_{t}$ is $\sum_{i=1}^{n}(\sum_{t=1}^{s}a_{ti})^{k}$, where, $k$ is a positive integer. This problem has been extensively studied in the data stream literature, starting with the work of . Known lower bounds for this problem from that literature rule out low communication algorithms when $k>2$ in the distributed setting when the number of servers grows as a power of $n$ (), or when there are only two servers and the entries are allowed to be negative . Here we overcome these lower bounds for smaller $s$ and indeed will develop algorithms and lower bounds for estimating $\sum_{i=1}^{n}f(\sum_{t=1}^{s}a_{ti})$, for a general class of functions $f:{\bf R}_{+}\rightarrow{\bf R}_{+}$.

We then extend these results to the following more general problem: there is a collection of vectors that is partitioned into $s$ parts - $W_{1},W_{2},\ldots,W_{s}$ - and server $t$ holds $W_{t}$. For each $t$ and each $i\in W_{t}$, there is an $n$-dimensional vector $v_{i}=(v_{i1},v_{i2},\ldots,v_{in})$ wholly residing on server $t$. Let $f:{\bf R}_{+}\rightarrow{\bf R}_{+}$ and $g:{\bf R}_{+}^{k}\rightarrow{\bf R}_{+}$ be functions. For a natural number $k$, define the $k$-th generalized moment $M(f,g,k)$ as

There are many applications of higher-order correlations, and we only mention several here. For a document collection, we seek statistics (second, third and higher moments) of the number of documents in which each trigram (triples of terms) occurs. For a bipartite graph $G(V_{1},V_{2},E)$ and constants $(r,u)$, we want to estimate the number of $K_{r,u}$ (complete bipartite graph) subgraphs. For a time series of many events, we want to estimate the number of tuples $(E_{1},E_{2},\ldots,E_{r};t_{1},t_{2},\ldots,t_{u})$ for which each of the events $E_{1},E_{2},\ldots,E_{r}$ occurs at each of the times $t_{1},t_{2},\ldots,t_{u}$.

Conceptually, for each $i$, we can think of a vector $a_{i}$ with ${n\choose k}$ components - one for each distinct tuple $(j_{1},j_{2},\ldots,j_{k})$. Suppose $a_{i;j_{1},j_{2},\ldots,j_{k}}=g(v_{i,j_{1}},v_{i,j_{2}},\ldots,v_{i,j_{k}}),$ and let $a_{t}=\sum_{i\in W_{t}}a_{i}.$ Our first theorem describes a way of estimating $M(f,g,k)$ up to a $(1+\varepsilon)$-factor, where, each server uses polynomial time and polynomial space, but we try to optimize total communication while keeping the number of rounds constant. For this algorithm, server $t$ explicitly constructs the vector $a_{t}$ first, so it uses $O(n^{k}|W_{t}|)$ space. Thereafter the space is linear in the total size of all the $a_{t}$. Our second theorem shows how to reduce space to linear in $n$. This algorithm does not construct $a_{t}$ explicitly, but instead performs a rejection sampling procedure.

Before stating our theorems, we need some notation. Let $c_{f,s}$ be the least positive real number such that

Note that for $f(x)=x^{k}$ (as in the $k$-th frequency moment), $c_{f,s}=s^{k-1}$, since for any non-negative real numbers $b_{1},b_{2},\ldots,b_{s}$, we have $(b_{1}+b_{2}+\cdots+b_{s})^{k}\leq s^{k-1}(b_{1}^{k}+b_{2}^{k}+\cdots+b_{s}^{k}),$ and taking $b_{t}=1$, we see that the factor $s^{k-1}$ cannot be improved.

The communication and computations are not assumed to be synchronous. We arbitrarily denote one of the $s$ servers as the Central Processor (CP). A round consists of the CP sending a message to each server and each server sending an arbitrary length message to the CP. A round is complete when the CP has received messages from all servers from which it is expecting a message in that round. All servers communicate only with the CP, which, up to a factor of two, is equivalent to the servers communicating directly with each other (provided they indicate in their message who the message is being sent to). For formal details of this model, we refer the reader to Section 3 of . Our algorithms take polynomial time, linear space and $O(1)$ rounds of communication.

Our Results

Low-rank matrix approximation and approximate PCA.

Our first set of results is for low-rank approximation: given an $n\times d$ matrix $A$, a positive integer $k$ and $\varepsilon>0$, find an $n\times d$ matrix $B$ of rank at most $k$ such that

Here, for a matrix $A$, the Frobenius norm $||A||_{F}^{2}$ is the sum of squares of the entries of $A$. A basis for the rowspace of $B$ provides an approximate $k$-dimensional subspace to project the rows of $A$ onto, and so is a form of approximate PCA. We focus on the frequently occurring case when $A$ is rectangular, that is, $n\gg d$.

Consider the arbitrary partition model where an $n\times d$ matrix $A_{t}$ resides in server $t$ and the data matrix $A=A^{1}+A^{2}+\cdots+A^{s}$. For any $1\geq\varepsilon>0$, there is an algorithm that, on termination, leaves a $n\times d$ matrix $C^{t}$ in server $t$ such that the matrix $C=C^{1}+C^{2}+\cdots+C^{s}$ is of rank $k$ and with arbitrarily large constant probability achieves $\|A-C\|_{F}\leq(1+\varepsilon)\min_{X:\text{rank}(X)\leq k}||A-X||_{F},$ using linear space, polynomial time and with total communication complexity $O(sdk/\varepsilon+sk^{2}/\varepsilon^{4})$ real numbers. Moreover, if the entries of each $A_{t}$ are $b$ bits each, then the total communication is $O(sdk/\varepsilon+sk^{2}/\varepsilon^{4})$ words each consisting of $O(b+\log(nd))$ bits.

In contrast to the guarantees in Theorem 1.1, in the streaming model even with multiple passes, a simple encoding argument formalized in Theorem 4.14 of shows the problem requires $\Omega(n+d)$ communication. We bypass this problem by allowing the $s$ different servers to locally output a matrix $C_{t}$ so that $\sum_{t}C_{t}$ is a $(1+\varepsilon)$-approximation to the best rank-$k$ approximation. We are not aware of any previous algorithms with less than $n$ communication in the arbitrary partition model.

In the row-partition model, in which each row of $A$ is held by a unique server, there is an $O(sdk/\varepsilon)$ word upper bound due to . This is also achievable by the algorithms of . As the row-partition model is a special case of our model in which for each row of $A$, there is a unique server with a non-zero vector on that row, our result implies their result up to the low order $O(sk^{2}/\varepsilon^{4})$ term, but in a stronger model. For example, consider the case in which a customer corresponds to a row of $A$, and a column to his/her purchases of a specific item. These purchases could be distributed across multiple servers corresponding to different vendors. Or in the case of search data, each column could correspond to a search term of a user, and the searches may be distributed across multiple servers for storage and processing considerations. These examples are captured by the arbitrary partition model but not by the row partition model.

The technique for our upper bound is based on a two-stage adaptive sketching process, and has played an important role in several followup works, including CUR Matrix Factorizations of and subspace embeddings for the polynomial kernel by .

Frequency moments and higher-order correlations.

Our next set of results are for estimating higher moments and higher-order correlations of distributed data.

Let $f:{\bf R}_{+}\rightarrow{\bf R}_{+}$ and $c_{f,s}$ be as in (1). There are $s$ polynomial time, linear space bounded servers, where server $t$ holds a non-negative $n$-vector $a_{t}=(a_{t1},a_{t2},\ldots,a_{tn})$. We can estimate $\sum_{i=1}^{n}f\left(\sum_{t=1}^{s}a_{ti}\right)$ up to a $(1+\varepsilon)$ factor by an algorithm using $O(s^{2}c_{f,s}/\varepsilon^{2})$ total words of communication (from all servers) in $O(1)$ rounds. Moreover, any estimation up to a $(1+\varepsilon)$ factor needs in the worst case $\Omega(c_{f,s}/\varepsilon)$ bits of communication.

We remark that the lower bound applies to any function $f$ with parameter $c_{f,s}$, not a specific family of such functions.

Let $f:{\bf R}_{+}\rightarrow{\bf R}_{+}$, $g:{\bf R}_{+}^{k}\rightarrow{\bf R}_{+}$ be monotone functions with $c_{f,s}$ as in (1). $k\in O(1)$ is a natural number and let $M(f,g,k)$ be the generalized moment. We can approximate $M(f,g,k)$ to relative error $\varepsilon$ by an algorithm with communication at most $O(s^{3}c_{f,s}/\varepsilon^{2})$ words in $O(1)$ rounds. Further, we use polynomial time and linear space.

A key feature of this algorithm, and our following ones, is worth noting: they involve no dependence on $n$ or $\ln n$, so they can be used when $a_{t}$ are implicitly specified and $n$ itself is very large, possibly infinite (provided, we can communicate each index $i$). In the theorem below $\Omega$ is the set of coordinates of each vector. It is analogous to $[n]$. We use $\sum_{x\in\Omega}$, which when $\Omega$ is infinite and the probabilities are densities, should be replaced with an integral; our theorem is also valid for the case when we have integrals.

Let $f:{\bf R}_{+}\rightarrow{\bf R}_{+}$, $g:{\bf R}_{+}^{k}\rightarrow{\bf R}_{+}$ be monotone functions with $c_{f,s}$ as in (1). Server $t$ is able to draw (in unit time) a sample $x\in\Omega$ according to a probability distribution $h_{t}$ on $\Omega$. Also, server $t$ can estimate $\sum_{x\in\Omega}f(h_{t}(x))$. Then with $O(s^{3}c_{f,s}/\varepsilon^{2})$ words of communication, CP can estimate $\sum_{x\in\Omega}f\left(\sum_{t=1}^{s}h_{t}(x)\right)$ to within relative error $\varepsilon$.

As a special case we consider the well-studied case of frequency moments. The best previous upper bound for the $k$-th frequency moment problem in the distributed setting is by who gave an algorithm that achieves $s^{k-1}\left(\frac{C\log n}{\varepsilon}\right)^{O(k)}$ communication, so the complexity still depends, albeit mildly, on $n$. Theorem 1.3 implies an algorithm with $O(s^{k+1}/\varepsilon^{2})$ words of communication. We further improve this:

There are $s$ servers, with server $t$ holding a non-negative vector $a_{t}=(a_{t1},a_{t2},\ldots,a_{tn})$.The vector $a_{t}$ need not be written down explicitly in server $t$. it just has to have the ability to (i) find $\sum_{i=1}^{n}a_{ti}^{k}$ to relative error $\varepsilon$ and draw a sample according to $\{a_{ti}^{k}/\sum_{j}a_{tj}^{k}\}$. Then, to estimate $A=\sum_{i=1}^{n}\left(\sum_{t=1}^{s}a_{ti}\right)^{k}$ to within relative error $\varepsilon$, there is an algorithm that communicates $O((s^{k-1}+s^{3})(\ln s/\varepsilon)^{3})$ words Each communicated word is either an index $i$ or a value $a_{ti}$. in $O(1)$ rounds.

Low-rank Approximation

For a matrix $A$, define $f_{k}(A)$ as: $f_{k}(A)=\min_{X:\text{rank}(X)\leq k}||A-X||_{F}.$ Recall that the rank-$k$ approximation problem is the following: Given an $n\times d$ matrix $A$, and $\varepsilon>0$, find an $n\times d$ matrix $B$ of rank at most $k$ such that $||A-B||_{F}\leq(1+\varepsilon)\cdot f_{k}(A)$.

We are mainly concerned with communication, so we omit the tradeoffs of different compositions and just use a composition for which $m=O(d/\varepsilon^{2})$, $PA$ is an $m\times d$ matrix of words each consisting of $O(b+\log(nd))$ bits, and $P$ can be specified using $O(d\log n)$ bits (using a $d$-wise independent hash function, as first shown in ), see Theorem 2.1 below. Since we will assume that $b$ is at least $\log n$, the $O(d\log n)$ bits to specify $P$ will be negligible, though we remark that the number of bits to specify $P$ can be further reduced using results of .

We will prove the following property about the top $k$ right singular vectors of $PA$ for a subspace embedding $P$.

We will combine this property with the following known property.

We can now state the algorithm, which we call AdaptiveCompress.

In AdaptiveCompress, the matrix $P$ is of size $O(k/\varepsilon^{3})\times n$.

(of Theorem 2.1.) Suppose $P$ is a subspace embedding for the column space of $A$. Form an orthonormal basis of ${\bf R}^{d}$ using the right singular vectors of $PA$. Let $v_{1},v_{2},\ldots,v_{d}$ be the basis.

Also, suppose now $u_{1},u_{2},\ldots,u_{d}$ is an orthonormal basis consisting of the singular vectors of $A$. Then, we have

Thus,$\|A-A\sum_{i=1}^{k}v_{i}v_{i}^{T}\|_{F}^{2}\leq(1+\varepsilon)^{4}f_{k}(A)^{2},$ as desired.

For the second part of the theorem, regarding the choice of $P$, fix attention on one particular $x\in{\bf R}^{d}$. We apply Theorem 2.2 of with $A,B$ of that theorem both set to $Ax$ of the current theorem and $m=O(d/\varepsilon^{2})$ in the notation of that theorem. This states that for $m=O(d/\varepsilon^{2})$, if $P$ is an $m\times n$ matrix with $O(d)$-wise independent entries uniform in $\{-1/\sqrt{m},+1/\sqrt{m}\}$, then for any fixed vector $x$, $\|PAx\|_{2}=(1\pm\varepsilon)\|Ax\|_{2}$ with probability $1-\exp(-d)$. We combine this with Lemma 4 in Appendix A of , based on , to conclude that for all vectors $x$, $\|PAx\|_{2}=(1\pm\varepsilon)\|Ax\|_{2}$ with probability $1-\exp(-d)$ (for a different constant in the $\exp()$ function).

Proof of Theorem 1.1: By definition of the AdaptiveCompress protocol, we have $\|A-C\|=\|A-AUVV^{T}U^{T}\|,$ where all norms in this proof are the Frobenius norm.

Notice that $UU^{T}$ and $I_{d}-UU^{T}$ are projections onto orthogonal subspaces, where $I_{d}$ is the $d\times d$ identity matrix. It follows by the Pythagorean theorem applied to each row that

where the second equality uses that $U^{T}U=I_{c}$, where $c$ is the number of columns of $U$.

Observe that the row spaces of $AUVV^{T}U^{T}$ and $AUU^{T}$ are both in the row space of $U^{T}$, and therefore in the column space of $U$. It follows that since $U$ has orthonormal columns, $\|AUVV^{T}U^{T}-AUU^{T}\|=\|(AUVV^{T}U^{T}-AUU^{T})U\|,$ and therefore

where the second equality uses that $U^{T}U=I_{c}$. Let $(AU)_{k}$ be the best rank-$k$ approximation to the matrix $AU$. By Theorem 2.1, with probability $1-o(1)$, $\|AUVV^{T}-AU\|^{2}\leq(1+O(\varepsilon))\|(AU)_{k}-AU\|^{2},$ and so

Notice that the row space of $(AU)_{k}$ is spanned by the top $k$ right singular vectors of $AU$, which are in the row space of $U$. Let us write $(AU)_{k}=B\cdot U$, where $B$ is a rank-$k$ matrix.

We apply the Pythagorean theorem to each row in the expression in (5), noting that the vectors $(B_{i}-A_{i})UU^{T}$ and $A_{i}UU^{T}-A_{i}$ are orthogonal, where $B_{i}$ and $A_{i}$ are the $i$-th rows of $B$ and $A$, respectively. Hence,

and the last equality uses the definition of $B$. By Theorem 2.2, with constant probability arbitrarily close to $1$, we have

It follows by combining (2), (3), (4), (5), (6), (7), that $\|AUVV^{T}U^{T}-A\|^{2}\leq(1+O(\varepsilon))\|A_{k}-A\|^{2}$, which shows the correctness property of AdaptiveCompress.

We now bound the communication. In the first step, by Theorem 2.2, $m$ can be set to $O(k/\varepsilon)$ and the matrix $S$ can be described using a random seed that is $O(k)$-wise independent. The communication of steps 1-3 is thus $O(sdk/\varepsilon)$ words. By Theorem 2.1, the remaining steps take $O(s(k/\varepsilon)^{2}/\varepsilon^{2})=O(sk^{2}/\varepsilon^{4})$ words of communication.

To obtain communication with $O(b+\log(nd))$-bit words if the entries of the matrices $A_{t}$ are specified by $b$ bits, Server 1 can instead send $SA$ to each of the servers. The $t$-th server then computes $PA_{t}(SA)^{T}$ and sends this to Server 1. Let $SA=RU^{T}$, where $U^{T}$ is an orthonormal basis for the row space of $SA$, and $R$ is an $O(k/\varepsilon)\times O(k/\varepsilon)$ change of basis matrix. Server 1 computes $\sum_{t}PA_{t}(SA)^{T}=PA(SA)^{T}$ and sends this to each of the servers. Then, since each of the servers knows $R$, it can compute $PA(SA)^{T}(R^{T})^{-1}=PAU$. It can then compute the SVD of this matrix, from which it obtains $VV^{T}$, the projection onto its top $k$ right singular vectors. Then, since Server $t$ knows $A_{t}$ and $U$, it can compute $A_{t}U(VV^{T})U^{T}$, as desired. Notice that in this variant of the algorithm what is sent is $SA_{t}$ and $PA_{t}(SA)^{T}$, which each can be specified with $O(b+\log(nd))$-bit words if the entries of the $A_{t}$ are specified by $b$ bits.

2 Lower bound for low-rank approximation

Our reduction is from the multiplayer SUM problem.

() Suppose each of $s$ players has a binary vector $a^{i}$ with $n$ bits and the first player wants to compute $\sum_{i=1}^{s}a^{i}$ mod $2$ with constant probability. Then the total communication needed is $\Omega(sn)$ bits.

(of Theorem 1.2.) We reduce from the $s-2$ player SUM problem, in which each player has a $k\times d$ binary matrix $A^{i}$ and the first player wants to learn their sum. By Theorem 2.3, this problem needs $\Omega(skd)$ communication, since even the mod $2$ version of the problem requires this amount of communication. Now consider the $s$-player problem $s$-RESTRICT-SUM in which the first player has $I_{d}$, the second player has $-I_{d}$, the remaining $s-2$ players have a $k\times d$ binary matrix $A^{i}$ and the first player wants to learn the sum of all inputs. This also requires $\Omega(skd)$ communication. This follows since if this problem could be solved with $o(skd)$ communication, then SUM with $s-2$ players would have $o(skd)$ communication by a simulation in which the first player of the $(s-2)$-SUM problem simulates the first three players of the $s$-RESTRICT-SUM problem.

In our $s$-player low-rank approximation problem, we give the first player $I_{d}$, the second player $-I_{d}$, and the remaining $s-2$ players each has a random $k\times d$ binary matrix $A^{i}$. Note that there is a unique rank-$k$ approximation to the sum of the $k$ player inputs, namely, it is the matrix $\sum_{i=3}^{s}A^{i}$. It follows that any algorithm which outputs a projection matrix $VV^{T}$ for which $\|A-AVV^{T}\|_{F}^{2}\leq(1+\varepsilon)\min_{X:\textrm{ rank}(X)\leq k}$, for any $\varepsilon\geq 0$, must be such that $VV^{T}$ is a projection onto the row space of $\sum_{i=3}^{s}A^{i}$. This follows because $(1+\varepsilon)\min_{X:\textrm{ rank}(X)\leq k}\|A-X\|_{F}=(1+\varepsilon)\cdot 0=0$.

Now, since the first player has $I_{d}$, his output is $I_{d}VV^{T}$, where the row space of $V^{T}$ equals the row space of $\sum_{i=3}^{s}A^{i}$. Suppose the total communication of our problem is $C$.

We use this to build a protocol for $s$-RESTRICT-SUM, which has the same inputs as in our $s$-player low rank approximation problem. Notice that $A=\sum_{i=3}^{s}A^{i}$ is a $k\times d$ matrix with rows in $\{0,1,2,...,s-2\}^{d}$.

Claim. The span of the rows of $A$ can intersect $\{0,1,\ldots,s-2\}^{d}$ in at most $(2s)^{k}$ distinct points.

Let ${\mbox{r}owspace}(A)$ denote the row space of $A$. We will bound the size of ${\mbox{r}owspace}(A)\cap GF(p)^{d}$ for prime $p$ with $s-2<p<2(s-2)$, where $GF(p)$ is the finite field with elements $\{0,1,2,\ldots,p-1\}$, and $GF(p)^{d}$ is the vector space over $GF(p)$. This will be an upper bound on the size of ${\mbox{r}owspace}(A)\cap\{0,1,\ldots,s-2\}^{d}$. Since ${\mbox{r}owspace}(A)$ is $k$-dimensional, so is ${\mbox{r}owspace}(A)\cap GF(p)^{d}$. Hence the intersection has at most $k$ linearly independent points. These $k$ linearly independent points can be used to generate the remaining points in ${\mbox{r}owspace}(A)\cap GF(p)^{d}$. The number of distinct combinations of these points is at most $p^{k}<(2s)^{k}$, bounding the intersection size.

Next, players $3,\ldots,s$ agree on random $\{+1,-1\}^{d}$ vectors $u^{1},\ldots u^{k^{\prime}}$ where $k^{\prime}=k\log 4s$ via a public coin. The entries of $u^{1},...,u^{k^{\prime}}$ need only be $O(k\log s)$-wise independent, and as such can be agreed upon by all the players using only $O(sk\log s)$ bits of communication. Each player $i$ then computes the inner products $A^{i}_{j}\cdot u^{1}$, …, $A^{i}_{j}\cdot u^{k^{\prime}}$ for each $j\in\{1,\ldots,k\}$. Here $A^{i}_{j}$ denotes the $j$’th row of a the $i$’th player’s matrix $A^{i}$.

Frequency Moments and Higher Order Correlations

In this section, we prove Theorems (1.3), (1.4) and (1.6).

We begin with some common notation. For $t\in[s],i\in[n]$:

Let $B=\sum_{i}B_{i}=\sum_{t}C_{t}\;;\;A=\sum_{i}A_{i}$. The task is to estimate $A$. We analyze the following algorithm. Let $l=100\frac{s\cdot c_{f,s}}{\varepsilon^{2}}$. The parameters in the algorithm will be specified presently.

To analyze the algorithm, we think of it differently: suppose CP picks $t$ for the first of its $l$ trials and asks that $t$ to pick $i$ according to its $\{f(a_{ti}/C_{t}\}$. Let $X$ be the random variable $BA_{i}/B_{i}$ for that $i$. Clearly the estimate made by the algorithm can be viewed as the average of $l$ i.i.d. copies of $X$. So it will suffice to show that (i) $X$ is unbiased : I.e., $E(X)=A$ and (ii) Var$(X)\leq c_{f,s}sA^{2}$ (whence, the variance of the average of $l$ i.i.d. copies of $X$ would have variance at most $\varepsilon^{2}A^{2}$ giving us the relative error bound.)

The first part is easy: Let $p_{i}$ be the probability that we pick $i$ by this process. Clearly, p_{i}=\sum_{t=1}^{s}{\sf Prob}(\text{ CP pickedt}){\sf Prob}(t\text{ picks }i)=\sum_{t}\frac{C_{t}}{B}\frac{f(a_{ti})}{C_{t}}=\frac{B_{i}}{B}. So, $E(X)=\sum_{i=1}^{n}B\frac{A_{i}}{B_{i}}\frac{B_{i}}{B}=A,$ proving (i). For (ii), we have $E(X^{2})=B^{2}\sum_{i}p_{i}\frac{A_{i}^{2}}{B_{i}^{2}}=B\sum_{i}\frac{A_{i}^{2}}{B_{i}}\leq ABc_{f,s}\leq c_{f,s}sA^{2},$ since, $A_{i}=f(\sum_{t}a_{ti})\leq c_{f,s}\sum_{t=1}^{s}f(a_{ti})$ by the definition of $c_{f,s}$ and by monotonicity of $f$, we have $B_{i}=\sum_{t}f(a_{ti})\leq sf(\sum_{t}a_{ti})$.

To prove the claimed resource bounds, note that polynomial time and linear space bounds are obvious, since, all that each server has to do is to compute all $f(a_{ti})$, sum them up and sample at most $l$ times. The communication is dominated by each of $s$ servers sending $\{a_{ti},i\in S\}$ to CP which is $sc_{f,s}/\varepsilon^{2}$ words per server giving us a total of $O(s^{2}c_{f,s}/\varepsilon^{2})$.

Now for the lower bound, we use (rather unsurprisingly) the set-disjointness problem. It is known ( ) that the following problem needs $\Omega(n)$ bits of communication even for a randomized algorithm: we distinguish between two situations: (a) Each of $s$ servers holds a subset of $[n]$ and the subsets are pairwise disjoint and (b) There is exactly one element common to all $s$ sets. We reduce this problem to ours. Let $S_{t}$ be the subset held by server $t$. By definition of $c_{f,s}$, there exist $x_{1},x_{2},\ldots,x_{s}\in{\bf R}_{+}$ such that $f(x_{1}+x_{2}+\cdots+x_{s})=c_{f,s}(f(x_{1})+f(x_{2})+\cdots+f(x_{s})).$ Let $n=\frac{c_{f,s}-1}{\varepsilon}.$ Let $a_{ti}$ be defined by: $a_{ti}=x_{t}\text{ if }i\in S_{t}$ and $a_{ti}=0$ otherwise. If the sets are disjoint, then $\sum_{i=1}^{n}f\left(\sum_{t=1}^{s}a_{ti}\right)=\sum_{t=1}^{s}|S_{t}|f(x_{t}).$ In the case (b) when the sets all share one element in common, $\sum_{i=1}^{n}f\left(\sum_{t=1}^{s}a_{ti}\right)=\sum_{t=1}^{s}(|S_{t}|-1)f(x_{t})+f(x_{1}+x_{2}+\cdots+x_{s})=\sum_{t=1}^{s}|S_{t}|f(x_{t})+(c_{f,s}-1)\sum_{t}f(x_{t})=\sum_{t=1}^{s}|S_{t}|f(x_{t})+\varepsilon n\sum_{t}f(x_{t}).$ Since $|S_{t}|\leq n$, it follows that if we can estimate $\sum_{i}f(\sum_{t}a_{ti})$ to relative error $\varepsilon$, then we can distinguish the two cases. But it is known that this requires $\Omega(n)$ bits of communication which is $\Omega(c_{f,s}/\varepsilon)$ proving the lower bound.

(of Theorem (1.4): The only change is in the sampling algorithm:

Order the $j=(j_{1},j_{2},\ldots,j_{k})$ lexicographically. Start with the first $j$ as the sample and compute $a_{tj}$ by making a pass through the entire data: For each $i\in W_{t}$, after $v_{i,j_{1}},v_{i,j_{2}},\ldots,v_{i,j_{k}}$ are read, compute $g(v_{i,j_{1}},v_{i,j_{2}},\ldots,v_{i,j_{k}})$ and sum over all $i\in W_{t}$.

Process the next $j$ similarly. After processing a $j$, say, $j=j_{0}$, compute $f(a_{tj_{0}})$ and keep a running total of $f(a_{tj})$ for all $j$ seen so far. Reject the old sample and replace it by the current $j_{0}$ with probability $\frac{f(a_{tj_{0}})}{\text{ Total of all }f(a_{tj})\text{ including }j_{0}}.$

If the old sampled $j$ is not rejected, just keep it as the sample and go to next $j$.

The proof of correctness and linear space bound follow straightforwardly by plugging in this sampling algorithm into Theorem (1.3).

We next turn to a more refined algorithm for estimating frequency moments with near-optimal communication, using the specific function $f(x)=x^{k}$. Here is the algorithm.

(of Theorem (1.6): Let $B_{i}=\sum_{t=1}^{s}a_{ti}^{k}$ and $A_{i}=\left(\sum_{t=1}^{s}a_{ti}\right)^{k}$ and $\rho_{i}=A_{i}/B_{i}$. Note: $1\leq\rho_{i}\leq s^{k-1}$. CP can arrange to pick $i\in[n]$ with probabilities $\{B_{i}/B\}$, where, $B=\sum_{i}B_{i}$ as we already saw. First pick $m=s^{k-2}/\varepsilon^{3}$ sample $i\in[n]$ according to $\{B_{i}/B\}$. Then, CP tells all servers all these $i$ and collects all $a_{ti}$ and thence all $A_{i},B_{i}$. Total communication is at most $ms\leq s^{k-1}/\varepsilon^{3}$.

Let $\rho_{i}=A_{i}/B_{i}$. Let $\beta$ range over $\{s^{k-1},e^{-\varepsilon}s^{k-1},e^{-2\varepsilon}s^{k-1},\ldots 1\}$, a total of $O(\ln s)$ values and let $S_{\beta}=\{i\in S:\rho_{i}\in[\beta e^{-\varepsilon},\beta)\}.$ Then, $\sum_{i\in S}\rho_{i}\approx\sum_{\beta}|S_{\beta}|\beta.$ Since each $\rho_{i}\geq 1$, we have $\sum_{i\in S}\rho_{i}\geq|S|$. So we need only accurate estimates of those $|S_{\beta}|$ with $|S_{\beta}|\geq\varepsilon|S|/\beta\ln s$. For each $S_{\beta}$, if we pick u.a.r. a subset $T$ of $S$ of cardinality $\Omega(\beta(\ln s)^{2}/\varepsilon^{3})$, then $\frac{|S|}{|T|}|T\cap S_{\beta}|$ estimates $|S_{\beta}|$ to within $(1\pm\varepsilon)$ for every $\beta$ satisfying $|S_{\beta}|\geq\varepsilon|S|/\beta\ln s$.

For each $\beta$, pick such a random subset $T$ from $S$. We have to recognize for each $i\in T$, whether it is in $S_{\beta}$. First, for each $i\in T$, we estimate $A_{i}$ as follows: We pick $l=s^{k-1}/\beta$ servers $t_{1},t_{2},\ldots,t_{l}$ u.a.r. and take $Z_{i}=\frac{s}{l}(a_{t_{1},i}+a_{t_{2},i}+\cdots a_{t_{l},i})$ as our estimate of $A_{i}^{1/k}=\sum_{t=1}^{s}a_{ti}$. If $Y$ is the r.v. based on just one random server (namely $Y=sa_{ti}$ for a u.a.r $t$), then ${\sf E}Y=A_{i}^{1/k}$ and

From this it follows by averaging over $l$ samples that

We do $\Omega(k\ln s+\ln(1/\varepsilon))$ such experiments and take the median of all of these to drive down the failure probability for a single $i$ to less than $\varepsilon^{2}/s^{k}$, whence, it is small by union bound for the failure of any of the at most $s^{k-1}/\varepsilon^{2}$ indices $i$ ’s. Thus all the $A_{i}^{1/k},i\in T$ are estimated to a factor of $(1+\epsilon)$ by this process whp. Since $k$ is a fixed constant, this also means $(1+O(\varepsilon))$ relative error in the estimate of $A_{i}$.

The subset that intersects $T$ is then $\{i\in T:\rho_{i}\geq\beta/s\}|\leq 20\frac{|T|}{|S|}\frac{s^{2}\ln s|S|}{\beta}=O(s^{2}\ln^{3}s/\varepsilon^{3})$. Now for these $i$ ’s in $T$, we collect all $a_{ti}$ and find $A_{i},B_{i}$ exactly. This costs us $O(s^{3}(\ln s)^{3}/\varepsilon^{2})$ communication. Thus the overall communication is bounded by

We now extend the above theorem and proof for a wide class of functions satisfying a weak Lipschitz condition (and generalizing the case of moments).

For a monotone function $f:\Re_{+}\rightarrow\Re_{+}$, define

Alternatively, $L_{f}$ is the Lipschitz constant of $f$ wrt the “distance” $d(x,y)=\log(x)-\log(y)$, i.e.,

For the function $f(x)=x^{k}$, we see that $L_{f}=k$.

For any function $f:\Re_{+}\rightarrow R_{+}$ with $L=L_{f}$,

Let $f$ be any nonnegative, superlinear real function with $L=L_{f}\geq 4$. Suppose there are $s$ servers, with server $t$ holding a non-negative vector $a_{t}=(a_{t1},a_{t2},\ldots,a_{tn})$ Then, to estimate $A=\sum_{i=1}^{n}f\left(\sum_{t=1}^{s}a_{ti}\right)$ to relative error $\varepsilon$, there is an algorithm that communicates $O(s^{L-1}(\ln s)^{3}/\varepsilon^{3})$ words in $O(1)$ rounds.

The algorithm is the following, essentially the same as in the moments case, with parameters defined in terms of $L_{f}$ for a general function $f(.)$ in place of $x^{k}$.

Pick an i.i.d. sample $S_{0}$ of $m=s^{L-2}/\varepsilon^{3}$ indices $i$, where each $i$ is picked according to $\{B_{i}/B\}$.

CP now gets an i.i.d. sample $S$ of $O(s^{L-1}\ln^{2}s/\varepsilon^{3})$ $i$ ’s, each according to $\{B_{i}/B\}$.

For each $\beta\in\{s^{L-1},e^{-\varepsilon}s^{L-1},e^{-2\varepsilon}s^{L-1},\ldots,1\}$, CP does the following:

Pick a subset $T$ of $S$ of cardinality $\Omega(\beta(\ln s)^{2}/\varepsilon^{3})$ u.a.r.

(of Thm. 3.3.) We point out the changes in the analysis from the special case of moments. Now we have $A_{i}=f(\sum_{t=1}^{s}a_{ti}),B_{i}=\sum_{t=1}^{s}f(a_{ti})$ and $A=\sum_{i=1}^{n}A_{i},B=\sum_{i=1}^{n}B_{i}$. Also $\rho_{i}=A_{i}/B_{i}$. The reader will have noticed that $k$ has been replaced by $L$ in the above algorithm. The first phase remains the same, and at the end we either have a good approximation for $A$ or we know that $A\leq sB$.

In the next phase we estimate $f(\sum_{t}a_{ti})^{1/L}$. To do this, we first estimate $\sum_{t=1}^{s}a_{ti}$, then apply $f$ to this estimate. We need to analyze the error of both parts. For the first part, let $Y=sa_{ti}$ as before. Then ${\sf E}(Y)=\sum_{t=1}^{s}a_{ti}$ and since the server used to define $Y$ is chosen uniformly at random, we have

The last phase for estimating $B_{i}$ and putting together the estimates for all the $\rho_{i}$ is again the same as in the case of moments.

Acknowledgements. We are grateful to Dick Karp and Jelani Nelson for helpful discussions, as well as the Simons Institute at Berkeley. Santosh Vempala was supported in part by NSF award CCF-1217793 and David Woodruff by the XDATA program of the Defense Advanced Research Projects Agency (DARPA), administered through Air Force Research Laboratory contract FA8750-12-C0323.

References