# Presentation: 159452

## The ZOOM Minimization Package

- David SachsMark FischlerFermi National Accelerator Laboratory

## Minimization Package

- Object-oriented C++ library for minimization
- Minuit capabilities, and more
- Suitable for use
- In stand-alone programs
- As part of Root applications
- In other frameworks
- The OO design allows the Minimizer to ‘play well with others’
- As central part of a fitter; provided as examples will be:
- an unbinned maximum likelihood fitter
- a binned chi-squared fitter,

**Notes:**
This package is suitable for use in your C++ application, in applications in Root or other frameworks, and of course as a component in other utilities.  An example of the latter is a fitter, and a couple of fitters will be included with the distribution because they make for fine illustrations of how to use the Minimizer.

## Benefits of OO Design

- User convenience
- Domain concepts reflected more directly in code
- Natural to encapsulate function to be minimized
- Avoidance of global variables allows
- multiple simultaneous minimizations
- use in multi-threaded applications
- Extensibility
- Separation of concepts allows for independent replacement
- Benefits from *any *revision of code:
- Find features that “should have” been there
- Sometimes you uncover flaws... yes, even in Minuit

**Notes:**
It takes a long time to decide what the natural concepts and interface for minimization are, and to read, understand, document, and re-implement the methods in a complex package.  But this forced re-inspection of methods has some advantages.  Sometimes, you find features that fit in well with the existing capabilities, and which a user might justly expect to be there, but which were overlooked in the original program.  An example of this in Minuit is half-open ranges for parameters.
Sometimes, you uncover flaws – yes, even in the well-tested Minuit.  In the Simplex method, we have stumbled on two simple coding errors which can slightly increase the number of function calls needed for a given problem, one spot where the original code is susceptible to division by zero, and an algorithmic vulnerability which can cause the method to converge to a wrong point for relatively simple problems.  
While most of these flaws must be innocuous (after all, Minuit is used successfully by thousands of physicists, albeit with a bit of “folklore” about how sometimes things go wrong but you just start over from a slightly different point and they usually work) it is better to fix them than to leave them alone.
Once you have the restructured product, all the usual advantages touted for object-oriented code apply.  A particular improvement is that since well-designed C++ code avoids global variables, it becomes possible to run multiple minimizations independently, or a minimization within a minimization.  And since the package makes it trivial to minimize a functor, a class which behaves as a function but can have internal state data and other methods, the user doesn’t need to use his own global variables to supply information to the target function.

## Brief Survey of  C++ Minimizers

- Other C++ packages are available, but do not use modern C++
- and so lack the OO advantages
- Other Projects include
- The MINUIT Project [SEAL/MathLib]
- Hand-translation of Minuit (James and Winkler)
- ROOT minimization/fitting
- Hand-smoothed automated translation of Minuit

**Notes:**
The situation with regards to commercial minimization packages is, to our knowledge, about the same as it was a year and a half ago at the last CHEP conference – none of them really exploit the capabilities of C++ to make them useful as a component in somebody’s application, most of them are inferior for HEP needs to good old Minuit, and certainly none of them are so good that they will overcome the HEP community’s distaste for paying for a minor software component.
On the front of C++ Minuit clones, however, the situation has changed quite a bit in the past year.  Root has a fitter, based on a translation of Minuit. There are other fitters around, mostly based, again, on Minuit or very similar code.  And Fred James and Matthias Winkler are well along in a project to provide a C++ version of Minuit that is sticking very close to the original interface and concepts.

## An Example of Use

**double f(const std::vector<double>& x) {**
** ****return x[0]*x[0]+ //... some function of **
**	****	****	****    ****//  x[0] thru x[4]**
**}**
**int main() {**** using namespace Minimization;**
** ****const int Ndimensions = 5;**
** ****Problem m(UserFunction(Ndimensions, f));**
** ****m.minimize();**
** ****cout << m.currentPoint() << endl;**
**}**

## Directions of Extensibility

- class Terminator
- When do we stop improving the minimum?
- class Domain
- How do we restrict the allowed values of  the parameters?
- class Algorithm
- How do we refine existing minimization methods?
- How do we add new minimization techniques?

**Notes:**
The Minimization package is built to accommodate extensions by the user community, and by algorithm experts.  Extensions are added by creating sub-classes of the classes of Terminator, Domain, and Algorithm.
The Terminator classes can represent criteria like “how small is my estimated distance to minimum,” how small is the largest remaining uncertainty in any of the parameters,” “how many times have I called the function.”  And the Minimizer will be able to combine several Terminators in a natural way. 
The Domain expresses the notion that not every combination of parameter values is permitted.  The Domain corresponding to that available in Minuit is RectilinearDomain, in which each parameter is mapped into a fixed range – in this package, the range can be open at both, neither, or either end.
And Algorithms -- the package allows for adding other algorithms, for example, that used in FUMILI.
These possibilities for adding to the package are shown here in the order of most common and requiring the least sophistication, to more rarely done by experts.  It might be typical for a user to decide on a better rule for when to terminate, and this requires little knowledge beyond that of where to look for the relevant data.  It would take a bit more attention to detail to change the mappings used in the Rectilinear domain or to implement a new type of domain.  And adding a new algorithm will require, of course, expertise in creating a the algorithm – the C++ issues arising to work that algorithm into this Minimization framework are small compared to the effort to create a superior minimization method.

## Terminators

- A *terminator* is an object that can be asked whether or not minimization is complete
- it must be an instance of a class derived from class Terminator
- Several standard *terminators* are provided
- These can be combined by ORs and ANDs
*Terminators* typically make use of the state of the minimization
- Users can easily create their own *terminator * classes

**Notes:**
A Terminator is a class derived from the base Terminator class.  This has a virtual method called “finished()” which is called at times when the algorithm decides this might be a valid stopping point, and wants the user’s opinion on the matter.  
Hopefully, most users can do what they want by using and combining in logical ways the pre-packaged Terminator classes.  But it’s not hard to create your own class derived from Terminator.  A typical example is on the next slide. This asks the question “Is the improvement in estimated distance to minimum per function call, averaged over the last N function calls, still big enough to be worthwhile continuing the minimization?” 

## An Example Terminator

**class NCalls : public Terminator {**** int n;**
**public:****  explicit NCalls(int limit):n(limit) {}****  ... **
**  ****TerminationType ****  finished (const ProblemState& p)  {****     return p.functionCallCount >= n****         ? TTStop****         : TTContinue;**
**  ****}**
**};**

**Notes:**
Here is our new Terminator class.  The familiar idiom of extension by writing a derived class is used.
The actual content is in the method finished( ), which returns one of several status codes.  Here, if we have gone a sufficient number of function calls since last applying the test, we see if the rate of improvement is small enough to decide to stop.

## Domains

- A *domain* is an object that maps *internal coordinates* to *external coordinates* and vice versa
- It must be an instance of a class derived from Domain
- Custom *domain *classes must provide
- Mapping from external *M*-space to “interior” unrestricted *N*-space, and vice versa
- Gradient calculation
- One *domain* is provided: RectilinearDomain
- independent [-, ], [a,b], [-, b], [a, ] intervals
- maps for [a,b] intervals match those in Minuit

**Notes:**
The next sort of extension flexibility is that of defining a new domain.  A domain expresses the notion that the function becomes meaningless outside some region in parameter space and must not be evaluated with  meaningless parameters.
The only domain provided as part of the package is that of the rectilinear domain, in which some of the parameters may be restricted to a range.  In this domain, the range can be of the form “no greater than K,” or, as in Minuit,
“from A to B” or unrestricted.
But let’s say you have a function of three variables, defined only on the interior of a sphere.  Or a function defined only for all parameters positive, and with a sum no greater than 1.  Here, a different Domain will be useful.

## Some Possible Custom Domains

- Orthogonal domain, but with sigmoid mappings instead of arcsine mappings
- Superior for some cases, worse for others
- Multiple probability space
- All parameters must be non-negative
- Parameters must sum to 1
- Interior (or surface) of N-Sphere
- Lots of possible mappings
- But watch out if minimum lies very near a mapping singularity

**Notes:**
One simple example of a custom domain is an orthogonal domain, but with constraints enforced by a different map than is used in Minuit, perhaps a sigmoid that gradually approaches the boundary as the interior coordinate goes to infinity.  This mapping is likely to be slightly better when an algorithm wants to take large steps which might shoot past a boundary, but  will have slower convergence if the minimum is right on a boundary.  

Another example appears when the parameters in a problem represent probabilities of various mutually exclusive possibilities.  Or there might be a problem on an N-sphere or ellipsoid.

For each of these possible cases, the best way to write the new class is to start with RectilinearDomain, substitute appropriate routines for the methods that let the user establish the limits, and replace the maps, inverse maps, gradients and inverse gradients with maps that work for the new restrictions.  The C++ coding sophistication needed is very slight; the main work goes into producing suitable differentiable maps to enforce the domain.

## Algorithms

- An *algorithm* is an object that encapsulates a minimization technique
- It must be an instance of a class derived from Algorithm
- Custom *algorithm* classes can be written by experts.
- Two algorithms are provided:
- Migrad
- Simplex
- Other algorithms will be provided:
- Seek and others completing Minuit collection
- Further algorithm implementation possibilities:
- Fumili
- biConjGradStab
- LEAMAX

**Notes:**
This motivates adding a new algorithm, a simplex algorithm but with the improved edm estimation plus some way to avoid the false convergence phenomenon in reasonable functions.   It is still going to take a bit of thought and perhaps cleverness to come up with improvements to avoid false convergence.  But the virtue of the object-oriented design is that once this mathematically challenging part has been done, the coding and integration of the new Algorithm will be straightforward.

## Adding a New Algorithm

- We might wish to add a Simplex-like *algorithm* but with:
- improved estimation of the distance to the minimum (edm)
- strategies to avoid false convergence
- The mathematics of finding the best strategy is still hard
- OO design makes adding the algorithm straightforward, once the proper strategy is discovered

**Notes:**
This motivates adding a new algorithm, a simplex algorithm but with the improved edm estimation plus some way to avoid the false convergence phenomenon in reasonable functions.   It is still going to take a bit of thought and perhaps cleverness to come up with improvements to avoid false convergence.  But the virtue of the object-oriented design is that once this mathematically challenging part has been done, the coding and integration of the new Algorithm will be straightforward.

## The Simplex Algorithm

- The algorithm:
- N+1 points in N-space form a simplex
- Repeatedly reflect the “worst” point about the centroid of the remaining points
- Sometimes shrink if that improves function
- Sometimes expand if reflected point is best point
- When estimated distance to minimum is small, stop
- Problem:  The edm estimate formed by simplex “is largely fantasy.”

**Notes:**
The last part of this presentation will discuss a known headache in the simplex method and one of the flaws which was found during the rewriting and testing.  I’ll use these to motivate wanting to add a new algorithm, and in fact this algorithm makes good use of the object-oriented nature of the package, because it involves a minimization within a minimization.
In the simplex algorithm, you have N+1 linearly independent points forming a tetrahedroid simplex in N dimensions.  You repeatedly reflect the point with the worst function value, thru the hyperplane defined by the other points.  Depending on whether this improves the function and by how much, the simplex may either shrink or grow in that direction as well.
The idea is that the simplex squirms its way toward lower points in the function, taking bigger steps if the function keeps improving, and getting thinner to fit along steep valleys and continue going downhill.  Once it surrounds a local minimum the simplex rapidly shrinks inward to converge.  When the estimated distance to minimum – the difference between the best and worst point in the simplex – is small enough, the method stops.  If you were presented with a function whose values lay along a collection of planar “shards,” then migrad would do a miserable job once its step size became comparable to the shard scale, while simplex should converge quite well.
The known headache is that the quantity used by simplex as an estimated distance to minimum is  nothing of the sort.  It is useful as a convergence criterion, but if one cares about the actual distance to the function minimum, it really gives no clue until such time – if ever – that the method has converged on the true minimum.

## Discovery of a Flaw in Minuit

- Simplex can converge to a non-minimum
- Even on well-behaved (quadratic) functions!
- When starting a large distance from the true minimum
- How does this happen?
- Simplex becomes thin, aligned with an equipotential
- There is a good direction to expand off that line but...
- simplex shrinks before it finds that direction and...
- triggers convergence criterion based on the false edm
- If the *edm* were big at that point, we would at least know something was wrong
- The fix is still not trivial but knowing a meaningful edm is a start

**Notes:**
The flaw in simplex was discovered while testing the correctness of our re-write, using a 5-variable nearly quadratic function.  This is quite well behaved, so one would not expect any problems.  However, if the correct minimum is near the origin, and the starting point is very far away, the simplex method can converge onto a point which is not even a local minimum! 
What happens is that the simplex, after wandering around for a while, tends to contract into a long thin shape, lying roughly aligned with an equipotential line, far from the origin.  When this happens, although there is a direction in which the simplex could expand and yet improve on the function, most reflections in fact make the value worse.  The simplex starts to shrink in on itself, much faster than its center is moving in parameter space.  And unless you are lucky enough to hit on the golden direction before the volume becomes too small, the edm criterion terminates the search.
If the estimated distance to minimum were more meaningful, we would at least have an idea that a tiny simplex volume is inconsistent with the large edm, and might be able to take action to fix the problem.

## Apparent Convergence

- MINUIT

- Minimization

- The simplex algorithm is stopped when it indicates it is near the function minimum.

**Notes:**
Here we see the simplex method’s convergence estimate: the estimated distance to the minimum. MINUIT believes it has reached the minimum in half the time that Minimization does...

But actually, the problem is that MINUIT’s EDM is badly in error; in reality, the minimum is still far away. Minimization avoids so large a mis-estimation, and converges correctly.

Note that as the minimum is approached, the EDM eventually becomes very accurate.

## Convergence Failure

- MINUIT’s simplex indicates convergence while still far from it.
- The error is fixed in Minimization

- MINUIT

- Minimization

## Conclusion

- Minimization package available from
- _[http://](http://cepa.fnal.gov/aps/minimization.shtml)__[cepa.fnal.gov/aps/minimization.shtml](http://cepa.fnal.gov/aps/minimization.shtml)_
- Standard build mechanism (*configure/make*)
- Modern, Object Oriented C++
- Suitable for some multi-threaded operation
- Designed to be extensible to new
- Termination criteria
- Domains (parameter restrictions)
- Algorithms
- Corrects a flaw in Minuit
- In the Simplex algorithm

**Notes:**
Since we are just now beginning to let users loose on this package, the term “conclusions” is a bit of an overbid. 
But to summarize this presentation:
The object-oriented minimization package is available at this link.  It is in modern C++, and designed to be extensible in several directions.  
And in the course of doing this re-write, we have uncovered several minor flaws (and one serious problem) in some of the less-heavily used Minuit algorithms.

## The ZOOM Minimization Package

- David SachsMark FischlerFermi National Accelerator Laboratory