?

Log in

Friday, January 22nd, 2016

Parallelizing

In 1991 Daniel Bernstein published a short paper on a computational technique for solving diophantine equations of the form p(a) + q(b) = r(c) + s(d): https://cr.yp.to/papers/sortedsums.pdf He describes how a heap can be used to visit all sums p(a) + q(b), with a,b < H, in increasing order, using only near-linear space and near-linear time (in H). You can then set up two heaps and "increment" whichever one has the lowest minimum value, outputting any cases of equality.

Can this technique be efficiently parallelized, on a SIMD machine (such as GPU)?

A simple approach is to divide up the elements in one of the heaps. That is, instead of, say, [math]r(d) = d^4[/math], you could have a family of functions [math]r_i(d) = (192d+i)^4[/math], which partition the search space. Each thread gets a different function to work with. This reduces one of the heaps by a factor of 192, but the other is just as large. Where previously we did H^2 + H^2 heap operations, now we have to do H^2 + H(H/192) heap operations, and the operations on the smaller heap are more efficient in the worst-case as well. Also, every thread needs its own pair of *both* heaps which is somewhat memory-intensive, limiting either the size of the problem that can be tackled or the number of available threads.

We could also partition *all* of the inputs as long as every partition is tested against every other combination of partitions. (This is a natural optimization when some restriction on the equations can be found for some moduli, thus eliminating a broad range of possible solutions.) For example, if we divided all the inputs in half by parity, then we need to check 8 different equations. This again requires a pair of heaps per thread.

We could perhaps continue to use a shared heap for the left-hand side of the equation, at the cost of threads becoming idle when their minimum value grows too large. I don't know how to evaluate what this does to the speedup. It also introduces some serialization as the heap operation itself can't be parallelized.

A variation on the shared heap would "unroll" the left-hand heap and have a shared array of the next K values. I can't decide whether this is helpful or not.

On a side note, perhaps the availability of cheap large-capacity flash devices may make it feasible to apply techniques previously dismissed as infeasible due to space requirements.
(Leave a comment)

Friday, January 1st, 2016

How I've been entertaining myself

I participated in Advent of Code, and Advent calendar of programming puzzles. Of the 25 days of puzzles, one I solved without doing any coding, one other by minor scripting of an existing application, and the rest in Python. It was a fun exercise. (Sadly, I did not submit a solution quickly enough on the last day to make the leaderboard.)

For Christmas I received a NVIDIA Jetson TK1, a development board for the Tegra K1 "system-on-a-chip". Most notably, it comes with a 192-core GPU built in, making it a 5"x5" supercomputer.

I thought a good project to get more familiar with programming it would be to reimplement some of the Advent of Code problems using CUDA (Nvidia's programming API for GPU computing.) I posted my first effort on github: https://github.com/mgritter/advent-of-cuda/tree/master/day4 (But, I haven't rerun my old Python solution to demonstrate the speedup!) The new code is capable of finding the solutions to the original problems in less than a second each, and of finding solutions to the "next-hardest" version of the problem in about 9 seconds.

I also have a bunch of new games to play: Mini Metro (a fun logistics game), Concrete Jungle (deck-building + city-building), Infinifactory (factory layout for evil aliens), and Offworld Trading Company (a logistic RTS.) Yes, even my games are really all about programming. :)
(Leave a comment)

Monday, November 2nd, 2015

Bones was right about the transporter

Star Trek canon records at least 11 transporter accidents:
* Splitting Kirk into good and evil halves (TOS: The Enemy Within)
* Transport to the mirror universe (TOS: Mirror, Mirror) --- not counting later deliberate action
* Transport to another universe (TOS: The Tholian Web)
* Two deaths in Star Trek: The Motion Picture
* Weyoun 5 dies in DS9: Treachery, Faith, and the Great River
* Clone of Riker (TNG: Second Chances)
* Reverting four individuals to younger versions of themselves (TNG: Rascals)
* Travel through time (DS9, Past Tense)
* Merging multiple individuals (VOY, Tuvix)
* Damaging the holographic doctor's mobile emitter (VOY, Drone)
* Rendering two crew members incorporeal (TNG, The Next Phase)

To calculate a malfunction rate, we need some estimate of how many times the transporter was used in total during this period. There are 79 TOS episodes, 178 TNG episodes, 176 DS9 episodes, and 172 VOY expisodes. That gives 605 episodes. If we generously estimate 4 transports per episode that gives a failure rate of 0.4%, unacceptably high.

Transporters appear to be able to cross planetary distances easily (no waiting for proper orbital alignment.) An estimate of 25,000 miles per transport --- that is, assuming starships are usually in high orbit --- and 6 people at a time gives an accident rate of around 30 per billion passenger-miles. Modern cars and light trucks have a fatality risk of 7.3 per billion passenger-miles. (However, it should be noted that many of the transporter accidents were reversible or otherwise non-fatal.)
(2 comments | Leave a comment)

Monday, July 20th, 2015

Fun with Bad Methodology

Here's a Bill Gross talk at TED about why startups fail. Mr. Gross identifies five factors that might influence a startup's success, analyzes 250 startups according to those factors, and analyzes the results.

This is possibly the worst methodology I have seen in a TED talk. Let's look at just the data presented in his slides. I could not find a fuller source for Mr. Gross's claims. (In fact, in his DLD conference talk he admits that he picked just 20 examples to work with.)

I did a multivariate linear regression--- this does not appear to be the analysis he performed. This produces a coefficient for each of the factors:

Idea	 0.021841942
Team	 0.033134009
Plan	 0.047012997
Funding	-0.03324002
Timing	 0.133669929


While this analysis agrees that "Timing" is the most important, it differs from Mr. Gross on what is second. It actually says that the business plan is a better predictor of success. That's strike one--- the same data admits multiple interpretations about importance. Note also that linear regression says more funding is actively harmful.

The second methodological fault is taking these numbers at face value in the first place. One might ask: how much of the variation in the dependent variable do these factors explain? For linear regression, the adjusted R^2 value is about 0.74. That means about 25% of the success is explained by some factor not listed here. What checks did Mr. Gross do to validate that his identified factors were the only ones?

While we're still on problems with statistical analysis, consider that the coefficients also come with error bars.

	Lower 95.0%	Upper 95.0%
Idea	-0.082270273	0.125954157
Team	-0.082753104	0.149021121
Plan	-0.049741296	0.143767289
Funding	-0.126226643	0.059746603
Timing	 0.034672747	0.232667111


The error bars are so wide that few useful statements can be made about relative ordering of importance. (Obviously this might change with 230 more data points. But it seems like Mr. Gross was operating on only 20 samples.)

Next let's transition to the data-gathering faults. Mr. Gross's sample of companies is obviously nonrandom. But it is nonrandom in a way that is particularly prone to bias. He lists startup companies that actually got big enough to attract sufficient attention that he'd heard of them! The even mix of success and failure when he's looking to explain the high rate of failure should be a huge red flag.

Suppose we add 15 more failed companies to the list that have an average timing of 7 (slightly better than the sample) but average on the other factors of 5.

Oops, now timing has slipped to third!

Idea	 0.079811441
Team	 0.067029361
Plan	-0.007925792
Funding	 0.033260711
Timing	 0.064723176


Not surprising, because I deliberately set things up that way. But Mr. Gross doesn't know what the "next" set of failures look like either. He has a complete list of IdeaLabs companies, but not the outside world, which is its own bias--- maybe it's Mr. Gross who is prone to timing failures, not the rest of the world!

Picking only large, well-known failures for your analysis is nearly the very definition of survivorship bias.

Finally, the inputs themselves are suspect even where they are complete. Mr. Gross already *knows* whether these companies are successes or failures when he's filling out his numbers. Does Pets.com really deserve a "3" in the idea category, or is that just retrospective thinking? Why are the successful IdeaLabs companies given 6's and 7's for Team, but the unsuccessful ones get two fours and a five? Did Mr. Gross really think he was giving those two ventures the "B" team at the time?

Even the Y axis is suspect! Pets.com had a successful IPO. Uber and AirBnB can't claim to have reached that stage--- they might yet implode. (Unlikely, admittedly, but possible. And their revenue numbers are better than Pets.com ever was.) As an investor, "Did I receive any return on my investment" is the measure of success.

To summarize,
* The data were generated from the opinions of the principal investigator and not subject to any cross-checks from other parties.
* The data exhibit survivorship bias and other selection bias.
* The analysis includes no confidence intervals or other measurements of the meaningfulness of the results.
* The results presented appear highly dependent upon the exact analysis performed, which was not disclosed.

Am I being unfair? Perhaps. This was not an academic conference talk. But if we are to take his advice seriously, then Mr. Gross should show that his analysis was serious too.
(Leave a comment)

Saturday, July 4th, 2015

Stupid NFS Tricks

One of the Tintri SE's asked whether the advice here about tuning NFS was relevant to us: http://docstore.mik.ua/orelly/networking_2ndEd/nfs/ch18_01.htm (It is not.)

Back in the days where CPU cycles and network bits were considered precious, NFS (the Network File System) was often run over UDP. UDP is a connectionless, best-effort "datagram" service. So occasionally packets get lost along the way, or corrupted and dropped (if you didn't disable checksumming due to aforesaid CPU costs.) But UDP doesn't provide any way to tell that this happened. Thus, the next layer up (SUN-RPC) put a transaction ID in each outgoing call, and verified that the corresponding response had arrived. If you don't see a response within a given time window, you assume that the packet got lost and retry.

This causes problems of two sorts. First, NFS itself must be robust against retransmitted packets, because sometimes it was the response that got lost, or you just didn't wait long enough for the file server to do its thing. This is not so easy for operations like "create a file" or "delete a file" where a second attempt would normally result in an error (respectively "a file with that name already exists" and "there's no such file.") So NFS servers started using what's called a "duplicate request cache" which tracked recently-seen operations and if an incoming request matched, the same response was echoed.

The second problem is figuring out what the appropriate timeout is. You want to keep the average cost of operations low, but not spend a lot of resources retransmitting packets. The latter could be expensive even if you don't have a crappy network--- say the server is slow because it's overloaded. You don't want to bombard it with extra traffic.

Say you're operating at 0.1% packet loss. A read (to disk) normally takes about 10ms when the system is not under load. If you set your timeout to 100ms, then the average read takes about 0.999 * 10ms + 0.000999 * 110ms + 0.000000999 * 210ms + so on, about 10.1ms. But if the timeout is a second, that becomes 11ms, and if the timeout's 10 seconds then we're talking 20ms.

So, at least in theory, this is worth tuning because a 2x difference between a bad setting and a good setting is worthwhile. Except that the whole setup is completely nuts.

In order to make a reasonable decision, the system administrator needs statistics on how long various NFS calls tend to make, and the client captures this information. But once you've done that, why does the system administrator need to get involved? Why shouldn't the NFS client automatically adapt to the observed latency, and dynamically calculate a timeout value in keeping with a higher-level policy? (For a concrete example, the timeout could be set to the 99th percentile of observed response times, for an expected 1% over-retransmit rate.) Why on earth is it better to provide a tuning guide rather than develop a protocol implementation which doesn't need tuning? This fix wouldn't require any agreement from the server, you could just do it right on your own!

Fortunately, in the modern world NFS mainly runs over TCP, which has mechanisms which can usually tell more quickly that a request or response has gone missing. (Not always, and some of its specified timeouts suffer the same problem of being hard-coded rather than adaptive.) This doesn't remove the first problem (you still need the duplicate request cache) but makes an entire class of tuning pretty much obsolete.

Nothing upsets me more in systems research and development, than a parameter that the user must set in order for the system to work correctly. If the correct value can be obtained observationally or experimentally, the system should do so. If the parameter must be set based on intuition, "workload type", expensive consultants, or the researcher's best guess then the real problem hasn't been solved yet.
(Leave a comment)

Friday, May 29th, 2015

Internet arguments about math I have participated in, because I'm on vacation

1. Is there such thing as an uncomputable integer? Or, even if a function F(x) is uncomputable, does, say F(10^100) still count as computable because, hey, some Turing machine out there somewhere could output it! Or does a number count as computable only if we (humans) can distinguish that machine from some other machine that also outputs a big number?

Consider that if F() is uncomputable, it seems extremely dubious that humans have some mechanism that lets them create Turing machines that produce arbitrarily-large values of F.

I'm probably wrong on the strict technical merits but those technical grounds are unduly Platonist for my taste.

2. Is trial division the Sieve of Eratosthenes? Maybe it's the same algorithm, just a particularly bad implementation!

I'm with R. A. Fisher on this one. If the idea was just some way to laboriously generate primes, it wouldn't be famous, even among Eratosthenes' contemporaries. The whole thing that makes the sieve notable is that you can apply arithmetic sequences (which translate into straight lines on a grid!) to cut down your work.

The counterargument seems to be that Eratosthenes might well have performed "skip, skip, skip, skip, mark" instead of jumping ahead by five.
(Leave a comment)

Saturday, May 2nd, 2015

Never Trust Any Published Algorithm

The slides from my Minnebar talk, "Never Trust Any Published Algorithm", can be found here.

I got a couple good stories in response.

A local person who works in computer vision talked about finding a new algorithm for shadow detection. It worked great on the graduate student's ten examples. It didn't work at all on his company's hundreds of examples in all different lighting and weather conditions.

A coworker shared:

A few years ago my brother was working at an aircraft instrument company. He was
looking for an algorithm that described max loading for a plane still able to take off
based on many variables.

He found what was supposed to be the definitive solution written by a grad student at
USC. He implemented the algorithm and started testing it against a large DB of data
from actual airframe tests. He quickly found that the algorithm worked just fine for
sea-level examples, but not as airport altitude rose. He looked through the algorithm
and finally found where the math was wrong. He fixed his code to match his new
algorithm and found that it now matched the data he had for testing.

He sent the updated algorithm to the grad student so he could update it on public
servers. He never heard back nor did he ever see the public code updated.


The example I put as a bonus slide wasn't previously mentioned in my series of blog posts is a good one too. In Clustering of Time Series Subsequences is Meaningless the authors showed that an entire big data technique that had been used dozens of times in published papers produced essentially random results. (The technique was to perform cluster analysis over sliding window slices of time series.)
(Leave a comment)

Friday, April 24th, 2015

Infinitesimals in the family of asymptotics

I answered this question about the Master Theorem (which is used to calculate the asymptotic behavior of divide-and-conquer algorithms): Why is the recurrence: T (n) = sqrt(2) *T (n/2) + log(n) a case 1 of master method?

The source of confusion (see the comments to my answer) seemed to be that the original poster did not really understand the asymptotic behavior of log(n) and n/log(n). In fact, he went so far as to post a graph "proving" that n^0.99 is larger than n/log(n). However, this is false for large enough numbers (large enough being somewhere around 10^100 in this case.) The logarithm grows more slowly than any positive power of n. As a consequence, n/log(n) grows asymptotically faster than any positive power of n less than 1.

What I realized is that this might be some student's first (and maybe only!) experience with infinitesimals! Normally we think of infinitesimals as numbers "smaller than any positive real number". For example, before epsilon-delta proofs and modern analysis, the differential and integral calculus informally treated dx as an infinitesimal value. But while students are told the "idea" is that dx is a very small value, they are repeatedly cautioned not to treat it as an actual number.

The surreal numbers (and thus the class of combinatoric Games too) have infinitesimals, but they are far outside the normal curriculum.

So what model is a student supposed to apply when confronted with O(log n) or Θ(log n)? It behaves exactly like an infinitesimal in the family of asymptotic limits. big-O = bounded above, Ω = bounded below, and Θ = bounded both above and below, thus:

log n = Ω( 1 )

log n = O( n^c ) for any c > 0, i.e., log n = O(n), log n = O(sqrt(n)), log n = O(n^0.0001)

n^k log n = O( n^c ) for any c > k.

n / log n = Ω( n^c ) for any c < 1

Nor does taking a power of the logarithm make any difference.

log^100 n = (log n)^100 = O( n^c ) for any c > 0

Crazy, right? But we can get Wolfram Alpha to calculate the crossover point, say, when does (log n)^100 = n^0.01? At about 4.3*10^50669.

That's what makes logarithm an infinitesimal. No matter what power we raise it to (think "multiply") it is still smaller in its asymptotic behavior than the smallest power of n (think "any positive real number.") And there's a whole family of infinitesimals hiding here. log log n is infinitesimally smaller than log n. Don't even ask about the inverse Ackermann function or the iterated-log function.

So it's not surprising that students might have difficulty if this is their first journey outside the real numbers. Everybody can handle the fact that O(n^2) and O(n^0.99) and O(n^0.5) are different, but the fact that none of these examples will be O(log n) is kind of perplexing, because O(log n) is obviously bigger than O(n^0). (Jumping to exponentials like O(2^n) seems like an easier stretch.)

What was your first encounter with an infinitesimal?
(Leave a comment)

Wednesday, March 18th, 2015

Sums of Cubes

Somebody on Quora asked for an example of w^3 + x^3 + y^3 = z^3 in positive integers, all relatively prime. Turns out they were looking for a counterexample of a theorem they thought they proved, which is sort of a passive-aggressive way to approach it.

Anyway, brute force solves the day. Noam Elkies came up with this monster set of equations to characterize solutions to the "Fermat Cubic Surface": http://www.math.harvard.edu/~elkies/4cubes.html, reducing the search space dramatically since we only need to look at triples (of both positive and negative integers) and filter the resulting cubic equations to those which have only 1 negative value and the relative-primeness condition. Here's Elkies' characterization in python:

def elkies( s, r, t ):
    w = -(s+r)*t*t + (s*s+2*r*r)*t - s*s*s + r*s*s - 2*r*r*s - r*r*r
    x = t*t*t - (s+r)*t*t + (s*s+2*r*r)*t + r*s*s - 2*r*r*s + r*r*r
    y = -t*t*t + (s+r)*t*t - (s*s+2*r*r)*t + 2*r*s*s - r*r*s + 2*r*r*r
    z = (s-2*r)*t*t + (r*r-s*s)*t + s*s*s - r*s*s + 2*r*r*s - 2*r*r*r
    return (w, x, y, z)


The smallest one I found was:

365^3 = 107^3 + 209^3 + 337^3


and the first with two different solutions:

67477^3 = 36919^3 + 46049^3 + 54205^3
          26627^3 + 46747^3 + 57103^3


and just for posterity, here are the first 1000:
Read more...Collapse )
(Leave a comment)

Saturday, February 28th, 2015

More on not trusting published algorithms

If you haven't read this paper on how not to implement the Sieve of Eratosthenes in a functional language, it's well worth your time: "The Genuine Sieve of Eratosthenes" by Melissa O'Neill of Harvey Mudd. It's about how a common simple functional programming example is completely wrong. (Her talk at Stanford on a new family of random number generators is worth the time too.. It's a cool hack to incorporate permutations into linear-congruential RNGs.)

I ran into a question on Quora about summing the first K prime numbers, and the person was asking whether there was a better way than trial division. When I pointed out that the Sieve was much better (even at relatively small scales!) his response was "then I need to fill all integers upto INTEGER.MAX_VALUE"--- i.e., how do we bridge the gap between "All the prime numbers up to N" and "The first K prime numbers".

There are three ways to tackle this and they're all interesting in their own way.

1) The Sieve is so much better than trial division that it almost doesn't matter. The prime number theorem tells us that there are about N/log(N) prime numbers less than N, and that the Kth prime is about K log K. So, trial division to find the Kth prime require work that's around O( (K log K)*K ) = O( K^2 log K ). We might be generous and note that most composite numbers don't require testing all previous primes in order to find a divisor, and bound below by Omega(K^2).

The sieve takes O(N log log N). If you're looking for, say, the 1000th prime, you can easily sieve up to 1000^2 = 1,000,000 and still come out ahead. (Particularly since, in practical terms, the Sieve doesn't require any division.)

2) The prime number theorem lets us come up with a pretty-good bound on how large an N we need to find the K'th prime. Run the sieve on, say, N = 2 K log K and that's a solid win.

3) Even if we could not predict where the Kth prime was likely to occur (or, say some subset of the primes of special interest), the Sieve can still be run to an indeterminate stopping point at only a constant-factor slowdown.

The paper linked above gives a sophisticated way of doing this, but brute force is fine too.

Say we sieve up to N=1000 and don't yet find the primes we're looking for. Then double N and start over from scratch. Do this until we find K primes; iteration J of the sieve has size N =(2^J)*1000. If we had an oracle to tell us J, we could have done one sieve pass. Instead we did sieves of size 1000, 2000, 4000, 8000, ... , 2^J*1000, but the sum of the sizes is just 2^(J+1)*1000 since it's a geometric series. So even if our answer was at 2^(J-1)*1000 + 1, the total amount of numbers sieved was at most 4x the minimum we needed to do.

And that's not the best we can do, if we re-use the first 1000 numbers rather than throwing them away and starting from scratch it gets somewhat better--- but that's just quibbling about constant factors. Even treating the Sieve as a black box gives us a better result than trial division.
(Leave a comment)
Previous 10 | Next 10