What is too big for standard linear algebra/optimization methods?

Question

Different numerical linear algebra and numerical optimization methods have different size regimes where they're a 'good idea', in addition to their own properties. For example, for very large optimization problems, gradient, stochastic gradient and coordinate descent methods are used instead of Newton or Interior Point methods because you don't have to deal with the Hessian. Similarly, dense linear solver methods stop being feasible after a certain size.

So given that both the algorithms and the computer hardware are changing constantly, what's a good way to know, and keep up with, how big is too big for standard linear algebra and optimization solvers?

(I'm thinking about this because when you're an end-user of numerics algorithms it's important to have a vague idea when those algorithms can be applied. Part of that is problem structure and the sort of solution desired, but part of it is also just the size of the problem.)

EDIT: For more concreteness, what got me thinking about this was varying rules of thumb on the upper bounds for how large a problem interior point algorithms could solve. Earlier papers said the dimensionality should be around 1000 while later papers had revised upwards to 5000 and even more recent papers allow for even larger depending on if you can take advantage of sparsity. That's a rather large range, so I'm curious what is large for state of the art interior point methods.

Answer 1

If sparsity is preserved, optimal preconditioners are available, and inequality constraints can be resolved by a multiscale method (or the number of active constraints is not too large), the overall algorithm can be $\mathcal O(n)$ time and space. Distributing across a parallel machine adds a logarithmic term to time. If enough sparsity is available, or if matrix-free methods are used, on the order of one million degrees of freedom can be solved per core. That puts the problem size for today's largest machines at around one trillion degrees of freedom. Several groups have run PDE simulations at this scale.

Note that it is still possible to use Newton-based optimization with large design spaces, you just need to solve iteratively with the Hessian. There are many approaches to doing so efficiently.

So it all depends how you define "standard methods". If your definition includes multilevel structure-preserving methods, then extremely large problems are tractable. If your definition is limited to unstructured dense methods, the feasible problem sizes are much smaller because the algorithms are not "scalable" in either time or space.

Answer 2

The limit is primarily given by the memory it takes to store the matrix representation, and the time to retrieve it from memory.

This makes a few thousand the limit for dense matrix method in simple environments. For sparse problems the limit for direct solvers is much higher, but depends on the sparsity pattern, as fill-in must be accomodated. The limit for iterative methods for linear solvers is essentiall the cost of a matrix vector multiply.

The limit for solving the linear subproblems directly translates into corresponding limits for local solvers for nonlinear systems of equations and optimization problems.

Global solvers have much more severe limits, as they are limited by the number of subproblems that need to be solved in a branch and bound framework, or by the curse of dimensionality in stochastic search methods.

Answer 3

For a given class of problem, good ways to figure out what constitutes how big is "too big" are to:

• search the literature in your field; understand what factors will limit execution time (these include memory usage and communication)
• look at test problem sets (better for serial algorithms than parallelized algorithms, since I don't know of many test problem sets for parallelized algorithms)
• break the algorithm down into smaller, constituent parts, and use what you know about the parts to make inferences about the whole (example: Krylov subspace methods boil down to loops of matrix-vector products, so then you have to figure out how many iterations is "too many" and what makes a matrix-vector product "too big")
• try it out yourself (on a serial machine, or if you have the resources, a parallel machine). See how large of an instance you can run, or run several instances of different sizes, and do an empirical scaling analysis (time vs. problem size).

To give an example, in global optimization, the concrete answer is extremely structure-dependent. As Arnold Neumaier notes, deterministic global optimization algorithms tend to be limited by the number of subproblems that must be solved in a branch-and-bound (or branch-and-cut) framework. I have solved mixed-integer linear programs (MILPs) containing thousands of binary variables, but I suspect that the reason I can solve such large problems (comparatively speaking, for MILPs) is because the problem structure was such that few subproblems were required to solve some critical set of binary variables, and the rest could be set to zero. I know that my problem is "large"; I've constructed other MILPs of similar size that solve tens of thousands of times more slowly.

There are global optimization test sets that give you an idea of what is "run-of-the-mill," and the literature can give you ideas of what problems are "large". Similar tactics exist for figuring out the state-of-the-art in problem sizes in other fields, which is how Jed Brown and Arnold Neumaier can quote these figures. It's great to get these numbers, but it's far more valuable to be able to figure out how to get them yourself when the time comes.