Matrix Sketching

Last time, I wrote about a clever algorithm for approximating the histogram for a stream using bounded memory. The post was motivated by this paper, which is an extension of that algorithm to a problem that seems unrelated at first glance, which is matrix sketching. The matrix sketching problem is as follows: given the rows of a large matrix $A \in \mathbb{R}^{n \times m}$ as a stream, produce a matrix $B \in \mathbb{R}^{l \times m}$ where $l << n$ which is a "good" approximation for $A$ when multiplied with vectors. Specifically, the algorithm in the paper achieves the following result: if $l = 1/\epsilon$, then it produces a matrix $B$ such that, for any unit vector $x$,

$||Ax||^2 \ge ||Bx||^2 \ge ||Ax||^2 - \epsilon ||A||_f^2$

where $||A||_f$ is the Frobenius norm of $A$. Here you can see the parallels with the frequency approximation algorithm; the error is a function of how "big" the stream is, which in this case is the Frobenius norm of the input matrix.

The algorithm works as follows: start with $B$ as a $l \times m$ matrix of all zeroes, and for each input row $A_i$ do the following update:

1. Set $B_l = A_i$ (the last row of $B$).
2. Compute the singular value decomposition (SVD) of $B$, so we obtain $B = U \Sigma V$ with the standard assumption that the diagonal values of $\Sigma$ are $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_l$.
3. "Remove" the smallest singular value from $\Sigma$ by letting
   
   $\bar{\Sigma} = \sqrt{\max(\Sigma - I_l \sigma_l^2, 0)}$
   
   where $I_l$ is the $l \times l$ identity matrix.
4. Then set $B = \bar{\Sigma}V$ (note that the last row of $B$ is all zeroes after this step because the last row of $\bar{\Sigma}$ is all zeroes by construction).

At the end, just output the current value of $B$. I won't go into any of the proof details (they can be found in the paper), but it's interesting to try to understand what the algorithm is doing intuitively. The SVD can be thought of (very loosely) as breaking down a matrix into three transformations applied to a multidimensional space: a rotation ($V$), followed by a scaling along the axes ($\Sigma$), and lastly another rotation ($U$). So the singular values are the scaling factors of the matrix in orthogonal directions, and we are removing the smallest one from each of these directions equally. As a result, we only lose a fraction of the accuracy in any particular direction (i.e. $Bx$ versus $Ax$ for a specific $x$) compared to the Frobenius norm of $A$ as a whole.

This would be pretty cool even if it was a purely theoretical result since it ties two problems together in a very elegant way, but it gets better. The paper goes on to explore the algorithm's performance with some experiments and observes that the accuracy is quite an improvement over existing techniques for matrix sketching. Moreover, computing the SVD is somewhat expensive, so the author describes how the algorithm can be parallelized as well as a way to reduce the computational cost of the SVD step and only slightly relaxing the guarantees. It's a very nice paper that spans both the theoretical and practical domains for the matrix sketching problem.

Streaming Frequency Approximation

It's been a while, but I finally motivated myself to write another blog post! Today's topic is approximating the frequency of items in a stream. With all the data that is generated today, streaming algorithms have become very popular as one of the efficient ways to process datasets that are far too large to fit in memory. One common data mining problem is determining the frequency of items in data, i.e. building a histogram. Consider a stream $X = x_1, x_2, \ldots, x_n$ of $n$ items, each of which belongs to one of $m$ classes $c_1, c_2, \ldots, c_m$. If $m$ is small, the problem is simple, as we simple maintain the entire histogram in $O(m)$ memory and output the exact result. The interesting case is when there are an extremely large number of distinct classes, the number of which is usually unknown. In this case, we will present an algorithm originally from this paper that approximates each value in the histogram within $\epsilon n$ using $O(1/\epsilon)$ memory.

At any point while processing the stream, the in-memory state is a partial histogram with at most $k = 1/\epsilon$ elements (class/count pairs). When processing an element $x_i$, we use the following procedure with three cases:

1. $x_i$ is already in the histogram: increment its count;
2. we have fewer than $k$ elements in the histogram: add $x_i$ with count $1$;
3. the histogram already has $k$ other elements: while the histogram has $k$ elements, decrease the count of each element by $1$ and evict any elements with zero count, and then add $x_i$ with count $1$.

That's the entire algorithm, and it's easy to prove that it achieves the desired approximation as well. The only case in which we "lose" information is case 3, where the count of each of $k$ elements is decreased by $1$. But we can clearly decrease counts at most $n/k = \epsilon n$ times, since that would necessarily reduce all counts to zero. Thus, using the partial histogram and outputting the frequency of any class not in the histogram as zero, we are guaranteed that each count is within $\epsilon n$ of the actual count. In particular, this means that the algorithm will be guaranteed to output a positive count for any class that occurs with frequency $\epsilon n$ or more in the stream.

I've always found this problem to be particularly beautiful because of the simplicity of the result, the algorithm, as well as the proof. Moreover, it's a very practical problem with many applications. To make things even more interesting, a recent paper shows how this concept can be generalized to do matrix sketching (i.e. approximating a matrix with a smaller one). I hope to get around to writing about that at some point as well.

Dremel Data Model

It is common knowledge that analyzing large datasets efficiently can benefit greatly from column-oriented storage. That is, instead of storing all the data for a single row together like a traditional database would, separate the columns out and store all of the data for each column together (see the Wikipedia link for an example). The benefits of this are twofold: (1) queries that only access a subset of the columns can reduce the amount of data that has to be retrieved, and (2) compression algorithms often perform better on homogeneous data, e.g. a column with many values of the same type. As such, data stores such as Cassandra and HBase have become widely used when dealing with records that have a large number of columns. Given that both of these projects were based on the original ideas of Google's BigTable, it should be no surprise that Google has continued leveraging the power of column-oriented storage to deal with the scale of data that it processes. Dremel is a project that has been used at Google for a number of years now; it takes the idea of column-oriented storage, generalizes it, and then optimizes it in order to perform queries over tens to hundreds of terabytes of data in seconds.

The key to Dremel's performance is how the data is represented and consequently laid out on disk. First of all, they generalize having a bunch of columns per record to a hierarchical structure of nested columns that supports repeated and optional fields. That may sound familiar because it is exactly the format of protocol buffers, which is the language-independent binary representation of data that Google generally uses, so it is a natural choice for the type of data Dremel should support. Column-oriented storage is simple in the normal case of one (potentially optional) value per column per record since you can easily figure out which value corresponds to which record, but it gets trickier when you allow for the full protocol buffer specification. Consider the following example schema (a simplified version of what is presented in the paper):


This represents a web document, which has an ID and a set of names associated with it; each name is a URL that points to that document and the languages associated with the URL. The generalization of column-oriented storage to handle nested columns is that each leaf field in the hierarchy is stored separately. For example, all values of "docId" are stored together, and all values of "name.language.country" are stored together since both are leaf fields in the schema for a document. Let's look at two example records which we will assume are stored in the order presented:


The "docId" field is simple and only needs to store the values "10" and "20" since it is a top-level, required field. The "name.language.country" field, on the other hand, will need to store the values "us", NULL, NULL, "gb", NULL. We will see why these NULL values are necessarily shortly, but they are basically placeholders for potential values of "name.language.country."

It should be clear that because of the repeated and optional values we cannot only store the values above without additional metadata, since doing so would lead to the loss of record boundaries. To solve this, Dremel introduces two metadata fields known as the "repetition level" and the "definition level." These are logically attached to every column value in order to preserve record boundary information. The repetition level handles repeated fields and indicates the depth in the hierarchy that was repeated between the current column value and the previous one (0 is the root, meaning a new document). The definition level handles optional fields and indicates the depth in the hierarchy that actually exists for the current value; this is only relevant for NULL values. For the "name.language.country" field, Dremel would store the following:

#	value	rep	def
1	"us"	0	3
2	NULL	2	2
3	NULL	1	1
4	"gb"	1	3
5	NULL	0	1

In row 1 of the table, rep = 0 because it is a new record, and def = 3 because all 3 levels of the hierarchy ("name", "language", and "country") are defined, which is why we have a non-NULL value. In row 2, rep = 2 because we are repeating "name.language" when going from the previous row to this one, and def = 2 because only "name" and "language" are defined so the value is NULL. In row 3, rep = 1 because we repeated at the "name" level, and def = 1 because only "name" is defined. In row 4, rep = 1 because we again repeated at the "name" level, and def = 3 because we have a real value. Lastly, row 5 has rep = 0 because it is a new record, and def = 1 because only "name" is defined. The meanings of the repetition and definition levels can be confusing at first, but if you understand how they are derived in the above table, you should be able to convince yourself that it is a lossless representation. Dremel then optimizes to the bit-level how much data they need to actually store, such as omitting representation and definition levels whenever possible and using as few bits as necessary to store the values.

It should be no surprise that Google's solution for interactive queries on huge datasets involves minimizing the amount of data that needs to be retrieved from disk, a concept that database engineers have employed for decades. Since the volume of data being collected today far outpaces amount that can be read off of disk in a few seconds, there seem to be few options in the analytics space for efficient queries other than designing to reduce what needs to be read. But one can certainly imagine more complex queries that involve a large number of fields or even joining two fields that would be extremely valuable yet essentially impossible given the limitations of today's technology. This area seems ripe for additional work in pushing the boundaries of what types of queries are possible on the huge datasets that companies are collecting, but Dremel is definitely a solid approach that Google has derived a lot of value from.

Anonymized Social Network Data

Social network analysis is a popular research area these days thanks to the rapid growth of networks like Facebook, Google+, and Twitter. Data for subsets of these networks is readily available, e.g. here, for researchers to analyze and discover interesting properties. Sometimes, these datasets are anonymized to protect the privacy of users while preserving the structure of the network; for example, instead of nodes being labelled with names, they are often replaced with a random ID. An interesting question that arises is whether this is secure, that is, can an attacker extract information from data anonymized in this way? Several years ago, it was shown that individuals could be identified from the anonymized Netflix Prize dataset, so this question is important for the research community to consider. In this paper, it is shown that an active attacker can quite easily defeat the naive random ID anonymization.

Let $G$ be the anonymized network with $n$ nodes that an attacker has access to. In this context, an active attacker is someone who has "user-level access" to the network prior to the creation of the anonymized dataset. For example, if $G$ is the anonymized Facebook network, an attacker can create accounts and friend people arbitrarily before $G$ is created such that those nodes and edges will be reflected in $G$. The resulting attack is quite simple. Let $H = \{ x_1, x_2, \ldots, x_k \}$ be the set of "fake" nodes created by the attacker and $W = \{ w_1, w_2, \ldots, w_b \}$ be the set of existing nodes that the attacker wants to compromise (in particular, learn about edges between $w_i$ and $w_j$). It is shown that with just $k = \Theta( \log{n} )$ attacker nodes, you can, with high probability, perform an attack efficiently on $b = \Theta( \log^2{n} )$ of the "real" nodes in the network. The attack is based on generating $H$ in such a way that it can easily be recovered from $G$ and connecting $H$ to $W$ so that each node in $W$ can be uniquely identified.

There are a few properties of $H$ that need to be guaranteed. There cannot be another induced subgraph in $G$ which is isomorphic to $H$, so that we can identify $H$, and there must be an efficient algorithm to find $H$ within $G$. To generate $H$, first include the edges $(x_i, x_{i+1})$ for $i = 1, 2, \ldots, k-1$ and then include every other edge with probability $1/2$. Then each $w_j$ is connected to a small, unique subset of nodes in $H$ so that we can recover all of the $w_j$'s once we find $H$. Furthermore, each $x_i$ has edges added randomly to nodes in $G - H$ until it reaches degree $\Delta_i = O(\log{n})$. Based on the structure of $H$ and the degrees, we can now recover it from $G$ using a search over short walks in $G$. Starting from each node, consider it as a potential $x_1$ and then, if that passes the degree requirements, consider all of its neighbors as potential $x_2$'s. For each neighbor that passes the degree requirements and the internal structure constraints of $H$, continue the walk, and repeat this until a walk with $k$ nodes is found. Through some delicate analysis, the authors show that, with high probability, there will be a unique walk that satisfies all the constraints (corresponding to $H$ itself) and that the search tree will have size $O(n^{1+\epsilon})$ so the algorithm is tractable on large datasets.

The paper describes experiments run on a LiveJournal dataset with over 4 million nodes and 70 million edges. With just $k = 7$ attacker nodes, they are able to successfully recover $H$ with roughly 90% probability, and it only goes up as $k$ increases. This means that it is simple to perform such attacks in practice and, perhaps more importantly, that the attacks have a sufficiently small effect in the network to avoid detection. The end result of the work is captured quite nicely by a quote from the paper, highlighting the observation that simple anonymization is not enough: "true safeguarding of privacy requires mathematical rigor, beginning with a clear description of what it means to compromise privacy, what are the computational and behavioral capabilities of the adversary, and to what information does it have access."

Collaborative Filtering

Collaborative filtering is a very popular technique being employed these days by recommendation systems such as those powering Amazon or Netflix. The idea is using existing user-item ratings (e.g. for movies or products) and learned relationships between users and items to predict what other items users might be interested in. It turns out that collaborative filtering is more accurate and scalable than other methods, such as profile-based recommendations, i.e. trying to recommend based on what the users describe their profile as. The most famous dataset for collaborative filtering comes from the Netflix Prize competition, which provided a set of 100 million user-movie ratings over a span of about six years; the goal was to have the most accurate predictions on a set of 1.4 million ratings that were not released. Although Netflix has evolved as a company since the competition, the dataset led to an incredible wealth of research in collaborative filtering which has certainly benefited the many other recommendation systems out there.

One interesting paper in this field that came from the winning team described how to integrate temporal dynamics into the predictions. It observes that the time when a user rates a movie is actually crucial for prediction. This piece of data was often ignored for the first 1-2 years of the content, but the methods used are able to capture real-life effects very well. To describe how temporal dynamics are incorporated, we'll first need some background on one of the most used modeling techniques for the Netflix Prize, latent factor models (also known as factor analysis or matrix factorization).

Consider the rating matrix $M$ where $M_{ui}$ is the rating that user $u$ would give to movie $i$, ignoring the time component for now. We assume that users and movies are characterized by $f$ latent factors $p_u, q_i \in \mathbb{R}^f$ so that $M_{ui} = \langle p_u, q_i \rangle$. For example, we can think of the $k$-th component of these factors as how much a user likes action movies and how much action a movie has, respectively. Thus when $p_{uk}$ and $q_{ik}$ correlate positively the rating will be be higher than when they do not. Note that the training set $S = \{ r_{ui}(t) \}$ (user-movie ratings with times) provided by Netflix is a very sparse subset of the entries in $M$ so the goal is to use that sparse subset to interpolate the rest of $M$, which in the latent factor model amounts to learning the appropriate $p_u$ and $q_i$ from the training set. This can be done through stochastic gradient descent on the regularized objective function to solve for

$$\min_{q, p} \sum_{r_{ui} \in S} (r_{ui} - \langle p_u, q_i \rangle)^2$$
There are a few additional details, including baseline predictors (to account for certain users or movies to have intrinsically higher or lower ratings) and regularization, but the crux of the predictive power comes from the latent factor representation. Once the $p$ and $q$ vectors are learned, prediction becomes trivial.

You'll notice, however, that the time domain is completely ignored in the described model. In fact, such a model makes the assumption that movie characteristics and user preferences do not change over time. While the former may be true, the latter most certainly is not, and accounting for that can improve upon the existing model greatly. So instead of $p_{uk}$ we have the time-dependent preference function $p_{uk}(t)$; by learning a function here, we allow for variation in a user's preferences over time. It's important, however, to have $p_{uk}(t)$ be structured in some way; otherwise, the number of variables is far too large for the size of the training set and we would most likely be overfitting. One of the time-based models from the paper is to let

$$p_{uk}(t) = p_{uk} + \alpha_{uk} \cdot \text{dev}_u(t)$$
where $\text{dev}_u(t) = \text{sign}(t - t_u) \cdot |t - t_u|^{\beta}$ and $t_u$ is the mean date of rating for user $u$. This assumes that a user's preference for one of the latent factors tends to either increase or decrease over time rather than oscillating. The $\alpha_{uk}$ and $\beta$ parameters determine the direction and shape of this variation. The authors use cross validation to choose the best value of $\beta$ (it turns out to be 0.4) and learn the parameters $p_{uk}$ and $\alpha_{uk}$. Since the model only has twice as many parameters per user than the non-temporal model, it is still efficient to learn. Again, there are additional details which can be found in the paper, but the main idea is to model variation in a user's preferences using the structured $p_{uk}(t)$ function.

Not only does adding temporal dynamics to the latent factor model increase the score on the Netflix Prize from 6.3% to 7.5% improvement over the Cinematch baseline (which, if you've worked on it at all, is an absolutely huge difference), but the author was able to draw several neat insights as a result of it. He was able to identify the potential cause behind the observed global increase in average ratings around 2004, which turns out to be improved recommendations by the existing Netflix system causing people to tend to watch movies that they are expected to enjoy. Moreover, he observed that the older a movie is, the higher its ratings tend to be, which is a product of the fact that people tend to select old movies more carefully than new ones. Temporal dynamics are very interesting because they are intuitive to people (as these observations show) but often difficult to model computationally. The work in the paper, however, shows that such dynamics are perhaps even more important than we would expect, and can be leveraged to significantly improve recommendation systems.