Are You Confident Your Optimization Is Really the Maximum?
How GAUSS 27’s Kernel Procs Made Convergence Checking Affordable
I’m not proud to admit this but I usually take an optimization solution as given when I estimate my models. I take a quick glance at gradients, make sure the estimates make theoretical sense, and call it day. If the gradients are small, the iterations have stopped, and the log likelihood looks reasonable, I write it up. Admitedly, this was part laziness and part praticality.
Robust optimization checking requires repeating the optimization from different starting points and comparing the results. That takes time, though. If one run takes several minutes and you have dozens of series to estimate, followed by a bootstrap, trying many starting values can become expensive. It’s tempting to try one or two and move on.
I ran into this problem while estimating a Markov-switching model across a panel series of 77 countries. I wanted to check how much the estimates depended on the starting values, but the time required made an extensive comparison difficult.
This is changing with the upcoming GAUSS 27 release and it’s new compiled kernel procedures. The tools is simple; adding kernel before proc allows GAUSS to compile the procedure to machine code. Combined with vectorization, this made a single likelihood evaluation roughly 15–20 times faster than in my original code.
With the faster code, I could repeat the estimation from many more starting points. In 58 of the 77 countries, those additional runs found higher likelihoods than the original run. The optimizer’s convergence report alone hadn’t revealed those better solutions. Not a result I was happy to see but an important one for improving my model.
Why starting values matter
I’m estimating a Markov-switching model of current account balances using an unbalanced panel of quarterly data for 77 countries. The model allows the average balance, volatility, and persistence of imbalances to change over time. That flexibility makes estimation more challenging, because different combinations of regimes can produce competing explanations of the data. I fit the model separately for each country using maximum likelihood in GAUSS.
Evaluating the likelihood means working through the observations in order. At each quarter, a Hamilton filter predicts the probability of each of the eight regimes from the previous quarter’s probabilities, then updates them using the new observation. Each quarter’s result feeds the next, so the recursion has to stay a loop even when everything else is vectorized, and that loop is where most of the time goes.
The likelihood can also have multiple local maxima. Different parameter combinations may explain the same observations in different ways. For example, a high-variance regime might capture a few unusual quarters or a longer period of volatility. These explanations can lead to different likelihood values, and a gradient-based optimizer can settle at different solutions depending on where it starts.
Simply exchanging regime labels can also produce equivalent solutions. For convergence checking, though, the concern is whether another starting point finds a higher likelihood.
A convergence report tells us that the optimizer met its stopping criteria. It doesn’t establish that a better solution couldn’t be found from a different starting point.
Checking additional starts requires many more likelihood evaluations. In this implementation, each numerical gradient requires about 21 evaluations, and an optimization can take hundreds of iterations. The bootstrap repeats the estimation roughly 1,000 times per country.
In the code I started with, all the calculations were inside the per-quarter loop. A full bootstrap was projected to take 18–22 days of serial computation, using one starting point per estimate. Repeating those estimates from multiple starting points would have extended an already lengthy run.
Making the likelihood faster
To make a broader search practical, I needed to reduce the time spent evaluating the likelihood. This could partly be done just by reorganizing the existing calculations. Beyond that, the new GAUSS kernel feature (coming in GAUSS 27) allowed me to compile the filtering loop and, eventually, the whole likelihood. In this case, the larger percentage gains game purely from vectorization but the kernel feature gave an extra boost and is going to be a game changer moving forward.
Step one: vectorize the density calculations
The original code calculated the eight state densities separately within each quarter of the filtering loop. Those calculations don’t depend on the filter’s updated regime probabilities, so I moved them outside the loop and computed the densities for all quarters at once.
That left only the sequential filtering operations inside the loop:
- Predict the regime probabilities.
- Update them using the current observation.
- Accumulate the log likelihood.
- Carry the updated probabilities into the next quarter.
Across the 77-country workload, this made each log-likelihood evaluation 6.4 times faster. The procedure returning per-observation likelihood contributions, used to calculate gradients, became 5.4 times faster.
As a robustness check, I compared the revised code against the original at 616 parameter vectors. These included the saved estimates for all 77 countries, randomly selected vectors, and deliberately chosen cases that exercise the density cutoff and penalty branches. The estimations were nearly identical with the largest scaled difference equal to 1.4 × 10⁻¹⁵.
Step two: compile the filtering loop
With the density calculations vectorized, I turned to the filtering loop that still had to run one quarter at a time. To use GAUSS’s new kernel feature, I moved that recursion into a self-contained procedure, passing in the matrices it needed rather than accessing global variables.
Once the loop was organized this way, I could compile it by adding kernel before proc. This gave me a direct comparison between ordinary and compiled versions of the same filtering code.
On a simulated example with 200 quarters and a fixed random seed, the results were:
Compiling the loop cut the total optimization time from 0.79 to 0.48 seconds, about 1.6 times faster than the vectorized version. I compared the two versions at 616 different sets of parameter values, and with GAUSS’s default exact mode their outputs were identical down to the last digit.
Get Eric Clower’s stories in your inbox
Join Medium for free to get updates from this writer.
To see where the time savings come from, think of the optimizer as a hiker climbing toward a summit in fog. After each step, the hiker stops to check the map and decide where to go next. Each map check means evaluating the likelihood dozens of times. The total time depends on two things: how many steps the hiker takes, and how long each map check takes.
- Vectorizing mostly made the map checks faster, cutting the total time from 5.32 to 0.79 seconds. It also happened to reach the summit in slightly fewer steps (98 instead of 106), but that was luck rather than design. The two versions agree to about 15 digits. Even differences that small can nudge the optimizer onto a slightly different route, and on real data I saw the step count go up as well as down. Nearly all of the gain came from faster map checks.
- Compiling didn’t change the route at all. Both versions took exactly 98 steps and reached the same summit, with a log likelihood of −269.717143. So the 1.6× speed-up from compiling comes entirely from faster map checks.
On the real 77-country workload, compiling the loop sped up the likelihood calculations by about 1.4 times. That’s less than in the demo because only part of each calculation was compiled: the loop that updates the regime probabilities quarter by quarter. The rest, which computes how likely each observation is under each regime, was still ordinary GAUSS code. Compiling only speeds up the part you compile.
Step three: compile the whole likelihood
An engineer at Aptech suggested going one step further: put everything, the regime calculations and the quarter-by-quarter loop, into a single compiled procedure, written with plain loops throughout.
That runs against the usual GAUSS advice. In step one, I got most of my speed-up by getting rid of loops and working on whole columns at once. But that advice exists because ordinary GAUSS loops are slower than vectorized computations. Once code is compiled, simple loops are fast, and they can even beat the vectorized version because they skip the overhead of building intermediate matrices.
On the real 77-country workload, the fully compiled version made each likelihood calculation about twice as fast as the version where only the loop was compiled, and each optimizer step 1.6 to 2.2 times faster. The benchmark on simulated data showed a larger gain, about 3.5 times, which is a reminder that speed-ups depend on the problem.
This time, the answers were not quite identical. At 548 of the 616 test points they matched exactly. At the rest, they differed only around the 16th digit, because the compiled loops add numbers up in a different order. That’s far too small to matter for any single calculation.
But remember the hiker: tiny differences can still send the optimizer down a different route. And on real data they did. Starting from the same point, the two versions sometimes reached the same summit by a different number of steps, and sometimes ended on different summits altogether. Neither version is wrong. The likelihood surface simply has many peaks, and small differences can decide which one the hiker reaches. That’s why the speed mattered so much in the end, as the next section shows.
What the faster code let me check
The speed-up mattered less for any single estimate than for what it let me check. A run covering all 77 countries with ten different starting points each took 288 seconds. The projected time for the full bootstrap fell from 18–22 days to roughly 1–1.5 days. (That’s still a projection; the bootstrap hasn’t been run yet.)
The ten-start run was sobering. In 58 of the 77 countries, only one of the ten starts reached the highest likelihood found. Going back to our hiker analogy: the fog hides many summits, and where the hiker ends up depends on where they start. So I widened the search to 100 random starts per country plus five “perturbation chains” that repeatedly nudge the best solutions and keep any improvement. That took three to four hours for all 77 countries, which would have been impractical with the original code.
Checking the results along the way turned up three surprises.
Some summits were fake. A few end points had the highest likelihoods of all, but their gradients were around 10⁰⁹: far from flat, so not genuine peaks. In these, one regime had shrunk to fit a handful of quarters almost perfectly, which can push the likelihood up without limit. Picking “the highest likelihood” blindly would have selected exactly these runs. So every end point has to pass the same convergence checks before it’s ranked, and I added a small penalty that stops a regime from shrinking onto a handful of points.
The same run could give different answers. In one country, repeating a run from the same starting point gave different results. The cause was OPTMUM’s recovery step: when a line search fails, it takes a small random step, drawn from GAUSS’s global random-number generator. Its outcome depended on whatever had run before. Setting a seed before each optimizer call, separately for each country and start, made every run exactly repeatable.
The lesson applies to any optimizer, not just OPTMUM: many draw random numbers somewhere you might not expect, such as a fallback step, a random restart or an initial search pattern. If they share a global generator, a run’s result can depend on what ran before it, so results change when you reorder runs, split them across sessions or resume after an interruption. A single seed at the top of a script doesn’t fix this. Seed each optimizer call on its own, and test repeatability by running the same start twice in different orders.
Different code could reach different summits. Even with that fixed, the original and vectorized likelihoods, which agree to about 15 digits, sometimes led the optimizer from the same start to different peaks. With many summits close together, a rounding-level difference early on can be enough to tip the hiker onto another path. So agreement between two individual runs isn’t the right test. The right test is whether a thorough search recovers the same best answer. That’s also why I’d rerun the comparison after any change that alters the arithmetic, such as the full kernel version or floating-point reordering.
Even the wider search isn’t a guarantee. In some countries, running it again from the best points found so far still turned up higher summits, so I now repeat the search until it stops improving.
None of this is specific to current accounts. Models with mixtures or switching regimes often have many plausible solutions, and an optimizer’s “converged” message doesn’t tell you which one you’ve found. Faster code bought me the time to find out.
A practical guide for evaluating convergence
Making the likelihood faster gave me time to investigate the estimates more thoroughly. To make those comparisons useful, I also needed a consistent way to assess each result. The following checks apply beyond this particular model, followed by a few GAUSS-specific considerations for compiling your own code.
Preparing and checking GAUSS kernel procedures
Separating the numerical calculations from data loading, file output, and other setup made the kernel conversion easier. These are the main considerations from the prerelease build I tested and subsequent guidance from Aptech.
Official kernel procedure documentation was still forthcoming when I tested the code. Check the release documentation for the supported inputs, functions, and compilation behavior in your installed version.
Conclusion
I started this project wanting to check whether my estimates depended on their starting values. Faster likelihood evaluation made that investigation practical: I could try more starts, compare the solutions, and examine results that an optimizer’s convergence report alone wouldn’t have led me to question. This connects with a question I explored in an earlier post on evaluating forecasts: what evidence do we need before relying on a model’s results?
Vectorization provided the largest initial gain. The new GAUSS kernel procedures reduced computation time further, first by compiling the filtering loop and then by compiling the whole likelihood. Together, those changes made a more thorough estimation workflow possible.
Time is money, especially when a calculation is repeated across countries, starting points, and bootstrap samples. For a workload like this, the time saved makes a strong case for the GAUSS 27 upgrade paying for itself. The value extends beyond finishing sooner: it gives me more opportunity to check the results before relying on them.
