Visit to DIMA at Technische Universitaet

I had a great visit today at the DIMA laboratory at TU in Berlin.  They are working on an interesting system called Stratosphere which provides an interesting generalization generalization of map-reduce.  Of particular interest is the run-time flexibility for adapting how the flow partitions or transfers data.

They accomplish this by having a lower level abstraction layer that supports a larger repertoire of basic options beyond just map and reduce.  These operations include match, cross product and co-group.  Having a wider range of operations and retaining some additional flow information at that level allows them to do on-the-fly selection of the detailed algorithm for different operations based on the statistics of the  data and the properties of the user-supplied functions.

Here's a pic of me answering questions about startups and log-likelihood ratio tests.

A Tour of the Multi-update For Zookeeper

I have recently been working on some new features for Apache Zookeeper that represent probably the largest change in the Zookeeper client API in several years.  Other recent changes such as observers and read-only clusters have changed the server-side semantics, but the client API is nearly unchanged since the original public release of Zookeeper several years ago.  The JIRA covering the overall issue is Zookeeper-965 and you can look at the code as it is being refined into committable state in my fork of Zookeeper that I keep on github.

(related announcement)

The Problem

Almost from the beginning of the Zookeeper project, there have been repeated questions on the mailing about how to update several nodes at once.  The answer has always been to consolidate the structures that you need to update atomically into the contents of a single znode and then update those contents atomically using standard methods.  This update problem is often exacerbated by a desire to use ephemeral nodes so that state disappears correctly when a service crashes.

This is a pretty reasonable answer and it handle many important cases fairly well.  For instance, in the common case of servers that are assigned certain roles and then in turn advertise which roles they have successfully taken on, this pattern can be used by giving each server a single file of assignments and a single file of advertised capabilities.  Each of these files can be ephemerally created by the server so they will both vanish when the server disappears.  Assignments can be added atomically and only the server ever modifies the advertised roles.  Zookeeper works quite well in these cases.

Keeping track of a graph full of nodes with directional connections to other nodes is a good example of where Zookeeper's update model is not so nice.   In this case, nodes have a list of connections to other nodes and a list of connections from other nodes.  If any node has a list of outgoing connections that are not matched by the corresponding list of incoming connections on other nodes, then the graph is referentially inconsistent.  Keeping the graph consistent under changes is not possible with normal Zookeeper updates unless you store the entire graph in a single file.

Lazy Operations

The new multi() method allows such a data structure to be updated with strong guarantees of consistency.  This is done by currying the normal database mutation operators and then passing all of the resulting closures to Zookeeper at once for execution.  Obviously this requires a bit of syntactic sugar in languages like Java and C which don't like to express closures directly.

The way that this works is that there is a new class call Op.  There are sub-classes of Op called Op.Create, Op.SetData, Op.Check and Op.Delete that correspond to the operations  normally invoked by calling the similarly named methods on a ZooKeeper object.  In essence, these sub-classes of Op represent reified methods or closures that can be frozen at one point in time and then executed later.  These sub-classes are instantiated using factory methods on Op class.  The names and signatures of these factory methods mimic the corresponding methods on ZooKeeper.

Once you have all of the operations you would like to perform, you can execute them all in a single transaction using ZooKeeper.multi().  This will either execute all of the Op's or none of them.

An Example

I have placed an example program for doing a graph-based computation over on github.  This program builds a graph consisting of 32 nodes in the form of a 5-dimensional hyper-cube and then uses a numerical technique called over-relaxation to solve for the voltages that would result if all the links in the hyper-cube were replaced by unit resistors.  A picture of the graph is just to the right. In this graph, the voltage for node 0x1F is labeled as $V_5$ and the voltage for node 0x0 is labeled as $V_0$.  There are lots of connections and solving this problem analytically is difficult unless you twig to the trick of using symmetry.  Real-world networks generally don't exhibit such congenial properties and can be nearly impossible to solve analytically.

Normally, of course, we wouldn't actually use Zookeeper for this computation, but the point here isn't so much a realistic computation as a test-bed for the new multi() method and a demonstration of how the closure based style associated with multi-ops makes code easier to read and understand.

The Code

To read the code, start with the class ZGraphTest.  This is a JUnit test that drives all the rest of the code.  In this test, a graph is constructed (named zg in the code), nodes are connected in the pattern of a hyper-cube which means that nodes are labeled with integers from 0 to 31 and a node $n_1$ is connected to another node $n_2$ if $n_1$ and $n_2$ differ in exactly one bit.

In each iteration of the code, the ZGraph.reweight() method of the graph is called on a randomly selected node.  This sets the value at that node to be the average of the values at the neighbors of that node.  This computation converges geometrically to the correct value with errors decreasing by a factor of 5-10 with every 1000 iterations.  As the actual computation proceeds the error versus the theoretically known correct values is shown every 1000 iterations and  you can see the errors go to zero with 6 significant figures at about 7000 iterations.

Internally, a ZGraph keeps a connection to Zookeeper and the name of the root directory for all of the nodes in the graph.  Methods like connect(), reweight() and update() all work by defining lazy update or version check operations on Node objects and then executing all of the update or version checks at one time in a transaction.  For instance, in the reweight() method, this code gets the weight of all of the neighbors of node g1:

        double mean = 0;

        Set neighbors = g1.getIn();

        List ops = Lists.newArrayList();

        for (Integer neighbor : neighbors) {

          Node n = getNode(neighbor);

          ops.add(n.checkOp(root));

          mean += n.getWeight();

        }

As a side effect, we also collect a list of version check operations into the list ops.  Then in this code we add one more operation to set the weight on g1 and add the update operation to the list ops:

        // set up the update to the node itself

        g1.setWeight(mean / neighbors.size());

        ops.add(g1.updateOp(root));

Finally, multi() is called to actually do all of the version checks and updates that we have collected,

        zk.multi(ops);

The essential point here is how the list of operations was collected a bit at a time and then executed all at once.  The versions of the nodes that were inputs into the operation were collected in the form of check operations.  When those check operations are executed along with the update of g1, we guarantee that g1's new weight will accurately reflect the average of its neighbors and if somebody changes the value of any of the neighbors in between the time that we read the value and the time we actually do the update, we will run the update again.

Similarly, and closer to the idea of maintaining referential integrity, the ZGraph.connect() collects update operations for the two nodes being connected and executes both updates together to avoid the possibility that anyone would see a situation where one node connects to a second but the second doesn't have the reverse connection.  This code does the job,

        Node g1 = Node.readNode(zk, root, id1);

        Node g2 = Node.readNode(zk, root, id2);

        g1.connectTo(id2);

        g2.connectFrom(id1);

        zk.multi(Arrays.asList(g1.updateOp(root), g2.updateOp(root)), results);

Again, the closure passing style makes this operation very easy.

One additional point to be noted is that having each Node return closures for updates that can be executed in any context the caller would like has the effect of removing all direct Zookeeper operations from the Node other than for reading or creating a node.  It also removes all references to serialization of a Node from any outside code.  This simplifies both the caller and the Node because the Node doesn't have to implement all plausible combinations of multiple updates and the caller doesn't have to worry about any aspect of serialization and can focus on the appropriate task of arranging the transactional semantics for updates.

Credit Where Credit is Due

Of course, this new Zookeeper capability wouldn't be possible if the Apache Zookeeper project team hadn't built Zookeeper in a very flexible way in the first place.  Kudos to Ben Reed and Mahadev and the other early committers on the project.

Also, the actual Zookeeper-965 project would have stalled out if Marshall McMullen and Camille Fournier hadn't stepped in with a bit of extra momentum.  I had completed the wire protocols and all of the client side work, but Marshall did all of server side work and Camille provided really excellent advice about how the changes should be done.

The final credit goes to MapR Technologies for supporting open source by allowing capabilities like multi to be contributed back.

Online algorithms for boxcar-ish moving averages

One problem with exponentially weighted moving averages is that the weight that older samples decays sharply even for very recent samples.

The impulse response of such an average shows this clearly. In the graph to the right, the red line is the impulse response of the exponential weighted average is shown by the red line. The impulse response of a different kind of moving average derived by John Canny in the early 80's is shown by the black line.

The Canny filter puts higher weight on the events in the recent past which makes it preferable when you don't want to forget things right away, but do want to forget them before long and also want an on-line algorithm.

The cool thing about the Canny filter is that it is only twice as much work as a simple exponential moving average. The idea is that if you take the difference between two exponentially weighted averages with different time constants, you can engineer things to give you an impulse response of the sort that you would like.

The formula for such a difference looks like this
\begin{eqnarray*}
w(x) &=& k_1 e^{-x/\alpha} - k_2 e^{-x/\beta}
\end{eqnarray*}
Here $k_1$ and $k_2$ scale the two component filters in magnitude and $\alpha$ and $\beta$ are the time constants for the filters.

It is nice to set the magnitude of the filter at delay $0$ to be exactly 1. We can use this to get a value for $k_2$ in terms of $k_1$
\begin{eqnarray*}
w(0) &=& k_1 - k_2 = 1 \\
k_2 &=& k_1 -1
\end{eqnarray*}
Similarly, we can constrain the slope of the impulse response to be $0$ at delay $0$. This gives us $\beta$ in terms $\alpha$
\begin{eqnarray*}
w'(0) &=& {k_1 \over \alpha} - {k_2 \over \beta} = 0\\
{k_1 \over \alpha} &=& {k_1-1 \over \beta} \\
\beta &=& \alpha \frac{ k_1-1} { k_1}
\end{eqnarray*}
The final result is this impulse response
\begin{eqnarray*}
w(x) = k \exp \left(-{x \over \alpha}\right)-(k-1) \exp\left(-{k x\over \alpha (k-1)}\right)
\end{eqnarray*}
We can do a bit better if we note that as $k \to \infty$, the shape of the impulse quickly converges to a constant shape with mean of $w(x) \to \frac{3a}{2}$ and a total volume of $2a$.

Exponentially weighted averaging for rates

In a previous post, I talked about how to produce exponentially weighted averages of time-embedded values.  This is fine for cases where you really want averages, but it doesn't work for rates.  Happily, the technique I presented for simple averages can be easily extended to work for rates as well.

To recap a bit, an on-line algorithm for computing a simple average works by accumulating the number of observations and the sum of the observations. There is no mention of time.  For each observation $(x, t)$ we update our state $(s, w)$ according to
\begin{eqnarray*}
s_{n} &=& s_{n-1} + x \\
w_{n} &=& w_{n-1} + 1
\end{eqnarray*}
The exponentially weighted average approach that I mentioned before uses the time of each observation to discount the previous state based on the time of the previous observation.  Thus, the state $(s,w,t)$ is updated according to
\begin{eqnarray*}
\pi &=& e^{-(t-t_{n-1})/\alpha} \\
s_{n} &=& \pi s_{n-1} + x \\
w_{n} &=& \pi w_{n-1} + 1 \\
t_{n} &=& t
\end{eqnarray*}
If we were to compute simple rates without discounting, then instead of interpreting $w$ as a count, we would interpret it as a time interval. Thus we would update state $(s,w,t)$ according to with an observation $(x, t)$ according to
\begin{eqnarray*}
s_{n} &=& s_{n-1} + x \\
w_{n} &=& w_{n-1} + (t - t_n) \\
t_n &=& t
\end{eqnarray*}
By analogy with the discounted version of the simple average, we can produce an exponentially weighted rate average by updating the state according to this
\begin{eqnarray*}
\pi &=& e^{-(t-t_{n-1})/\alpha} \\
s_{n} &=& \pi s_{n-1} + x \\
w_{n} &=& \pi w_{n-1} + (t-t_n) \\
t_{n} &=& t
\end{eqnarray*}
This approach has a defect in that each data pont should really be considered to be a report of events spread out in an uncertain way over the time since the last report. As such, the interval and the events should probably both be discounted as if they had been reported uniformly over the period in question. In practice, this correction does not matter much since the averaging time constant $\alpha$ is typically large compared to the average reporting interval. Another consideration comes up when multiple sources are reporting concurrently. If we want to do proper discounting of each observed number, then we have to keep track of the time since last report for each of the sources. Since this probably won't matter and since it considerably complicates the implementation, I would rather not do it.

See https://issues.apache.org/jira/browse/MAHOUT-634 for code. This will be in Mahout before long.

Update on exponential time-embedded averages

I will be adding this code to Mahout shortly.

See https://issues.apache.org/jira/browse/MAHOUT-634 for status.

Also, if you are measuring rates, then it is common for rates to be reported from multiple sources independently.  Such an average can be computed pretty easily using this same framework if the sources report often relative to the averaging time constant.  This simple implementation just attributes each reported count as if they occurred in the interval since the most recent report from any reporter.  If the time constant is relatively long, this can work out reasonably well as long as we are careful.

If reporting intervals are longer, then the averaging is a bit trickier because we really would like to attribute the reported counts over the entire interval from the last report from the same source.  This means that we have to discount some of the counts because they are effectively kind of old.

More details shortly.

Exponential weighted averages with irregular sampling

Ted Yu on the hbase list wants to compute an exponential weighted average for the number of queries per second or other statistics that characterize the performance or state of an hbase system.

The trick is that the measurements are only available at irregular intervals. If they were sampled regularly, then the standard mixing trick would work:
\[
m_{n+1} = \mu m_n + (1-\mu) x_{n+1}
\]
where $m$ is our current estimate of the mean, $x_n$ is the $n$-th sample and $\mu$ determines how much history to use.

With unequal sample times, things become a bit more complicated. If we get lots of measurements all at once, we want to give them nearly equal weight but if we have a long gap, we want to weight the very old samples much less.

In fact, we want to weight old samples according to how old they are with exponentially decreasing weight. If we sample values $\left \lbrace x_1 \ldots x_n \right \rbrace$ at times $t_1 \ldots t_n$ then we want the weighted mean defined by
\[
m_n = {\sum_{i=1}^n x_i e^{-(t_n - t_i)/\alpha} \over \sum_{i=1}^n e^{-(t_n - t_i)/\alpha} }
\]
Here $\alpha$ plays the same role as $\mu$ did before, but on a different scale. If the evenly sampled data comes at time intervals $\Delta t$ then $\mu = e^{\Delta t / \alpha}$.

Happily, there is a very simple recurrence relationship that allows us to keep only two intermediate values while computing the value of $m_1 \ldots m_n$ in an entirely on-line fashion as the $x_i$ values arrive.

To see this, define

\begin{eqnarray*}
\pi_n &=& e^{-(t_{n+1}-t_n)/\alpha} \\
w_{n+1} &=&
\sum_{i=1}^{n+1} e^{-(t_{n+1} - t_i)/\alpha} =
1+e^{-(t_{n+1}-t_n)/\alpha} \sum_{i=1}^{n} e^{-(t_{n} - t_i)/\alpha} \\
& =& 1 + \pi w_n\\
s_{n+1} &=&
\sum_{i=1}^{n+1} x_i e^{-(t_{n+1} - t_i)/\alpha} =
x_{n+1}+e^{-(t_{n+1}-t_n)/\alpha} \sum_{i=1}^{n} x_i e^{-(t_{n} - t_i)/\alpha} \\
&=& x_{n+1} + \pi_n s_n
\end{eqnarray*}


Then note that
\[
m_{n+1} = {s_{n+1} \over w_{n+1}}
\]

This leads naturally to a procedure that has state consisting of $t, w, m$ which are updated with using new values of $t_n, x_n$ according to
\begin{eqnarray*}
\pi &=& e^{-(t_{n}-t)/\alpha} \\
w &=& 1 + \pi w \\
s &=& x_n + \pi s \\
m &=& {s \over w} \\
t &=& t_{n}
\end{eqnarray*}
Isn't that a kick!

To do this right, however, we need a test. Here are some data vectors computed for $\alpha=5$:


         t          x        pi        w          s          m
1 11.35718  1.5992071 1.0000000 1.000000  1.5992071  1.5992071
2 21.54637 -1.3577032 0.1303100 1.130310 -1.1493105 -1.0168100
3 28.91061 -0.3405638 0.2292718 1.259148 -0.6040683 -0.4797436
4 33.03586  0.7048632 0.4382129 1.551775  0.4401527  0.2836447
5 39.57767  0.3020558 0.2702621 1.419386  0.4210124  0.2966159

Writing a proper unit test is an exercise left to the reader. (but I will be happy to help)

Recent lecture on Mahout for SDForum

I gave a lecture last night on recent additions to Mahout.

The slides are here: http://www.slideshare.net/tdunning/sdforum-11042010

I can amplify this post with answers to any questions that anybody puts in the comments.