I am new to linear programming and have formulated a linear program (LP) with order of $10^{13}$ variables and $10^{13}$ constraints, although the constraint matrix is extremely sparse.
I wanted to know if an LP of this scale is tractable or not?
$\begingroup$
$\endgroup$
-
$\begingroup$ I presume you mean $10^{13}$ variables and $10^{13}$ constraints? $\endgroup$ – Geoff Oxberry Jun 15 '12 at 11:22
-
$\begingroup$ Yes, sorry about the bad formatting. $\endgroup$ – stressed_geek Jun 15 '12 at 11:26
$\begingroup$
$\endgroup$
$10^{13}$ is beyond tractable on today's largest supercomputers due mostly to memory limits. The largest problems I've seen practically solved have had on the order of $10^5$ rows and columns, but the most important factor tends to be the number of nonzeros, where we are just crossing into solving problems with $10^6$ nonzeros. See Mittelman's parallel benchmarks page to get a feel for what the best freely available and commercial solvers can do on a range of problems at this size.
answered Jun 16 '12 at 12:40
$\begingroup$
$\endgroup$
A huge number of inequality constraints is generally tractable by constraint generation techniques, processing at each time only a limited number of constraints, and adding constraints violated at the current solution to the constraint set (while deleting strongly satisfied ones). But in the cases I have seen, this requires that the number of non-slack variables is limited, too.
So the question is whether you can reformulate your problem in such a way that the problem itself or its dual has many fewer variables than constraints.
Edit: See also recent work by Nesterov on ''Subgradient methods for huge-scale optimization problems'', http://dial.academielouvain.be/vital/access/services/Download/boreal:107876/PDF_01 . His technique works under favorable circumstances with a complexity of $O(n\log n)$ if the accuracy requirements are moderate.
$\begingroup$
$\endgroup$
My comment on Aron's answer got too long, but I would like to augment his answer:
I like bringing up the example of parallel benchmarks. Several comments to add here. All four of the solvers tested are commercial, but have free academic licenses available. Furthermore, the test cuts off the run time at 25000s, or ~8h, which is arbitrary, and depends heavily on external resource constraints (i.e., in a company, people may be unwilling to wait more than a day for results; in academia, some people can run their simulations for months). The test is run on a single quad-core machine, which is not indicative, to me, of bleeding edge performance.
When I was hunting around for data to answer this question, I found papers 10+ years ago that were solving problems of roughly the same size, which suggests to me that we might be able to do better now with the infrastructure we have. Certainly not $10^{13}$ variables, but based on the $n^{3.5}$ scaling of interior point methods, if the linear algebra and parallelism were implemented well, and you had the time and a modestly sized cluster, I don't see why you couldn't attempt to solve a problem with $10^{7}$ or maybe $10^{8}$ variables (only if you had special structure you could exploit with decomposition methods like Benders' decomposition or Dantzig-Wolfe, plus cutting-plane generation algorithms). (I will add that I am ignoring the effect of constraints, which complicate matters depending on how many bits are stored; this effect only makes the estimates below more pessimistic.)
I believe GAMS has a parallel implementation, and since it uses solvers like CPLEX, Gurobi, MOSEK, and Xpress (i.e., the four solvers in the benchmark Aron cites), I don't see why one couldn't run parallel instances of those solvers (in fact, I know this is possible for CPLEX and Gurobi). The limiting factors will be memory (because memory is expensive) and communication more than processing power, since a linear program reduces, ultimately, to constructing and solving a system of linear equations repeatedly (a massive oversimplification, yes, but linear algebra is something we know how to parallelize).
But $10^{13}$ variables is too much. Assuming that memory and communication were no object, you'd need to take the largest problem in the benchmark, and scale the run time on that machine by a factor of roughly $10^{24}$ before possibly exploiting special structure that your particular problem may have. That's not to say that you couldn't try to solve it approximately using the methods that Professor Neumaier has suggested, but a solution to optimality is likely impossible without waiting a really long time, using a huge computer, and having a scalable implementation of an LP solver tuned to that huge computer.
As a very rough order of magnitude estimate, the Intel Core 2 Quad used in the benchmark Aron cites can operate at a peak speed of 40 gigaflops. Assuming you were to get on Jaguar, Oak Ridge National Lab's supercomputer, and you could use the whole machine (extremely unlikely, but let's dream big), you'd have roughly 2 petaflops at your fingertips (based on the TOP 500 numbers here, or roughly 50000 times the computing power, not counting penalties due to communication, memory limitations, or the fact that no one ever runs at peak speed (not even the LINPACK benchmark).
Going from $10^{6}$ to $10^{7}$ variables means roughly a factor of 10,000 increase, which you could conceivably split among a cluster of 50-100 machines, and waiting a month (assuming you're willing to wait, you have the machines, and again, memory and communication aren't limiting, all of which are big "ifs"). Going from $10^{6}$ to $10^{13}$ variables means going from your desktop to using all of Jaguar, and waiting roughly $10^{17}$ to $10^{18}$ years. (And again, these effects also ignore the fact that you're going to have more constraints!)
answered Jun 16 '12 at 15:58
-
$\begingroup$ I was going to start off with an estimate of how long it would take on a supercomputer, but one of the problems is that many of the faster LP solvers rely on a Cholesky decomposition, and that will eventually blow up memory usage due to fill-in. Also, these machines implement shared-memory parallelism, but true distributed-memory parallelism is, for whatever reason, somewhat rare in commercial LP solvers right now. $\endgroup$ – Aron Ahmadia Jun 17 '12 at 5:21
-
$\begingroup$ The $O(n^{3.5})$ complexity of interior point method is for dense problems only. Commercial sparse LP solvers solve routinely LPs with several million of variables if these have a favorable sparsity pattern. $\endgroup$ – Arnold Neumaier Jun 17 '12 at 11:37
-
$\begingroup$ @ArnoldNeumaier: Agreed. I have heard secondhand that airlines solve problems of that scale. The analysis above is an order-of-magnitude estimate, of course, and I've only ever seen the computational complexity analysis for the dense case. I suspect that the sparse case has a factor of at least $O(n^{2.5})$, in which case the dire estimates of time required above would decrease accordingly, but would still remain unacceptably high. $\endgroup$ – Geoff Oxberry Jun 17 '12 at 12:47
-
$\begingroup$ @AronAhmadia: Agreed. All I see in solver manuals are dense or sparse Cholesky factorizations used to solve the barrier method system. No one uses Krylov subspace methods at all, nor have I seen distributed memory (i.e., MPI-based) parallelism. I don't see optimization problems (with the possible exception of PDE-constrained problems) solved on supercomputers. Maybe I'm not looking at the right papers. $\endgroup$ – Geoff Oxberry Jun 17 '12 at 12:50
-
$\begingroup$ The worst case complexity of interior point methods is $O(\sqrt{n})$ iterations, but in practice the number of iterations is typically about 30-50, independent of the dimension (or perhaps growing with its log). Thus in practice, dense IP methods scale with $O(n^3)$, and sparse IP methods with about $O(ne^2)$, where $e$ is the average number of entries in a row of the Cholesky factor. - But of course problems with $n=10^{13}$ are far bigger than anything that has been solved so far. $\endgroup$ – Arnold Neumaier Jun 17 '12 at 13:14
$\begingroup$
$\endgroup$
It's an old question, but hey, 25+ years ago somebody could already solve a 1.1M constraint, 2.6M variable problem in 3h on a PC cluster. http://www.maths.ed.ac.uk/~gondzio/CV/finance.pdf
I would like to see the generating equations, perhaps it's smart to do some serious decomposition before you throw this problem to the algorithms, and I would like to say as a practitioner that maybe it is smart to chew on it before feeding the hardware. Also sounds like the size that would induce to numerical errors in the formulation given limited computer memory and precision.
