stupid + stupid = brilliant in scientific computing

Question

I'm interested in examples of very effective methods in scientific computing that are the sum or naive combination of very ineffective or bad ones.

Answer 1

In W. Kahan, "Interval arithmetic options in the proposed IEEE floating point arithmetic standard". In: Karl L. E. Nickel (ed.), Interval Mathematics 1980, New York: Academic Press 1980, pp. 99-128, the author addresses some common fallacies with regard to floating-point computation, which he calls anti-theorems. One of these is the notion that if some intermediate results suffer from significant error compared to the corresponding mathematical results, then any final result strongly dependent on these must be wrong as well. He provides the following counterexample:

Assume we wish to compute $f(x) = (\exp(x) - 1) / x$. A novice might compute f (x) := if x = 0 then 1 else (exp (x) - 1) / x. Somebody who has taken an introductory course in numerics would realize immediately that this suffers from subtractive cancellation near zero. Kahan then proposes to fix this as follows: x := exp (x); if x ≠ 1 then x := (x - 1) / ln (x); f(x) := x (one can tell that Kahan comes from an era when both register and memory space were at a premium).

Most people's gut reaction upon seeing this for the first time is that this supposed fix is preposterous, as it compounds a bad idea with an even worse idea. Yet it works like a charm for $x$ near $0$. An example using IEEE-754 binary32 format (single precision), using $x = 3 \cdot 2^{-24} \approx$ $1.78813934 \cdot 10^{-7}.$ exp(x) - 1 evaluates to $2^{-22} \approx 2.38418579 \cdot 10^{-7}$ because of catastrophic cancellation. log (exp (x)) evaluates to $2.38418551 \cdot 10^{-7}$. (exp (x) - 1) / (log (exp (x)) thus produces a final result of $1.00000012$, which is basically accurate to single precision. The naive computation (exp(x) - 1) / x would instead produce a result of $1.33333337$.

Kahan simply explains that the trick consists in the cancellation of errors. But it took deep insight to come up with this solution. It is therefore not surprising that William Kahan subsequently became the “father of IEEE floating-point arithmetic”. A detailed numerical analysis of this algorithm is provided by Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, Second Edition, SIAM 2002, section 1.14.1 on pp. 19-21.

The practical relevance of Kahan's counterexample is that it extends trivially to the accurate computation of $\exp(x)-1$ in the form of a function expm1(), which is frequently needed in finite-precision floating-point computation to avoid cases of subtractive cancellation. See example C code at the end. In the 1980s exp() and log() implementations were not necessarily faithfully rounded, but they were usually among the most accurate, well-tested, and fastest elementary math functions provided by computer math libraries, making expm1() implementations based on Kahan's algorithm sufficiently robust (compared to computing expm1 (x) = 2 * sinh (x/2) * exp (x/2), for example).

The ISO-C99 standard added the expm1() function to the standard C math library, and this was subsequently incorporated into the ISO-C++11 standard. A perusal of the Fortran 2008 standard shows that this function is not provided. There is a similar trick for the accurate computation of $\ln(1+x)$, usually called log1p(), from log(). It first appeared in Hewlett-Packard HP-15C Advanced Functions Handbook, Hewlett-Packard 1982, Appendix: Accuracy of Numerical Calculations, p. 193. Style and contents of this appendix strongly suggest that Kahan provided it.

/*
  The following is based on: W. Kahan, "Interval arithmetic options in the 
  proposed IEEE floating point arithmetic standard". In: Karl L. E. Nickel 
  (ed.), "Interval Mathematics 1980", New York: Academic Press 1980, pp. 
  99-128. See pages 110-111. For a detailed explanation see Nicholas J. 
  Higham, "Accuracy and Stability of Numerical Algorithms, Second Edition",
  SIAM 2002. In section 1.14.1 on pages 19-21.
*/
double my_expm1 (double x)
{
    double u, m;

    u = exp (x);
    m = u - 1.0;
    if (m == 0.0) {
        // x very close to zero
        m = x;
    } else if (fabs (x) < 1.0) {
        // x somewhat close zero 
        u = log (u);
        m = m * x;
        m = m / u;
    }
    return m;
}

Answer 2

Consider the constrained optimization problem $$\min_x f(x) \quad\text{s.t. } g(x) = 0$$ where, to make things nice, we'll assume $f$ is convex and $g$ has convex level sets. There are two bad ways to do this: the quadratic penalty method and the Lagrange multiplier method. The quadratic penalty method seeks an unconstrained minimizer $x_\alpha$ of the objective $$f_\alpha(x) = f(x) + \frac{\alpha}{2}\|g(x)\|^2$$ and takes the limit as $\alpha \to \infty$. Provided that $\alpha$ is large enough, this problem is usually convex, which is nice. But this approach suffers from very bad conditioning in the limit as $\alpha$ grows large though. Alternatively, one can use the Lagrange multiplier method, where we instead seek a critical point of the Lagrangian $$L(x, \lambda) = f(x) - \langle g(x), \lambda\rangle.$$ We've now got a saddle point problem, it can be difficult to prove any estimates on the norm of the multipliers $\lambda$, or to guarantee that the updates reduce the objective or deviation from the constraint.

The augmented Lagrangian method is the rather naive combination of the penalty and Lagrange multiplier methods: $$L_\alpha(x, \lambda) = f(x) - \langle g(x), \lambda\rangle + \frac{\alpha}{2}\|g(x)\|^2.$$ Unlike the pure penalty method, we can still obtain convergence to a solution without taking $\alpha \to \infty$. We can do so using a very silly gradient ascent update for $\lambda$ and, for sufficiently large $\alpha$, we get explicit guarantees on the convergence rate of the Lagrange multipliers. The reasons for this are explained in Bertsekas (1982).

In short, penalty = stupid, Lagrange multipliers = stupid, augmented Lagrangian = penalty + Lagrange multipliers = brilliant.

Answer 3

Disclaimer: stupid is in the eyes of the beholder.

For decades, it seemed like what would now be known as deep learning was combining stupid with stupid.

Stupid 1: over parameterization. A rule of thumb is that the more parameters you put in your model, the more unnecessary noise your model can pick up and the less reliable it will become. Deep learning laughs at this rule and in some cases has more parameters than data points.

Stupid 2: gradient descent. This is the first optimization algorithm you will learn in a numerical analysis course and is comically slow in many common problems. The fact that these methods are so popular in deep learning is confusing to many who are classically trained in numerical analysis (including myself for a good period of time!). Switching to stochastic gradient descent was better justified approach, but the question "have you tried gradient descent?" felt like "I watched my first online course and now I think I know it all".

And yet...these approaches have completely revolutionized the problem of prediction from data, so much so that they are being used to revolutionize almost every other field of science.

Answer 4

B-tree based data storage:

1. Stupid: Binary search tree. Paths are too long, need to fetch too many nodes on lookup.

2. Stupid: Linear array. Array gets too long, need to fetch too many nodes on lookup.

1+2: Search tree of linear arrays. Large branching factor constrains depth, short array length allows quick loading of each node.

Answer 5

The usage of random numbers in numerical calculations can be considered as something in that category. A naive view of "random" is that it is something "lacking a definite plan, purpose, or pattern" (Webster). On the other hand, Monte-Carlo methods, based on random numbers, have become a huge and highly successful field, in statistical physics, risk analysis etc.

Answer 6

This one comes from polynomial interpolation, however, I don't know whether the ingredients fully meet the assumption of being stupid.

1. Polynomial interpolation of a function $f(x)$ in the monomial basis is usually a stupid idea. The task there is to find an approximation $\tilde f(x) = \sum_{k=1}^N \alpha_k x^k$ such that $f(x)=\tilde f(x)$ holds at $N$ selected and distinct gridpoints $\{x_k\}_{k=1}^N$. The direct solution of this problem leads to Vandermonde matrices which are in typically exponentially ill-conditioned, and this translates into the fact that one can't find accurate parameters $\alpha_k$.

2. The Lanczos-Arnoldi algorithm and other Krylov subspace methods are actually ingenious, as they allow the approximate orthogonalization of large matrices, but for the present thread it could be considered stupid (at least when compared to other dense methods such as QR, SVD, etc.): Given a quadratic matrix $\mathbf G$, a (usually) random vector $\mathbf v$ and a positive integer $M$ that is usually much smaller than the matrix dimension, it produces a from-the-first random basis $$ \Big\{\mathbf v, \mathbf G \mathbf v, \ldots , \mathbf G^M \mathbf v \Big\}, $$ Gram-Schmidt orthogonalizes these vectors on the fly, and then diagonalizes the matrix $\mathbf G$ in this basis.

Now, the brilliant magic happens when you combine these two approaches, as has been discovered in a 2019 paper by Trefethen and coworkers. By reformulating the polynomial interpolation problem as an Arnoldi algorithm (using $\mathbf v = \mathbf 1$ and $\mathbf G = \mathbf X$, the diagonal matrix with the gridpoints $x_k$ on the diagonal), they show that the previously ill-conditioned problem becomes stable, yielding highly accurate solutions to the interpolation problem. The basic interpretation is that the Arnoldi process with its Gram-Schmidt orthogonalization provides a convenient basis of orthogonal polynomials that give similar benefits as the classical orthogonal polynomials. For details, please see the linked paper.

Answer 7

Boosting in machine learning consists of iteratively learning weak (i.e. "stupid") models and combining them all into a stronger final model. If performed correctly, the final model is typically far superior to any of the individual weak models (which individually may only be slightly better than random guessing).

Answer 8

Ensemble learning relies on the notion that combining even slightly useful models can produce a near-perfect one. Suppose you have a binary classification problem and a classifier with 51% accuracy in a balanced test set. If you combine many different classifiers that each have 51% accuracy and take a majority vote among them, you will have accuracy higher than 51%. As you use more and more models, the chance that over half of them are wrong becomes arbitrarily small. You can use lots of "stupid" models that perform barely better than random guessing, and find that together, they are almost always correct.

Answer 9

Another example comes from mixed finite elements. Let $V$ be a Hilbert space, $u$ an element of $V$, $Q$ another Hilbert space and $p$ an element of $Q$. Finally, let $A : V \to V^*$ be self-adjoint and positive-definite, and let $B : V \to Q^*$ be bounded. The problem is to find a saddle point of the functional $$L(u, p) = \frac{1}{2}\langle Au, u\rangle - \langle f, u\rangle + \langle Bu - g, p\rangle - \frac{\epsilon}{2}\langle p, p\rangle$$ where $f$ and $g$ are the input data and $\epsilon$ is small (possibly 0). Examples include the incompressible Stokes problem of viscous fluid flow and nearly incompressible elasticity.

The conventional approach to solving such problems is to find a pair $V_h$, $Q_h$ of finite-dimensional spaces satisfying a discrete inf-sup condition: $$\inf_{q\in Q_h}\sup_{v\in V_h}\frac{\langle Bv, q\rangle}{\|v\|\|q\|} \ge \beta$$ for some $\beta$ that does not depend on the fineness $h$ of the discretization. We then use a conventional Galerkin discretization of the original problem.

Finding good element pairs is difficult, and many reasonable-looking elements are unstable (see Boffi, Brezzi, and Fortin). Some element pairs result in linear systems that have no solution, while others give systems that are underdetermined. For example, the CG1 velocity / DG0 pressure pair for the Stokes equations is not inf-sup stable. There are circumstances, however, where an unstable element pair can be resuscitated by integrating the constraint terms ($\langle Bu - g, p\rangle$, $\epsilon\langle p, p\rangle$) inexactly. In other words, a sensible person would choose a quadrature rule of high enough accuracy to integrate all terms in the variational problem exactly, which you can do if you know what degree of finite element basis was being used and something about the operator $B$. Choosing to intentionally under-integrate those terms can, with some choices of basis that are not stable when integrated exactly, result in a stable discretization. See $\S$8.12 of Boffi, Brezzi, and Fortin, or Oden and Kikuchi (1982) for more details.

In short, unstable basis = stupid, inexact quadrature = stupid, unstable basis + inexact quadrature = brilliant (but only sometimes).

Answer 10

Along the lines of linear algebra, there are a couple related things worth mentioning. First, I often describe Jacobi's method as typically fitting the prompt "the sum or naive combination of very ineffective or bad ones." Why? Because Jacobi's method involves

1. Solve a multivariable equation for a single variable
2. Plug a result back into itself as if looking for a fixed point

but, it turns out that if your system is diagonally dominant, then you really have a contraction, and it's provably convergent. It's also one of the most commonly used solvers/preconditioners because it's so simple.

With that said, Jacobi's method is brilliant in one way, but also simple in another way. If you take

1. Jacobi's method, which is simple because it converges quite slowly with respect to iterations
2. Interpolation, which on its own does very little for solving linear systems

and if you combine them, you end up with the multigrid method, which has an optimal computational complexity for convergence and is one of the most popular methods for solving extremely large linear systems.

Answer 11

Gauss-Seidel solvers can be presented as Jacobi iteration solvers where one has made the "bug" of overwriting into a single vector instead of using two separate buffers as in a correct Jacobi iteration solver. The basic Jacobi iteration only guarantees convergence with diagonal dominance, but somehow miraculously the Gauss-Seidel solver (a buggy Jacobi solver) now converges for arbitrary symmetric positive definite matrices.