Analytics

Showing posts with label math. Show all posts
Showing posts with label math. Show all posts

Sunday, March 2, 2014

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.

Saturday, February 15, 2014

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.

Saturday, April 27, 2013

Blotto Games

I'm currently taking the Model Thinking class on Coursera, which is about applying mathematical models to real-world scenarios in order to derive interesting properties about how people or organizations make choices and why certain things happen. One model that was recently discussed in the class is Blotto games, which has applications in understanding presidential elections, terrorism, trials, hiring, etc. It has a nice mathematical formulation, which is why I decided to learn a bit more about it. The specific Blotto game I want to briefly discuss is as follows: there are two players in the game and they must allocate resources across three different "fronts." They each have an equal number of resources (assumed to be one without loss of generality), and whoever has more resources allocated to a front wins that front. The winner of the game is the one who wins two of the three fronts.

We can think of each player's allocation as a point in Euclidean space lying on the triangle given by $x + y + z = 1$ and $x, y, z \ge 0$. It has vertices at $(1, 0, 0)$, $(0, 1, 0)$, and $(0, 0, 1)$. If we visualize the plane of the triangle we get something like this:


Here, suppose the red dot is some allocation (or strategy). Then we can divide the rest of the triangle into two regions: the green region is the set of strategies which beat the red strategy and the blue region is the set of strategies which the red strategy beats. This follows from the fact that the strategies in the green region are closer to two of the triangle's vertices than the red strategy, and the opposite holds true for the blue strategies. From this picture we can immediately see that the regions close to the vertices are very bad, which is intuitively obvious because you are almost certainly going to lose the two fronts you have not allocated much to. It is also the case that any strategy has a "counter" strategy, so if the two players alternately picked pure strategies there would be no equilibrium.

It turns out, however, that there are mixed strategies for which you are guaranteed to win with probability at least $1/2$ regardless of what the other player chooses, so that allows for an equilibrium. To visualize one such strategy, consider the following picture:
The strategy involves only choosing points inside this hexagon, which has vertices at $(0, 1/3, 2/3)$ and permutations of that. It is not uniform, however, as the strategies near the edge of the hexagon are better than those inside. The exact distribution is constant along the edge of the hexagon, zero at the center, and increasing linearly from the center to the edge. This strategy, along with a number of others, can be proved "optimal" in the sense that your probability of winning is always at least $1/2$. So while none of these strategies is individually optimal, randomizing the choice makes it so that the other player cannot effectively counter your strategy. Since the game is symmetric, the existence of such a strategy solves the game; both players will opt to play one of these mixed strategies, and they will each of a $1/2$ chance of winning.

Mathematical games are often pretty interesting to think about, and even this simple Blotto game has unexpected properties and equilibria. Mixed strategies, in particular, are not very intuitive, yet quite valuable in understanding these games and sometimes even in understanding how the world works.

Sunday, April 14, 2013

Distributed Matrix Computation

Anyone who has worked with MapReduce knows that, while it suits many distributed computations well, it is very difficult to implement matrix computations using the framework. It is possible to represent matrix-based algorithms such as PageRank using MapReduce (e.g. here) but it definitely isn't the most natural or efficient way. Moreover, PageRank is pretty simple as far as matrix computations go, and once you want to do things like Cholesky decompositions for solving systems of linear equations or singular value decompositions it gets very tricky. MapReduce was simply not built for these use cases, so there is nothing to fault it for, but since it is the standard go-to distributed computing framework people have tried to do just about everything with it. Different problems require different solutions, though, and the creators of MadLINQ, a system coming out of Microsoft Research, are trying to overcome the difficulty of implementing complex, distributed matrix algorithms with their new system.

The MadLINQ system shares a lot in common with its cousin DryadLINQ that it integrates with. They both have use the directed acyclic graph (DAG) model of computation where each node represents some unit of computation that depends on the outputs of some other nodes. For example, MapReduce itself can be represented in this model by a set of "map" nodes and a set of "reduce" nodes where an edge exists from each of the map nodes to each of the reduce nodes. The key then lies in how MadLINQ maps a matrix computation into a DAG that can be efficiently executed in parallel. To do this, the authors use the familiar tiling abstraction applied to many parallel matrix operations (e.g. tiled matrix multiplication) where the matrix is broken into disjoint, square sub-matrices to operate on. Each node in the DAG is then an optimized, parallel matrix operation on a set of input tiles from the original matrix or intermediate results. One nice thing is that they abstract away the underlying matrix library for the individual node computations so that improvements in that space (which many people work on) will result in improvements in MadLINQ.

One interesting optimization implemented in MadLINQ is what the authors call "fine-grained pipelining." In both MapReduce and DryadLINQ, a node in the DAG must wait for its predecessors to completely finish processing before reading all of their output and executing. This has two negative consequences: first is the fact that during periods where there is low parallelism in the DAG (you can only break down the units of execution so far before the overhead exceeds the benefits) the cluster will be idle, and second is the fact that network utilization will be more bursty because a node decides to read the entirety of its inputs at once. MadLINQ allows individual nodes to produce partial output and consume partial input, enabling pipelining at the sub-node level (hence the term "fine-grained"). This is particularly useful in the context of matrix computations because many of the individual nodes also execute tile-based algorithms which lend themselves to this granular pipelining. The authors went through quite a bit of trouble to optimize the system as much as possible to really "solve" the use case of distributed matrix computations.

MadLINQ really speaks to how hard problems in distributed systems are. MapReduce has been a huge step forward in the last 10 years, but it has its limitations, and as more and more people and companies are turning to large-scale computations and processing, they will need to invent new solutions to solve the wide array of problems out there. The fundamental design decisions that must be made when building a new distributed system inherently limit its applications in a way that, for example, modern processors are not limited. People (generally) do not design specific processors for their applications anymore, and one would expect distributed systems to eventually reach that point as well, but there is certainly a long road ahead.

Thursday, January 24, 2013

Sampling a Normal

The normal distribution is one of the most commonly used probability distributions because it turns out that many probabilities in the real world are approximately normal as a result of the central limit theorem. The standard normal distribution has a probability density function of

$$f(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x^2}$$
For something so common, its density function is not particularly nice to deal with; its cumulative distribution cannot be expressed in terms of elementary functions. This makes it a bit tricky to sample from as well, since the easiest way to sample a random variable is to sample a value from the uniform distribution on $[0, 1]$ and then invert the cumulative distribution function. Fortunately, people have come up with very clever ways of sampling a normal, and I'll provide a derivation of one, called the Box-Muller transform, here.

Consider the joint distribution of two independent normal variables:

$$f(x, y) = f(x)f(y) = \frac{1}{2\pi} e^{-\frac{1}{2}(x^2+y^2)}$$
This is a distribution across all points in the two-dimensional plane given by their $(x, y)$ coordinates. If we instead consider the polar representation of points, we can transform the given joint distribution via the change of variables formula. Using $x = r\cos{\theta}$ and $y = r\sin{\theta}$ gives (most calculations omitted for brevity)

$$f(r, \theta) = f(x, y) \left| \begin{array}{cc} \cos{\theta} & -r\sin{\theta} \\ \sin{\theta} & r\cos{\theta} \end{array} \right| = f(x, y) \cdot r = \frac{1}{2\pi} r e^{-\frac{1}{2} r^2}$$
Integrating out $\theta$ over the range $[0, 2\pi)$ leaves us with something that we can obtain the cumulative distribution for:

$$f(r) = r e^{-\frac{1}{2} r^2} \Rightarrow F(r) = \int_0^r f(t) dt = \int_0^r t e^{-\frac{1}{2} t^2} dt = 1 - e^{-\frac{1}{2} r^2}$$
Now we can apply the inversion technique to sample $r$; let $U_1$ be a random variable drawn from a uniform distribution on $[0, 1]$. We can compute $r = F^{-1}(1-U_1) = \sqrt{-2\ln{U_1}}$ (we decide to swap $U_1$ with $1-U_1$ to make the result nicer since the two have the same distribution). Notice that we can sample $\theta$ as well because it's actually uniformly distributed on $[0, 2\pi)$, so if $U_2$ is another uniform random variable on $[0, 1]$ then we can take $\theta = 2 \pi U_2$. Putting it all together, we get two (independent) sampled normals

$$x = r\cos{\theta} = \sqrt{-2\ln{U_1}} \cos{(2 \pi U_2)} \\ y = r\sin{\theta} = \sqrt{-2\ln{U_1}} \sin{(2 \pi U_2)}$$
To wrap things up and make it a bit more practical, here's a Java snippet of this in action.

Sunday, January 20, 2013

Erasure Coding

Durability is a critical feature of distributed storage systems. Because these storage systems are typically built on a large number of machines using commodity hardware, individual component failures are expected and common; thus durability comes from how fault tolerance is designed into the system. The original Google File System (GFS) was among the earliest systems built at a scale where these problems started to manifest, and one of the key features was that all data in GFS was redundantly stored in three copies. Combined with the proper distribution of these three copies so they don't fail simultaneously, this design allowed Google to scale out their storage system without losing data. While GFS was quite a landmark piece of technology, we've come a long way since then, and people are constantly devising new ways to achieve durability while optimizing other important features of distributed storage systems like latency, storage overhead, and reconstruction cost.

Storing identical copies of data is a simple strategy for redundancy, but more involved techniques such as erasure coding can greatly improve the storage overhead while providing the same durability guarantees. The idea behind erasure coding is to encode a message of $K$ symbols using $N > K$ symbols such that you can reconstruct the original message using only a subset of the $N$ symbols. For example, you can think of storing three copies of data as erasure coding with $K = 1$ and $N = 3$, where as long as you have one of the three copies you can reconstruct the original data. Another common example is encoding $K$ bits using $N = K+1$ bits where the extra bit stores the parity/XOR of the first $K$ bits; losing any single bit will still allow you to recover the original message. These examples demonstrate how erasure codes are redundant and prevent losing the message even when individual symbols are lost (e.g. a machine storing data dies).

The most common erasure codes are Reed-Solomon codes, which are present in CDs, DVDs, and barcodes to ensure that these media are resistant to partial damage. Reed-Solomon codes extend the idea of the parity bit by leveraging polynomial interpolation. Given the message of $K$ symbols, consider the polynomial of degree (at most) $K-1$ with the symbols as the coefficients. Evaluate the polynomial at $N > K$ points known to both the sender and the receiver, and send the results of those evaluations. As long as $K$ of the resulting symbols are received, the polynomial can be interpolated and thus the original message reconstructed. The details of the algebra on polynomial rings and an efficient decoding algorithm are beyond the scope of this post, but suffice it to say that people have figured out how to use Reed-Solomon codes very effectively. It shouldn't come as a surprise then that the new GFS (discussed here) uses them and the Hadoop File System (HDFS) has a module which adds the functionality.

The story doesn't end there, though. Last summer, a paper came out of Microsoft Research and the Windows Azure Storage (WAS) team that describes a new technique called local reconstruction codes (LRC) used in WAS, Microsoft's equivalent of GFS and HDFS. One of the downsides of Reed-Solomon codes in storage systems is that, if one of your original $K$ symbols is lost, you must read a full $K$ other symbols in order to recover that one symbol. When these symbols are chunks distributed across machines, it can be costly to do that many accesses. LRC can be thought of as an extension of Reed-Solomon codes which adds an additional parameter that determines how many local groups there are within a message. It facilitates the reconstruction of lost symbols by leveraging additional redundancy introduced at the local group level. The work demonstrates that the technique is more efficient both in terms of storage overhead and reconstruction cost than Reed-Solomon codes when both methods are required to provide the same fault tolerance as storing three copies of the data (which is still the standard for systems today).

The theory behind LRC is really cool (and my two-sentence summary doesn't do it justice, so read the paper if you're at all interested), but the additional practical discussions in the paper are equally interesting. Some of the key points that I found particularly enlightening were the following:
  • When retrieving a fragment (unit of a file which is erasure coded) from WAS, a request to reconstruct the fragment is issued simultaneously, so that if the fragment lives on a machine that is under heavy load it's possible that the reconstruction is more efficient and returns first. I imagine that this has a nice effect on the long tail of retrieval times.
  • As files are being written, they are initially stored as three redundant copies, but once a file is "sealed" and can no longer be written to, the system will asynchronously erasure code it and remove the redundant copies. This is so the erasure coding doesn't have any performance impact on the "critical path of client writes."
  • Even the process of erasure coding must be fault tolerant. The WAS system stores progress information for files that are being erasure coded so that fragments which are complete don't need to be encoded again after a failure. Furthermore, they employ many cyclic redundancy checks during the process and perform extensive randomized reconstruction tests prior to committing to the erasure-coded version and deleting the original copies.
Erasure coding is a powerful technique for providing durability in distributed storage systems and seems to be becoming the standard. LRC is a neat innovation in this space, and it will be exciting to see how these technologies evolve as these types of systems become more and more prevalent.

Sunday, January 13, 2013

Fantasy Proleague

Fantasy Proleague is a fantasy league for StarCraft 2 that I've been participating in for the last month or so. I played fantasy baseball for a few seasons in the past, and the concept has always been an interesting optimization problem to me. Recently, the problem of finding the best possible fantasy team for the previous round was posed, and I decided to take a stab at it. I'll briefly describe the rules and then outline the approach I took.

Rules

A fantasy Proleague team consists of six players and one Proleague team. Each player and team you can acquire has an initial cost at the beginning of the round. You begin the season with 30 points to allocate to those six players and one team however you like, with one restriction: you must have one player of each race. Each week, for the four weeks in the round, players and teams accumulate points based on their performance and their costs adjust accordingly. At the end of each week, you are permitted to make two trades, which means taking a player from your team and swapping it for a different player from the pool. A trade is valid so long as the current cost of the player you are trading is at least the current cost of the player you are acquiring and the race requirement is still satisfied. One of the trades may be used on the team as well, with the same cost restriction. You lose one point for every trade you make. The goal is, naturally, to maximize the number of points scored by your team over the course of the round.

Model

Modeling the data is simple; each team has a cost and a point gain for each week (eight total values). The costs represent the cost of the team at the beginning of the week, and the point gains correspond to how many points the team gained in the week that followed. Players are represented in the same way except they additionally have a race.

I'll note that this problem is essentially boils down to search optimization; even without trading it is a knapsack problem, and trading only complicates things further. Luckily, there are only 89 players, eight teams, and four weeks, so a handful of optimizations can be applied to make the search space manageable. The remainder of this post breaks down into two sections: optimizing the initial players/team and optimizing the trades. For most of this discussion, I'm going to ignore the race requirement as it makes the problem much more difficult; fortunately, we still end up with a provably optimal solution (although not in the most elegant way).

Initial Team

Reducing the set of possible initial teams greatly limits the search space, so any pruning we can do before even considering trades is extremely helpful. The key here is identifying which players are strictly superior to which others. We'll define player A to be strictly superior to player B if the following conditions hold (defined similarly for teams):
  • player A's initial cost is at most player B's initial cost;
  • player A's cumulative points are at least player B's cumulative points at the end of each week;
  • player A's cost is at least player B's cost at the end of each week; and
  • there is some week for which player A's cumulative points or cost is strictly greater than player B's cumulative points or cost, respectively.
If all of those conditions hold, it's pretty easy to see that you should never choose player B over player A in your initial team. Whatever you do with player B in the future can be done with player A, and player A has some strict advantage over player B that you may be able to leverage (note that this is not true with the race requirement, because player B may have a race that is necessary for your team). So when choosing initial teams, we enforce that we never choose a player unless all players who are strictly superior to the player in question are on the team already, e.g. if there are six or more other players strictly superior to him, he will never be chosen on the initial team. This relatively simple optimization limits the initial space of a few billion teams to on the order of 100,000.

Trading

To reduce the space of trading, we need a very similar concept to strict superiority. Player A is superior from week w to player B if, from week w onward, the following conditions hold:
  • player A's cumulative points from week w are at least player B's cumulative points from week w at the end of each week; and
  • player A's cost is at least player B's cost at the end of each week.
Then, when deciding whether to trade player A for player B at the beginning of week w, we only consider it if player A's cost that week is at least player B's (required), player B is not already on the team (required), and player A is not superior from week w to player B after subtracting one from B's points (to account for the trade cost). So for each week we can build the list of potential trades for each player after applying this filter.

We can further reduce the search space by the following observation. Suppose we are considering trading player A at the beginning of week w. If players B and C both satisfy the above criteria but player B is superior from week w to player C, there is no reason to explore the space resulting from trading player A for player C since trading for player B is always a better choice. For example, this means that when making the final trade, you only need to consider the player that scores the highest in that last week among all possible players you can trade for, which is a pretty intuitive result.

Result

After applying the optimizations described, I was able to exhaustively search for all potential optimal teams in less than an hour. The best team (for those who actually played) can be found in this forum post. Note that the team there does indeed satisfy the race requirement, which is the last point to discuss. I implemented all of the optimizations ignoring the race requirement because it's much more difficult to determine when a player should always be chosen over another when the race factors in. It turns out, however, that the best team ignoring the race requirement altogether scores the same number of points as the best team I found which satisfies the race requirement. That proves that the team is indeed optimal, although not in the most satisfying way. I've continued to think about a good way of modeling the race requirement without eliminating most of the benefit from the optimizations but have come up empty so far. Regardless, I found interesting theoretical value in the process as well as a practical result; we'll see how it deals with the next round of data once it's in.

Tuesday, December 25, 2012

Almost Surely

The concept of probability when dealing with infinite spaces can be a little confusing and counter-intuitive. One "paradox" that came up while I was watching lectures for a computational finance course on Coursera was the idea that an event which has probability zero of occurring can actually occur. For example, if you choose a random real number between 0 and 1, the event of choosing 0.5 exactly can certainly occur, but any reasonable calculation of the probability will lead to a value of zero. Naturally, I assumed people had formalized this notion, so I did a bit of searching online.

This led me to the concept of events that happen almost surely, which is to say they occur with probability one; alternatively, events that happen almost never have zero probability of occurring. There's nothing too crazy here; an event that happens almost surely simply has an infinitesimally small chance of not happening. So for mathematical purposes it only makes sense to assign the probability as one. To quote the Wikipedia article, "the difference between an event being almost sure and sure is the the same as the subtle difference between something happening with probability one and happening always."

Examples of events that happen almost surely are: flipping a coin an infinite number of times and getting tails at least once, picking a random real number between 0 and 1 and getting an irrational number. Examples of events that happen that happen almost never are: a random variable drawn from a continuous probability distribution taking a particular value, picking a random real number between 0 and 1 and getting 0.5. Moreover, the complement of an almost sure event happens almost never and vice versa. So, to sum things up, an event that happens with probability zero might not never happen, just almost never happen.

Sunday, December 23, 2012

Sieve of Eratosthenes in Scala

The Sieve of Eratosthenes is one of the most efficient ways to find all of the primes numbers below some limit (e.g. 100 million) and an essential tool for any programmer interested in various algorithms contests. Recently, I decided to start doing some Project Euler problems in Scala (my new favorite language) and before long I had to write the sieve.

One tricky thing about working with Scala for those of us coming from Java/C++ backgrounds is that some of the nice constructs can lead to pretty bad performance. More details on Scala performance will be saved for another day, but in this case I was able to come up with a pretty concise and "idiomatic" (quotes because I'm still a Scala newbie) way of implementing the sieve which has quite good performance.



Here, we take advantage of Scala's Stream and Range classes for our implementation of the sieve. The outer Stream foreach does the job of picking every prime that we should check divisibility against, basically combining a while construct (takeWhile) with a conditional (filter). The inner Range foreach is a simple replacement for a for loop, marking the multiples as composite. Note that Streams in Scala are lazily evaluated, which means that the filter actually takes into account the current state of the prime array at each step rather than filtering uselessly upfront.

The foreach construct is one example of a language feature that can sometimes give you unexpectedly bad performance, but the developers of the language/compiler have done a good job to make sure that foreach on a Range is very fast (the overhead on Stream is nontrivial, though, as far as I can tell) which means that our inner loop still performs adequately. All in all, this example runs about 50% slower than a Java implementation using for loops, clocking in at about 1.5 seconds on my laptop, which is not bad at all.

This is probably the first of many posts about Scala since it's a really interesting and expressive language that I also happen to use for my job, so look forward to more interesting discussions on the language!