# stupid + stupid = brilliant in scientific computing

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.

• Not really an example for "brilliant", but regularly seen: solving Hamiltonian problems using the Euler method, implementing the updates wrong and getting the symplectic Euler method, which effectively is the second order Verlet method with slightly different initial points. This gives, e.g., closed planetary orbits with small energy oscillations while the correct implementation of Euler would give rapidly diverging spirals. Sep 16 at 20:21
• I'm not sure if you consider this scientific computing and the results aren't that brilliant, but a neural network XML-CNN performed well, despite having a humongous bug in the code, due to a series of extremely fortunate (for them) coincidences. The original paper was highly cited and the bug wasn't discovered for many years. See this. Sep 17 at 5:46
• A few years ago there was a paper I think about a neural network that had no actual non-linear activation function, but it used values close to 0, and the irregularities of what is representable as a floating-point and what isn't, close to 0, was enough nonlinearity for the neural network to perform well.
– Stef
Sep 17 at 18:08
• Also, a random forest typically performs much better than a decision tree, even though a random forest is a collection of suboptimal decision trees.
– Stef
Sep 17 at 18:09
• @Stef more recent but I think this matches what you describe (in both paper and hourlong video format) tom7.org/grad/murphy2023grad.pdf youtube.com/watch?v=Ae9EKCyI1xU Sep 18 at 16:42

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;
}

• This is a long shot, but are there libraries/compilers that will do this automatically? i.e. take a function as an input (perhaps as a string), identify values where catastrophic cancellation or other numerical issues can occur, and then return a more stable form of the function considering the rules of floating point arithmetic?
– NNN
Sep 20 at 3:31
• @NNN Not to my knowledge. There is an open source tool called Herbie that can re-arrange floating-point expressions in numerically beneficial ways. Unfortunately, despite the fact that the authors have put in a lot of work over the years, the solutions it generates are usually mediocre and not competitive with expert human efforts, sometimes comically inept, and have matched my manual efforts only twice over five years and several dozen test cases. Your mileage may vary, so give it a try. Sep 20 at 4:38

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.

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.

• "gradient descent [...] is comically slow in many common problems" << Nothing is slow in the absolute sense; gradient descent is just slower than other approaches that solve the same problems. For instance, using gradient descent to find a minimum of a differentiable function f when you could instead solve equation f'(x) = 0 directly is pretty slow. But there aren't really any rival approaches to optimise the weights of a neural network, gradient descent is pretty much the only approach we have.
– Stef
Sep 18 at 15:27
• @Stef: added a citation for the second claim to the answer. I don't think the first claim ("Deep learning laughs...") since I don't think anyone believes I'm being literal. Sep 18 at 16:13
• @Stef: I don't want to restart the flame war that I cited in my link about "The fact that these methods are so popular in deep learning is confusing...". And while I will say I'm on your side that specialized 1st order (and quasi 2nd order) methods are currently SOTA, this is surprising to folks who think that almost any 1st order full gradient method can be effortlessly improved via LBFGS. So it's not the case that there are no rivals to gradient descent (although stochastic gradient methods are a different beast). Sep 18 at 16:19
• I'm just one data point, but I learned modeling and curve fitting the classical way and had the same reaction as @CliffAB when I learned about deep learning. "There are way too many parameters! Gradient descent is naïve and gets stuck in shallow local minima!" Sep 18 at 17:17
• I would argue that a naive combination of both of these methods would indeed be stupid on many problems, and that regularization and better training procedures are the only things that make neural networks do anything other than overfit Sep 18 at 17:51

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.

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.

• The game of go was considered too complex for computers because, as opposed to chess, there wasn't a simple heuristic evaluation function to tell how favourable a given position is for the player, so even if you explore the tree of possible moves, you can't compare branches and choose the best. And then in 2006 someone introduced "Monte Carlo tree search", which is basically random exploration of the tree. And suddenly computers weren't that bad at go.
– Stef
Sep 18 at 21:10

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).

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).

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.

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.

• Oh, so that's how democracy works. Sep 19 at 17:37
• To make things worse, the "slightly useful models" are sometimes referred to as "experts" in the jargon.
– Stef
Sep 20 at 8:54
• @JoshuaJurgensmeier Of course, an ensemble of worse-than-random agents will perform much worse than any of them individually, so there's that... Sep 20 at 13:07
• IMHO, this is a risky simplification of the reason why ensemble learning works Sep 25 at 20:41

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.