A Journey Through Multiple Global Minima
Introduction
Optimization problems don’t always have a single “best” answer. Sometimes, a function has multiple points that are all equally optimal — and this is exactly the fascinating property of Himmelblau’s function, a classic benchmark in numerical optimization.
In this post, we’ll explore Himmelblau’s function, find all four of its global minima using Python, and visualize the results with a stunning 3D surface plot.
What is Himmelblau’s Function?
Himmelblau’s function is defined as:
$$
f(x, y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2
$$
This function is widely used to test optimization algorithms because it has four identical global minima, all with a function value of exactly 0. The minima are located at approximately:
$$
(3.0, 2.0), \quad (-2.805118, 3.131312), \quad (-3.779310, -3.283186), \quad (3.584428, -1.848126)
$$
Because there are multiple equally good solutions, a naive gradient-descent-style optimizer starting from a single point will only find one of these minima — the one closest to its starting position. To find all four, we need to try multiple starting points across the search space.
Strategy
Our approach:
- Define Himmelblau’s function in Python.
- Use
scipy.optimize.minimizewith a multi-start strategy — launching the optimizer from many different initial points across the domain. - Cluster the resulting solutions to identify the distinct global minima (removing duplicates found from different starting points).
- Visualize the function as a 3D surface and mark the discovered minima.
- Also show a 2D contour plot for a clearer top-down view.
Since running the optimizer from a single point is fast, but running it from hundreds of starting points could be slow if done naively, we vectorize the initial point generation with NumPy and use SciPy’s efficient BFGS-based solver (L-BFGS-B) for speed, while keeping the total number of starts modest (a grid of 100 points) so it finishes almost instantly.
Full Source Code
1 | # ========================================================== |
Number of distinct global minima found: 4 Coordinates of global minima: Minimum 1: x = 3.000000, y = 2.000000, f(x,y) = 0.00000000 Minimum 2: x = -2.805118, y = 3.131313, f(x,y) = 0.00000000 Minimum 3: x = -3.779310, y = -3.283186, f(x,y) = 0.00000000 Minimum 4: x = 3.584428, y = -1.848127, f(x,y) = 0.00000000
Code Walkthrough
1. Defining the function and its gradient
The himmelblau() function directly implements the mathematical formula shown earlier. We also manually derived the gradient (himmelblau_grad) — the vector of partial derivatives with respect to $x$ and $y$:
$$
\frac{\partial f}{\partial x} = 4x(x^2 + y - 11) + 2(x + y^2 - 7)
$$
$$
\frac{\partial f}{\partial y} = 2(x^2 + y - 11) + 4y(x + y^2 - 7)
$$
Providing the exact gradient (instead of letting SciPy estimate it numerically) makes the optimizer converge faster and more accurately, since it avoids the overhead of finite-difference approximation.
2. Multi-start optimization
Since Himmelblau’s function has four global minima, starting the optimizer from just one point would only ever find one of them. To solve this, we:
- Create a 10×10 grid of starting points spanning the domain $[-6, 6] \times [-6, 6]$ — 100 starting points in total.
- Run
scipy.optimize.minimizewith the L-BFGS-B method (a fast quasi-Newton algorithm well suited to smooth, bounded problems) from each starting point. - Keep only results where the function value is essentially zero (
result.fun < 1e-6), confirming we’ve truly hit a global minimum rather than some other stationary point. - Deduplicate results: since many nearby starting points converge to the same minimum, we check the Euclidean distance between new solutions and previously found ones, discarding near-duplicates within a
1e-4tolerance.
This grid-based multi-start approach is a simple but effective way to perform global optimization using a fundamentally local optimizer — and because L-BFGS-B is very efficient, all 100 optimization runs complete in a fraction of a second.
3. Building the surface data
We create a fine mesh grid (200 × 200 points) over the same domain and evaluate the function at every point using vectorized NumPy operations. This gives us the Z array needed to draw a smooth 3D surface and contour map.
4. Visualization
- Left panel (3D surface): Shows the overall “landscape” of the function, with two tall peaks and four valley-like basins where the function dips to zero. The red dots mark the discovered global minima, sitting exactly at the bottom of each basin.
- Right panel (2D contour): A bird’s-eye view of the same landscape using color gradients (dark = low value, bright = high value). The four star markers, each labeled with its coordinates, make it immediately clear where all four minima are located relative to each other.
Together, these two views make it intuitive to see why Himmelblau’s function is such a popular test case: the four minima are spread across very different regions of the search space, forcing any global optimization algorithm to genuinely explore rather than just “roll downhill” from a single guess.
Conclusion
Himmelblau’s function beautifully illustrates a key challenge in optimization: not all problems have a unique answer. By using a multi-start strategy with SciPy’s L-BFGS-B solver and carefully deduplicating results, we successfully located all four global minima efficiently. The combination of 3D surface and 2D contour plots gives a complete, intuitive picture of the function’s structure — turning an abstract equation into something you can literally see and understand at a glance.
