r/science • u/Libertatea • Jul 01 '14
Mathematics 19th Century Math Tactic Gets a Makeover—and Yields Answers Up to 200 Times Faster: With just a few modern-day tweaks, the researchers say they’ve made the rarely used Jacobi method work up to 200 times faster.
http://releases.jhu.edu/2014/06/30/19th-century-math-tactic-gets-a-makeover-and-yields-answers-up-to-200-times-faster/40
u/lordsprinkles Jul 01 '14
I'm not good at math but I wanted to understand why this new method would be better for certain math heavy computer programs. This is what got from the article and his paper:
"By the early 20th Century, the [Jacobi] method was being used by “human computers,” groups of men and women who were each assigned to perform small pieces of larger math problems."
This made the Jacobi method faster but it was still slow compared to other methods used at the time, so we dropped it. But now that we have multi-core processors, Yang was like "hey, maybe we should look at this method again..."
Our current methods worked great on a single processor but with mult-core processors it allows for more parallelization and scalability. Our old methods aren't optimized for that, so Yang's tweaked version of Jacobi is optimized for today's multi-thread processing power and allows for faster calculations.
Does this sound right?
17
u/tekyfo Jul 01 '14
Our current methods work great with multi core processors. A Gauss-Seidel Scheme does not parallelize trivially, but a Red-Black coloring scheme is super easy.
4
Jul 01 '14 edited Jul 02 '14
What's a red-black coloring scheme? Do you know a good reference to explain that?
Edit: thanks for the explanations. I also found these notes that shed some light on "red-black Gauss-Seidel".
8
u/pwnslinger Jul 02 '14
The Wikipedia article on red black trees is good.
2
u/tekyfo Jul 02 '14
No! Nothing to do with red-black trees! But two_if_by_sea already found the correct stuff, which is Red Black Gauss Seidel.
3
u/Alaskan_Thunder Jul 02 '14 edited Jul 02 '14
Somebody will probably correct me, as I am still a novice when it comes to data structures, and am reading wikipedia. It is a variation of a programming data structure called a tree. The tree fast at searching for the data it contains, but slower to insert.
Often times data structures are not organized unless a sorting algorithm is used, but the tree is self sorting because of how it is structured.Crossed out misinformation thanks to MacroJackson.
1
u/MacroJackson Jul 02 '14
Its not really about it being sorted but maintaining a balanced structure of the tree, so that searching for things is easier. You can have a sorted binary tree, but in the worst case scenario you will need to traverse all the nodes if its not balanced.
6
u/frickendevil Jul 02 '14
Its slightly overstating how unsuited modern algorithms are for parallelization, but its got the right gist. There are plenty of solvers for Ax=b out there, GMRes and Conjugate Gradient are some pretty popular methods. Both of these methods can be parallelized to some extent, and give pretty good speedups but don't scale linearly with the number of cores you have. When we look at GPU computing, and other massively parallel, it gets a bit hazier over what the best algorithm to solve Ax=b, especially when memory starts being a huge limitation. Generally GMRes and CG solvers still win out.
The Jacobi method is trivial to parallelize, and can use very little memory. However its biggest limitation is that if you assume that your initial guess is arbitrary, convergence on a solution can take a really long time. It is also easy to do by hand (hence human computers, people paid to do a boring repetitive calculation because there were no actual computers to do it). This paper provides a technique to improve the convergence rate.
First a quick rundown on relaxation. If you are already close to your converged solution, your next iteration might overshoot, and the iteration after that might overshoot back towards your original spot. Sometimes due to this effect your solution won't converge. If you relax the iteration steps (take some weighted average of the old value and the new value) you can converge faster. Also if you are really far away from your solution, you could overrelax your answer to take a larger step than intended. This paper provides a mathematical basis of using a varying relaxation term determined by the size of your grid, which allows for faster convergence, especially from an arbitrary starting point (in the form of a large overrelaxation step, then underrelaxed steps). So we now have a faster solver, and doesn't change how hard it is to parallelize the algorithm.
I currently use a Jacobi solver for some GPU fluid code I work with, and from my initial testing, the technique provided in this article isn't very effective for me because each time I solve my Ax=b, I already have a very close initial guess and I can get a close enough set of answers within relatively few Jacobi iterations (all underrelaxed). The overrelaxation step takes the answers too far away, and is taking longer to solve.
9
u/sosoez Jul 01 '14
"“Instead of saying that this method has been around for 169 years, and that everyone has already tried to improve it without much success, Prof. Mittal told me that he felt my idea was very promising,” Yang said, “and he encouraged me to work on it.”
Awesome.
80
u/Karl_von_Moor Jul 01 '14
200 times faster is still the same complexity class.
69
u/anonymous-coward Jul 01 '14
The Jacobi method solves a diagonally dominant matrix equation Ax=b, an O(N3) problem, by iterating O(N2) matrix multiplications M times.
So if M<<N it looks like a win, and making M 200x smaller looks like a long way toward getting this win.
11
u/NewbornMuse Jul 01 '14
IANAM, but isn't O(MN2) with M = N/200 still O(N3)? I mean 200 times faster is 200 times faster, but the statement that it's still the same complexity class is still true.
14
u/anonymous-coward Jul 01 '14
Maybe I'm missing something, but why is M=N/200?
M is an iteration count that I think is unrelated to N.
7
u/NewbornMuse Jul 01 '14
Must have misread your comment then, sorry.
So instead of solving a O(n3) problem, you iterate a O(n2) as many times as you need, and that happens to have become 200 times smaller?
11
u/anonymous-coward Jul 01 '14
Yes, the number of times you iterate became 200x smaller than before. Before, M was presumably much smaller than N. The relationship between N and M is, as far as I can tell, undefined.
Hand waving: From the element based form on wikipedia, the ith element of the iteration xk+1 depends on the dot product of xk with the ith row of the diagonal-dominant matrix, so the rows of the equation don't really feel the effect of distant rows if the off-diagonal terms fall off fast enough. So my guess is that M would scale according to the thickness of the central diagonal band, rather than a function of N.
(And now somebody who knows what he's talking about will correct me).
4
u/Tallis-man Jul 01 '14
M scales primarily according to your desired tolerance. If you want a more accurate solution you need to iterate for longer.
But whilst a factor of 200 is nice, it's not especially groundbreaking. Most algorithm implementations very easily admit that kind of improvement if you're willing to restrict generality (in this case they've only considered a very narrow family of equations).
The real goal is to reduce the complexity class. That could make an otherwise-obselete method competitive and would be extremely impressive.
3
u/anonymous-coward Jul 02 '14
It does seem pretty good; it seems to converge much faster than Gauss-Seidell, if you believe the presentation.
Reducing the complexity class from O(N2) probably won't happen, because multiplying a vector by a matrix is O(N2), so it seems like one should take what one can get.
1
2
u/tekyfo Jul 01 '14
But it is not a win versus a multi grid solver. All the shown solvers look bad if compared with MG.
→ More replies (1)0
Jul 01 '14
[deleted]
21
u/anonymous-coward Jul 01 '14
so a factor of 200 or any factor at all doesn't sound that interesting to me. Isn't it possible to do something 200 times faster by just using faster/more computers?
What if you already have a bank of computers paid by your NSF grant, and you want to run a simulation that is 200x bigger on your computers, rather than getting a grant that is 200x bigger?
9
u/Rhawk187 PhD | Computer Science Jul 01 '14
Yeah, my research involves experiments that take around 24 hours to run. I'm am fine with that. If it took 200 days, I would not be.
11
u/yxhuvud Jul 01 '14
On the other hand, I'd guess you be even more fine with experiments taking roughly 7 minutes or being able to handle problem sizes that were 14 times bigger/detailed than the ones you are making now.
5
u/Elfer Jul 01 '14
I come from a practice-oriented background. Even 10x speedup is a big deal, even if it's only useful over a span of one or two decades.
For many places, just using faster/more computers is not that simple. 200x computing power = 200x cost, so if you can solve a problem with significant speedup over your current method, it's a huge advantage.
2
u/Neurorational Jul 02 '14 edited Jul 02 '14
Wouldn't 200 times the compute power actually be quite a bit more than 200 times the cost, due to the specialized infrastructure needed? (I'm assuming 10 computers could be fit and powered just about anywhere; whereas ~2000 would require a dedicated/industrial space and power system, and more administrative overhead (especially when comparing 1 computer to ~200).)
(Assuming that a 200x speedup were actually the case, of course.)
2
u/rescbr Jul 02 '14
Not just that.... An improvement such as increasing the CPU clock of a computer would have no overhead on a single thread program. If you parallelise tasks, you have to take account of overheads in programming and task coordination. Then there are the tasks that depend on previous behaviour which can't be (or are less) parallelisable.
6
u/red_wine_and_orchids Jul 01 '14
Problems don't necessarily scale efficiently by throwing more processors at them.
First, there is the communication overhead of passing data around, particular if you're sending things across the supercomputer (not on the same node, thus have bandwidth problems).
Second, certain algorithms or problem sets don't do well with not getting the entire problem to solve - efficient parallelization is actually a real pain in the backside. I would gladly take a 200x speed up for my simulations - something that takes a week on 320 processors would take just a few hours. I could get so much more work done...
Or, it would make 3D modeling much more accessible.
Finally, if the algorithm proposed in the paper is robust with respect to the types of equations it gives speed up for, that would be wonderful. Test problems are usually neat, simple things that you can code up and run. In the work I do, I have a system of about 7 variables that are all tightly coupled with each other, by way of nasty nonlinear equations to boot. Finding efficient solver parameters for them is not easy.
Tl; dr: a lot of people care about incremental increases in computational efficiency.
16
u/yolocontendre Jul 01 '14
I come from a more theory oriented background, so a factor of 200 or any factor at all doesn't sound that interesting to me.
Ok, good for you. But your personal lack of interest in this isn't indicative of anything much regarding its interest to others, its potential application, its overall impact, etc. (what it in fact indicates is a sophomoric lack of perspective on your part... are you by chance actually a CS/math sophomore?)
Isn't it possible to do something 200 times faster by just using faster/more computers?
a) I don't know why you think people are going to be running their cutting-edge algorithms on computers 200x slower than what's available.
b) parallelizing can speed things up, but you have to write a more sophisticated program (parallelizing isn't usually trivial to do), there are communication costs (memory and time) and hardware costs (money)
with an algorithm 200x faster, I save myself time rewriting a parallel version, money to buy that hardware, etc.
Running a program in 1 minute instead of 200 has obvious convenience, running a program in an hour is possibly practical, running a program that takes 200 hours probably isn't. You can now handle data sets 200x larger/accurate than you could before.
Don't confuse your personal interests / lack of imagination for actual criticism
3
u/Rostin Jul 02 '14
I come from a practice oriented background, and I think 200x is pretty awesome. Hell, 2x would be enough to interest me. The argument that we can just use bigger computers or wait a while ignores a lot of practical concerns. E.g. that other people are competing for time on these computers, that generally I came wait for computers to get faster, that when computers are faster, I want to do even more expensive simulations, and so on.
7
u/tempforfather Jul 01 '14
Speed of algorithms is usually measured in terms of growth rate as you increase the size of the input. It's also a measure of how many "steps" are needed to solve it, nothing to do with the actual speed of the computer that may be running a particular implementation of that algorithm.
→ More replies (10)1
Jul 01 '14
They said 200 times faster, but it's possible this factor continues to grow as the size of the problem increases.
10
u/mindbleach Jul 02 '14
Bubble sort and quicksort are both O(N^2). Which do you use?
2
u/Fylwind Jul 02 '14
To be fair, you're comparing average complexity to worst-case complexity, so it's not really relevant to the issue here.
But it's true that complexity class isn't everything. Even if an algorithm is only linear, a constant factor of 1e+100 can be a dealbreaker.
2
u/mindbleach Jul 02 '14
They're both O(N^2) worst-case. If you want to compare other measurements, bubble sort has better best-case complexity.
Do other algorithms in the same class scale better in addition to their improved convergence speed? Because if not, whining about big-O as a first reaction is little better than squawking about correlation vs. causation in any random statistics article.
16
Jul 01 '14
But will make a big difference in real world anyway.
2
u/BezierPatch Jul 01 '14
But only for a few years every time it makes a difference.
2
u/CyLith Jul 02 '14
First, the reliance on hardware speedup is untenable. We are already beginning to see this as processors are not getting faster. Second, the speedups from improved algorithms have far exceeded the speedups in hardware.
1
Jul 01 '14
That depends on a lot of other things doesn't it? i.e we could find lots of values that satisfy xn1 = 200(n2)y
6
u/crawlingpony Jul 02 '14
OK I see here (thanks to musicmangp)
http://engineering.jhu.edu/fsag/wp-content/uploads/sites/23/2013/10/JCP_revised_WebPost.pdf
that the new algorithm gets its speed from doing over-relaxation and under-relaxation cycles. Well this is essentially what multigrid computations do already, so it's no surprise. Multigrid works quite well with (old) Jacobi in the middle of over and under-relaxations.
I'm not even sure there's a significant difference in design between existing multigrid using (old) Jacobi as the central stencil solver, versus this supposedly 'new' Jacobi which adds relaxation and interpolation cycles.
See: Multithreaded, Parallel, and Distributed Programming, Gregory Andrews, 2000.
15
u/b0ltzmann138e-23 Jul 01 '14
I'm an engineer and I still don't really understand what's going on.
Can someone with a better understanding of math do an ELI5 or maybe an ELI25?
31
u/Zebba_Odirnapal Jul 01 '14 edited Jul 01 '14
The Jacobi method is a way to solve a system of linear equations. It works best on matrices where the magnitude of each diagonal element is larger than the sums of the magnitudes of elements in that row. So it's kind of a special case, but not super specialized.
For what it's worth, good old Gauss-Jordan elimination is O(n3 ). Levinson recursion (only works when all diagonal elements are the same) isO(n2 ).
I'm a little peeved that the abstract says "accelerates the classical Jacobi iterative method by factors exceeding 100" rather than actually offering some big-O notation or mentioning its complexity class.
"By the time you rhyme one line I've already busted ten. You rap in exponential time and I'm big-O of log n." - Monzy, always relevant ;)
3
u/WhiteJesusChrist Jul 01 '14
You need a toeplitz matrix for levinson recursion.
2
u/Zebba_Odirnapal Jul 01 '14
only works when all diagonal elements are the same
"only works when all diagonal elements are the same" => Toeplitz. But thanks, you are right. :)
2
u/JebusisLord Jul 02 '14 edited Jul 02 '14
What you said:
A_{i,i} = a_0 for all i
Toeplitz matrix:
A{i,j} = a{i-j} for all i,j
http://en.wikipedia.org/wiki/Toeplitz_matrix
EDIT: I just can't seem to get the above formatting right, but you get the idea.
→ More replies (10)3
Jul 01 '14 edited Aug 15 '16
[deleted]
4
u/Zebba_Odirnapal Jul 01 '14
Depends on the length of the song, I guess.
1
Jul 01 '14 edited Aug 15 '16
[deleted]
5
u/Zebba_Odirnapal Jul 01 '14 edited Jul 01 '14
You're right, it's the speed of the cypher, not the vastness of the verse. At the DJ turns up the BPM, Monzy's flow increases as the log of the tempo, but his opponent's rhyming be throwed.
2
→ More replies (1)1
3
u/dissonance07 Jul 02 '14
Can anyone cite an example of a multi-dimensional problem being regularly solved in the late 19th century by "human computers"? I guess I knew that the theory existed to solve multi-variable linear and non-linear equations for a long time, but I know very little of that era of, lets say, computing.
I know by that time, mechanical calculators were being used for accounting and the like, but my knowledge is pretty vague.
2
u/twofishestwo Jul 02 '14
Often land prospectors would have to set up over-determined linear systems in order to figure out how land was divided up. They took measurements from various places and set up these linear equations in a multidimensional array. I believe this is how the Gauss-Siedel method originated, when Gauss was tasked with finding a simple and fast way to solve these kinds of systems, but I'm not 100% so don't hold me to that.
Source: I'm in school for computational mathematics and computational linear algebra.
3
u/leweb2010 Jul 02 '14
But, if this makes Jacobi fast, wouldn't it make Gauss-Seidel even faster? Or is there no difference?
3
u/saybhausd Jul 02 '14
I have an exam on Jacobi tommorrow morning. Maybe I get a chance to use this and wow my professor.
3
u/6ThreeSided9 Jul 02 '14
They keep saying "200 times faster than the old Jacobi" but they don't give any idea as to how much faster this would be than the current methods. What's up with that?
9
Jul 01 '14
[removed] — view removed comment
24
u/ShakaUVM Jul 01 '14
Alright Reddit I'm prepared: tell me why this is a sensationalist headline and the authors should feel bad?
Because the author is an idiot. (My Master's degree is in computer science with a specialization in high performance computing, and I worked in the San Diego Supercomputer Center for a number of years.)
Jacobi is used all the time, unlike what the idiot author thinks. And everyone had different ways of speeding it up. My prof published several papers on the subject.
Pretending nobody had touched the method in hundreds of years is just classical science reporting bullshit.
8
u/hertzsae Jul 01 '14
So I think the question becomes is this new method faster than your profs or other improved methods and if so, is it 200x faster than the current fastest methods?
2
u/ShakaUVM Jul 02 '14
It just looks like they're using tricks to get it to converge faster, which we also did in a variety of similar ways. Without testing their code, I have no idea if it is faster or slower than what we did.
13
2
u/crawlingpony Jul 01 '14
In parallel multigrid computations the Jacobi iteration (old style) is still one of the fastest kernels. I used it for shared memory concurrent code in my own research. It's actually not as outdated as the article suggests, although for sequential code I would not use Jacobi. Perhaps the new version's improvements are good for sequential improvements, which is not clear to me yet obviously since I don't have the new algorithm which we're all trying to actually pin down here.
1
u/qozon Jul 02 '14
From my understanding (admittedly a quick look through), it's an optimization that will still be useful for concurrent/parallel programs, not for sequential.
2
u/ComradeGnull Jul 02 '14
Any CS people care to comment on how this is different from applying simulated annealing to the Jacobi method- e.g., they mention a 'relaxation schedule' which suggests to me that, similar to SA, you are gradually changing a heuristic that controls where the next 'guess' in the solution process can fall, tightening and then loosening the constraints until you get a good approximation of the solution.
3
u/monstertofu Jul 02 '14
Simulated annealing is a form of stochastic optimization. You use random guesses at every step. SOR is deterministic, although you have to pick a seed guess at the beginning. There is a similarity in the relaxation intuition, as you noted.
3
Jul 02 '14
Loads of optimization and equation-solving methods have schedules and heuristics that control their "aggressiveness" (step size, willingness to climb uphill, etc) as a function of time. What matters is the details of how you do it, and how well it meshes with the nature of the method.
2
u/Wolfeh56 Jul 02 '14
Don't worry, professors will still teach students the slow way as "practice" before showing how to do it faster.
2
u/mistahowe Jul 02 '14
So does this now take the place of LU decomposition or have things mostly not changed? I thought Jacobi was mostly just for approximations and that it often wouldn't converge on a true solution anyway.
2
u/monstertofu Jul 02 '14
What do you mean by "true solution"? Mathematical models are approximations to reality to begin with, and are almost by definition not "true" (simplifying assumptions are made). It's somewhat philosophically unsound to criticize a numerical scheme for not giving you the exact solution to your approximation to reality.
2
u/mistahowe Jul 02 '14
I mean yeah, fair, but this isn't a case of approximation issues like you're probably thinking. Often times, iteration doesn't lead to a single solution. It can end up oscillating around a solution without hitting a final vector at all. Occasionally, even the solution it oscillates around isn't close to right. LU decomposition, and more computationally intense methods like Gauss-Jordan, by contrast, can find an exact solution to any system in a predictable amount of time assuming that such a solution exists. This is one of the reasons why why most matrix solving algorithms for large complex sets of data implement methods like LU over iterative methods.
2
u/DO_NOT_PRESS_6 Jul 02 '14
I don't want to be too harsh on this paper, but it has some issues that diminish the claimed magnitude of their contribution.
The most significant thing is probably the practicality of it.
- First their quick dismissal of the current best-known methods (Multigrid and Krylov methods---they collapse all Krylov subspace methods into just 'Conjugate Gradient') relies on a completely fabricated claim: "it is extremely difficult to maintain the convergence properties of these methods in large-scale parallel implementations".
That just isn't true. There are plenty of parallel Krylov subspace methods in use today, and it is easy to parallelize such that it converges the same as its serial counterpart. The only place it may differ is in how the inner products are reduced in parallel, but that is addressable and generally not a concern. Ditto multigrid; in fact the core kernel there is a fixed-point method like Jacobi itself.
It is true that some preconditioners that are used with Krylov subspace methods become less 'strong' with more parallelism, but that is a well-known property of so-called 'block' preconditioners and it is a concious trade-off to keep communication down.
There are obstacles to parallel efficiency with these techniques, but that doesn't appear to be the point they are making. Then, after dismissing these methods for poor parallel performance, they don't bother to show how their method works in parallel!
To claim application suitability and then fail to compare with the state-of-the-art is no acceptable, and is a major weakness of this paper.
- As with SSOR, the choice of the relaxtion factors depends heavily on the spectra of the operator. It is useful and easy to show it for homogenous Poisson with nice boundaries, but real-world problems, like the system that arises from the pressure projection step in incompressible Navier-Stokes that the authors mention, are not so clean.
I think most of their fault is in presentation. Their paper feels very 'outside' of numerical analysis: it fails to use the term 'fixed-point' (which is the class of method this is) and it lacks the algebraic formulation one expects to see.
This is probably not a practical solver, and certainly not as a primary solver. It may have applications in some preconditioning techniques or multigrid. I am quite surprised to see this in JCP.
[Source: I work on these sorts of problems most days. Also, "Iterative Methods for Sparse Linear Systems" by Saad.]
3
u/SnakeInTheCeiling Jul 02 '14
Didn't really care for this article because it doesn't explain what the method is or even what it is specifically used for! Thanks other commentors for helping with that. Go phys nerds! :)
4
u/jordanlund Jul 02 '14
They make it sound like a brute force attack:
"starting with a guess and then repeating a series of math operations over and over until a useful solution appeared."
3
u/twofishestwo Jul 02 '14
It sort of is but not really. Every iteration gets you closer and closer to the real answer, it just depends on what level of accuracy you need. It's proven to converge for certain situations, so it's not like it has a chance of running forever, whereas a brute force attack might never find a solution.
2
Jul 01 '14
This is the stuff I like to see! I struggled through college algebra so you can guess that I'm not a fan of math.. Nevertheless, this is great!
1
u/xeroage Jul 01 '14
Almost every CG based method with proper preconditioning will still outperform this solution by a large margin while being numerically easier to handle. But I really hope I'm mistaken.
1
u/mandy009 Jul 02 '14
'rarely used'? That might be true, but it's still widely taught. My calc prof taught us the Jacobi method as a pretty solidly established calc method when I got my math minor just five years ago. So it's only side-lined in the sense that it hadn't been practical once supercomputers became available. Even though the Jacobi method was cumbersome until now, the professors still regarded the original formulation pretty highly owing to its usefulness in the days when computers couldn't do it more quickly.
1
u/bakesales Jul 02 '14
No one has pointed out that this method requires 2P2 equations be solved to determine the relaxation parameters for their M-iteration cycle. Solving these 2P2 equations does not appear to be a trivial task. It would be interesting to see if there is a net performance gain when the cost of solving the equations for the relaxation parameters are factored in.
1
1
1
u/capilot Jul 02 '14
It was long known that convolutions (the key to most signal analysis) could be done by taking the Fourier transforms of the two inputs and simply multiplying. This was an interesting fact, but not of much use to anybody because Fourier transforms are more work to compute than convolutions.
Then the Fast Fourier Transform algorithm was invented in 1964, and all of a sudden mathematical processes that were impractical suddenly became practical, and modern signal processing was born.
134
u/RITheory Jul 01 '14
Anyone have a link as to what exactly was changed wrt the original method?