COMPUTATIONAL SCIENCE

SIMULATION TRACK SOME PRACTICAL OPTIMIZATION

METHODS

Geoffrey Fox

Syracuse University

Northeast Parallel Architectures Center

111 College Place Syracuse, New York 13244-4100

(315) 443-2163 gcf@npac.syr.edu

http://www.npac.syr.edu

November 25, 1996

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Abstract of Optimization Presentation

• We review problems including TSP, linear programming, and data analysis (χ^{2}).

• We discuss relevance of central limit theorem, least squares, maximum likelihood, histograms, and scatterplots.

• We use simple length measurement, as well as full scale physics experiment to discuss χ^{2}, maximum likelihood, and their contrast.

• Computer implementations are discussed for Nonlinear Minimization, linear programming, and histogramming.

• The foils end with notes for a set of CPS713 projects in this area.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Optimization Examples—I

• Given f(\underline{x}), find minimum value of f as a function of \underline{x}—for example,

• f(\underline{x}) = x_{1}^{2} + x_{2}^{2} (continuous variables) minimum x_{1} = x_{2} = 0

Travelling Salesman Problem

• Let

i = 1 … N  label cities
η_{ip} = 1 if traveller is at city i at time step
p = 1 … N
d_{ij} = distance from city i to city j.

• Minimize

D = N ; ∑ ; i=1 N ; ∑ ; j=1 N−1 ; ∑ ; p=1 d_{ij} η_{ip} η_{jp+1}

This is a (NP Complete) Discrete Variables example. The concept of NP completeness, discussed in detail by computer science, appears of little practical significance.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Examples – II

• Given a collection of molecules, find its ground state (minimize its energy as a function of atomic positions) and see if it will make a nice drug.

• Given a set of pixels, forming an image (as seen by your eyes), optimize your life based on information—pattern recognition/(artificial) intelligence (e.g., if image contains a lion, run away; while if image contains football thrown by Jim Kelly, catch it and run to goal line).

• Given a set of ELINT (electronic intelligence = signals as a function of position), decide on best steps to win the battle—spatial reasoning or situation assessment.

• Minimize ∫ (∇ϕ)^{2} d^{3}x which is variational principle that derives Finite Element method

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Examples – III

• Given position on chess board, decide on best move. Write down a tree

and expand several “plies” (two plies—one by a white, and one by block make a move). Assign values to final position. Chosen move is one that leads to highest (best) value of worst possible scenario—minimax problem.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Examples IV – Linear Programming

Transportation

• Take M warehouses and N retail outlets. Consider one product. We wish to minimize cost of shipping products from warehouses to outlets.

• a_{i} units at warehouse i

• b_{j} units needed at outlet j

• x_{ij} units shipped from i to j at cost c_{ij} each

constraints \{ \begin{matrix}N ; ∑ ; j=1 x_{ij} & ≤ a_{i} \\ M ; ∑ ; i=1 x_{ij} & = b_{j} \end{matrix}

• Minimize N ; ∑ ; i=1 M ; ∑ ; j=1 c_{ij} x_{ij} objective function.

• This is Linear Programming—objective function and constraints are linear.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Optimization Examples – V

Interpolation

• Given a set of values y = y_{i} at x = x_{i}, find the best estimate of y = ỹ at x ≥ x̃.

• Typically do by assuming y(x) is a polynomial—maybe of degree 3 in this example.

• Note that interpolation, like many optimization problems, is mathematically undefined. One can find polynomials of degree 4 taking any given values of x_{1}, x_{2}, x_{3}, x_{4}, and x̃. We can only determine value at x̃ by \underline{assuming} polynomial of degree 3 (or other limiting assumption).

CPS615 Opt G.C.Fox, Dec 95  6  gcf@npac.syr.edu

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

A Collection of Methods for Optimization

Maximum Likelihood χ^{2} Maximum Entropy Combinational Optimization α-β Pruning Branch and Bound Branch and Cut Linear Programming Nonlinear Programming Integer Programming Heuristic Greedy Algorithm Conjugate Gradient Simulated Annealing Simulated Tempering Neural Networks Deterministic Annealing Genetic Algorithms

We will start with Maximum Likelihood and χ^{2} illustrated by data analysis problem

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Optimization Example – V: Analysis of Experimental Data

• Interpolation involves finding function values given function values at other points. This naturally extends to case where given function values have errors.

• We consider the general case where we have a set of observations x_{1}, x_{2} … x_{n} and we wish to deduce something from it.

• a) x_{i} are separate measurements of the length of a rod which differ due to uncertainities in measurement process (measurer drunk, rod quivering in earthquake). \underline{Goal}: What is best estimate for length of rod and error therein.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Analysis of Experimental Data (cont.)

• b) x_{i} are measurements of particle positions in a detector at the proposed SSC—Superconducting Supercollidor \underline{Goal}: Deduce which is the correct theory of the fundamental interactions. What particles exist?

• c) x_{i} could be number of people claiming they will vote for Clinton and i could label either geographic position or time sample taken. \underline{Goal}: Deduce who will win Presidential election and what strategy changes are appropriate to affect result.

• In each case, the x_{i} ’s are related to “goal” (which can be qualitative or quantitative) by probability distributions.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Data Analysis Examples—I

• Consider three examples defined already. Their characteristics are:

• (a) Length of Rod This is simplest case. Here, the length is well defined and uncertainties come from variations in measurement process itself.

• (b) High Energy Physics The events recorded are fundamentally probabilistic because quantum theory describes events this way, \underline{not} deterministically. There are no “hidden variables,” such as \underline{α} in example (c).

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Data Analysis Examples—II

• (c) Voting

• We can associate with the voting process, a function V

V = 0  doesn’t vote
= 1  votes for Clinton
= 2  votes for Bush
= 3  votes for Perot, etc.

• Each person is described by a set of attributes \underline{α} α_{1} = age, α_{2} = income, α_{3} = place they live, etc. with a big enough size vector \underline{α}

V ≡ V (\underline{α})

• Each person → \underline{α}_{i} a particular vector.

• If we were mathematicians, we could find V (\underline{α}) from survey (determine parameters in it) and so predict election result.

• But society is a “complex system.” Impossible to write down some function V even though actions might be described by some sort of function.

⇒ use probabilistic approach.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

The Computational and Numerical Issues

• 1. The underlying mathematics is that of random variables. The key basic theorem is the

• CENTRAL LIMIT THEOREM—this is at the level of CPS615.

• 2. There is a set of methods where we believe that we KNOW something precise about associated probability distributions.

• LEAST SQUARES METHOD

• MAXIMUM LIKELIHOOD METHOD

• 3. There are a set of techniques for data exploration where we search a large dataset for expected or unexpected correlations and clusterings which are qualitative and usually not describable by mathematical functions.

• HISTOGRAMS

• SCATTER PLOTS

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

SSC Scenario—I

• Data Analysis at the “SSC”

1. A typical problem which is rich enough to exhibit all these issues.

p \overline{p}  → N particles
→ ← N = 2 → 1000?

Particles labelled by

type,  π, k, p, \overline{p} …
momentum,  \underline{p}, E
and energy, mass^{2}  = E^{2} − \underline{p}^{2}
↑

essentially determined by type

• A typical experiment might record 10 million such events on some storage medium (e.g., ∼ 3 events per second for a year).

• The events might be “unbiased”, i.e., all possible events recorded or they may be selected by a “trigger.”

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

SSC Scenario II

• The events are specified by a set of parameters.

• N = number of particles, and for each particle, parameters defining them. This is a vector \underline{y}.

• Each event described by \underline{y}_{i} produced with probability p (\underline{y}_{i}).

• Any realistic detector is imperfect. “Image Processing” is used to find tracks of charged particles. Magnets used to bend charged particles where size of bend ⇒ sign of charge, magnitude of momentum.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Uncertainties in SSC Data Analysis

• Flaws in Image Processing cf. SDI problem of tracking a bunch of missiles “Mathematically ill-posed problem” use neural networks, etc.

• Geometric problems ⇒ In some regions of momentum \underline{p} particles are “lost” with some probability.

• Neutral particles only make tracks in dense material ⇒ harder to measure and “destructive”

• So first measure charged particles

• then measure photons (π^{0})

• then measure μ^{±}

• cannot measure neutrinos

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

SSC Data Analysis Problem Statement

• What you see = probability from Science convoluted (integrated) with probability from measurement probability from trigger.

• Probability from Science contains functional dependences which are understood together with those that are not. ⇒ this is also a “complex system”.

• The observed features =

P (\underline{y}_{i}) = ∫ known dependencies

convoluted with unknown dependencies

together with totally unknown probabilities.

• Can MODEL dependencies and probabilities (use suitably inspired = theoretically motivated “guesswork”) ⇒ PARAMETERIZE ignorance or use approaches which do NOT involve parameterization.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Maximum Likelihood Principle—Example

Rod Quivering in an Earthquake—I

• x_{1}, x_{2} … x_{N} are N measurements of the length of a rod.

• Owing to a Gaussian earthquake (or a Gaussian i.e., normal drunkenness) occurring during measurement, the x_{i} are scattered about the true value:

• You visit your favorite theorist (or do yourself after sobering up) who \underline{believes} that

• the probability distribution of x_{i} is Gaussian with mean α_{0} = true length and standard deviation σ_{0}.

• We may write down the (relative) probability of each measurement x_{i} as:

P_{i} (x_{i}) = c exp \{−\frac{1}{2σ_{0}^{2}} (x_{i} − α_{0})^{2}\}

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Rod Quivering in an Earthquake–II

• The Likelihood function for whole experiment of N events is defined as:

L (x_{1} … x_{N}) = N ; ∏ ; i=1 c exp \{−\frac{1}{2σ_{0}^{2}} (x_{i} − α_{0})^{2}\}
= c̃ exp \{−\frac{1}{2σ_{0}^{2}} _{i}^{∑} (x_{i} − α_{0})^{2}\}

• This is—up to a constant—the probability of the x_{i} \underline{given} the true length α_{0} and certain theoretical assumptions we represent as T.

L (x_{1} … x_{N}) = Pr {x_{1} … x_{n}|α_{0} , T }

in conditional probability notation.

• T is, for instance, the notion that distribution Gaussian and the value of σ_{0} taken (the latter being calculable from the strength of earthquake \underline{or} amount of beer drunk).

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Bayes Formula

• Use of Bayes Formula Now what we really want is the probability distribution of an estimate α of α_{0}. We can formally get this by use of Bayes law:

• First extend P (x_{i}) = exp \{−\frac{1}{2σ_{0}^{2}} (x_{i} − α)^{2}\} to be defined for any value of α: not necessarily α = true value α_{0}. Introduce:

• P (α|x_{i} , T) = Probability true length is α \underline{given} measurements x_{i} and predjudice T. Then Bayes law of conditional probability is

P (α|x_{i} , T) = P (x_{i}|α , T) \frac{P (α|T)}{P (x_{i}|T)}  (∗)

• This is the essence of the maximum likelihood principle. (∗) converts the likelihood L (x_{1} … x_{n}) = P (x_{i}|α , T) which is a function probability of x_{i} for fixed α , T into P (α|x_{i} , T) which is the desired probability distribution of α for given measurements x_{i} and predjudice T.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Interpretation of Bayes Relation – I

• Interpretation of equation (∗) on previous page

• (α) We can ignore P (x_{i}|T) because it is independent of α and fixed for our given x_{i} and T. It has the practical point that Likelihood is never normalized sensibly.

• (β) P (α|T) philosophically represents a priori knowledge of α before experiment (measurement of x_{i}) performed.

• This is shaky—not because we don’t really know how to define what |T (“given T ”) means—but, for instance, suppose we wish to parameterize “complete ignorance” of α before experiment performed.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Interpretation of Bayes Relation – II

• Natural is P (α|T) = constant

• But suppose we put β = α^{2} and then β has a priori probability φ(β|T)

φ(β|T)dβ = P (α|T)dα

i.e., φ(β|T) = \frac{1}{2\sqrt{β}} ≠ constant.

• So the concept of “complete ignorance” is not well defined. Jeffrey (a big bad English-Oxford-book) spends 100’s of pages trying to avoid this difficulty. His efforts are so much garbage.

• Fortunately, there is no trouble for large n because of both an obvious numerical reason (see next section (i) (d)) and a not so obvious mathematical reason.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Basic Likelihood Formula

• The Likelihood Maximizes Consider, again, the explicit form of the Likelihood of a rod quivering in the earthquake.

L (x_{1} … x_{N}) = c exp \{−\frac{1}{2σ^{2}} N ; ∑ ; i=1 (x_{i} − α)^{2}\}

• According to Bayes Principle (*), we now regard L as a function of α and so rewrite it:

L(α) = c exp \{−\frac{N}{2σ^{2}} (α − \frac{1}{N} n ; ∑ ; i=1 x_{i})^{2}
− \frac{1}{2σ^{2}} [_{i}^{∑} x_{i}^{2} − \frac{1}{N} (_{i}^{∑} x_{i})^{2}]\}

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Pictorial View of Maximum Likelihood

• The last term is a constant independent of α, and so as a function of α, L is just a Gaussian which maximizes at α = \frac{1}{N} N ; ∑ ; i=1 x_{i} and has width \frac{σ}{\sqrt{N}}.

• So as L = Pr (α|x_{i} , T) we say that maximum likelihood estimate of α is α^{∗} = \frac{1}{N} N ; ∑ ; i=1 x_{i} and its standard deviation is σ/\sqrt{N}.

• Note that whatever P (α) you put in the shape of LP (α|T) is essentially independent of P as N → ∞ and we get a Gaussian peaked at α^{∗} riding on top of constant value P (α^{∗}|T).

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

General Maximum Likelihood Method

• The above is glib in many ways—especially in its treatment of errors. We will repair this later with real live theorems. First, we state the obvious generalization.

• The general maximum likelihood method is to put P (α|T) = 1 in Bayes Principle (*) and form L(α) = ∏_{i=1}^{N} P (x_{i} , α) for any measurement x_{i} whose probability P (x_{i} , α) depends on value x_{i} and theoretical parameter(s) α.

• The maximum likelihood estimate of α is α = α^{∗} where L (α^{∗}) is a maximum. The error Δα is the standard deviation of the distribution of L as a function of α, i.e.,

L(α) = exp \{const. − \frac{1}{2Δα^{2}} (α − α^{∗})^{2}\}

• Note one usually manipulates

− ln L = − N ; ∑ ; i=1 ln P (\underline{x}_{i}, \underline{α})

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Lifetime Example for Maximum Likelihood—I

• Physics texts give the example of determining the lifetime τ_{0} of a particle from the observation of decay times t_{1} … t_{N} of N decays.

• The probability distribution is Poisson:

p(t)dt = \frac{1}{τ} exp (−\frac{t}{τ} ) dt

• Now dt is meaningless—it is independent of τ and only affects normalization of likelihood. The unknown P (α|T), P (x_{i}|T) terms in (B1) made this arbitrary anyway!

• Likelihood L(τ) ∝ \frac{1}{τ}e^{−\frac{t_{1}}{τ}} … \frac{1}{τ}e^{−\frac{t_{N}}{τ}} or L(τ) = exp [−\frac{1}{τ} N ; ∑ ; i=1 t_{i} − N ln τ] where we lazily convert proportionality into an = sign. Maximum is

\frac{dL}{dτ} = 0  or
\frac{1}{τ^{∗2}} N ; ∑ ; i=1 t_{i} − \frac{N}{τ^{∗}} = 0
\underline{or}  τ^{∗} = \frac{1}{N} N ; ∑ ; i=1 t_{i}

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Lifetime Example for Maximum Likelihood – II

• i.e., the maximum likelihood estimate is just the mean: such a simple result is only true because of particular Poisson example.

• In this special case, we can illustrate that maximum likelihood method does give results which “converge” as N → ∞.

• According to the central limit theorem, the variable y = \frac{1}{N} N ; ∑ ; i=1 t_{i} is a sum over independent random variables and has mean ⟨y⟩ = τ_{0} and standard deviation σ

σ^{2} = \frac{1}{N} \{⟨t_{i}^{2}⟩ − ⟨t_{i}⟩^{2}\}
= \frac{1}{N} \{\frac{1}{τ_{0}} ∫ t^{2}e^{−\frac{t}{τ_{0}}}dt − τ_{0}^{2}\}
or σ = \frac{τ_{0}}{\sqrt{N}}

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Multiplication of Experiments

• If we have two experiments that are independent but depend on same theoretical parameter α, then

• L_{1} = ∏_{i=1}^{N_{1}} p (x_{i} , α) is likelihood of

experiment #1

• L_{2} = ∏_{i=1}^{N_{2}} q (y_{i} , α) is likelihood of

experiment #2

• and clearly the combined likelihood is

L = L_{1} L_{2}

• i.e., one should multiply likelihoods when combining experiments. We will apply this in the next section (χ^{2}) to case when each L_{i} \underline{is} Gaussian in α.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

General Maximum Likelihood

• Note that one can prove that Maximum Likelihood is asymptotically the “best algorithm,” i.e., in limit N → ∞ (N total events, there is NO method which will give a smaller standard deviation).

• Also, Maximum Likelihood method will give one the “right answer,” i.e., Bias = Correct Estimate—Maximum Likelihood Estimate → 0 (like 1/N) as N → ∞.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

\underline{χ^{2} Method}

• This is standard and fully, correct (best) method for combining results from several independent experiments. It is in fact used (because it’s easy conceptually and for computer) even when it’s inapplicable (e.g., if there are correlations between supposed “independent” random variables or if the basic experiments have too few events to be Gaussian).

• The problem is: Given N experimental measurements x_{i} and errors σ_{i} plus estimates ξ_{i} of the x_{i} which are functions of N theoretical parameters. Then the probability of each observation is

P_{i} = \frac{1}{\sqrt{2π}σ_{i}} exp \{−\frac{(x_{i} − ξ_{i})^{2}}{2σ_{i}^{2}}\}

CPS615 Opt G.C.Fox, Dec 95  29  gcf@npac.syr.edu

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

χ^{2} Method Derived from Likelihood

• The Likelihood

L = N ; ∏ ; i=1 P_{i} ∝ exp \{− N ; ∑ ; i=1 \frac{(x_{i} − ξ_{i})^{2}}{2σ_{i}^{2}}\}

• so, the principle of maximizing the likelihood is equivalent to minimizing χ^{2} where

χ^{2} = N ; ∑ ; i=1 \frac{(x_{i} − ξ_{i})^{2}}{σ_{i}^{2}}

with respect to N theoretical parameters in ξ_{i}.

• (i) Notice that χ^{2} is the sum N ; ∑ ; i=1 v_{i}^{2} where v_{i} are independent Gaussian random variables of unit standard deviation and zero mean.

• We \underline{can} calculate χ^{2} distribution exactly, but in practice this is rarely necessary. Thus, we can easily show that—for fixed ξ_{i}—that

⟨χ^{2}⟩ = N
σ_{χ^{2}}^{2} = 2N (∗∗)

• and this is sufficient for large N as central limit theorem says χ^{2} will be Gaussian with above mean and standard deviation.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

χ^{2} Standard Deviation

(ii) Now (∗∗) is not quite correct because if we determine the theoretical parameters to minimize χ^{2} then implicitly ξ_{i} are also random variables that are functions of x_{j}.

• It is not too hard to show that in this case, one should just replace N by N − n (number of degrees of freedom) in (∗∗).

• In any case, (∗∗) is very important because it allows an easy criterion for the goodness of fit.

• Thus, if value of χ^{2} is many (say > 2) standard deviations away from its mean, then the theoretical form is in doubt.

• The usual maximum likelihood method does not have this advantage:

• Remember, L was unnormalized and so \underline{its} value had no significance for the goodness of fit. This is quite a difficult problem.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Scattering Experiment Example

Examples:

1) A physics experiment measures a scattering cross-section σ as a function of angle θ.

In the case of decay (of nucleus), we know rigorously that angular distribution σ(θ) is a polynomial in cos θ with highest degree ≤ 2J (spin of nucleus)

σ(θ) = a_{0} + a_{1} cos θ + a_{2} cos^{2} θ + a_{3} cos^{3} θ

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Individual Data Points

• There are two typical forms of data—firstly, isolated points.

• The measurements x_{i}’s are cross sections and the theoretical predictions ξ_{i} are linear in unknown parameters a_{j}.

ξ_{i} = n ; ∑ ; m=1 c_{im} a_{m}

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Histogrammed Data

• Now some experiments measure as suggested in previous picture at a set of discrete angles.

• More common is to measure a more or less complete range of θ but bin these, so that as in picture, basic measurement is number of events N_{i} in i’th “bin.”

• Here, we have nine bins in 0 ≤ θ ≤ 180.

x_{i} = N_{i} (up to constant)
error in x_{i} is \sqrt{N_{i}}

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Scattering Experiment—χ^{2} Minimization

• Under these circumstances, the χ^{2} method becomes:

minimize: N ; ∑ ; i=1 \frac{(x_{i} − ξ_{i})^{2}}{σ_{i}^{2}}  − (χ^{2})
x_{i} = const. N_{i}
σ_{i} = const. \sqrt{N_{i}}  σ_{i} = \frac{x_{i}}{\sqrt{N_{i}}}

where ith bin has N_{i} events. Note estimated error comes from measurement itself.

• Minimizing (χ^{2}) with respect to parameters a_{m} gives

N ; ∑ ; i=1 \frac{(x_{i} − ξ_{i})}{σ_{i}^{2}} \frac{∂ξ_{i}}{∂a_{m}} = 0 or
_{l}^{∑} A_{ml} a_{l}^{∗} = b_{m}

where ^{∗} denotes χ^{2} (or maximum likelihood) estimate.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

\underline{χ^{2} Formulae}

• The χ^{2} matrix

a_{ml} = _{i}^{∑} \frac{C_{im} C_{il}}{σ_{i}^{2}}

• which is a Positive Symmetric Matrix

= _{i}^{∑} \frac{(∂ξ_{i}/∂a_{m}) (∂ξ_{i}/∂a_{l})}{σ_{i}^{2}}

• The right hand side

b_{m} = _{i}^{∑} \frac{x_{i}C_{im}}{σ_{i}^{2}} = _{i}^{∑} \frac{x_{i}}{σ_{i}^{2}} \frac{∂ξ_{i}}{∂a_{m}}

• Therefore, \underline{a}^{∗} = A^{−1}\underline{b} is our best estimate for the parameters \underline{a}.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

\underline{Errors in χ^{2}}

• We can find errors in \underline{a} in either of two ways:

• 1. regard \underline{a}^{∗} as a function of x_{i} and calculate standard deviation from central limit theorem.

• 2. shape of likelihood

L = exp (−\frac{1}{2}χ^{2})

• We have − as far as shape of L is concerned. L = exp (Independent of a plus no terms linear in a − a^{∗} due to “maximum” condition −1/2 (\underline{a} − \underline{a}^{∗})^{T} A (\underline{a} − \underline{a}^{∗}) ).

• i.e., this above formula implies that A is error matrix for a.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Example II – Regression Analysis

• In typical regression analysis the ξ_{i} are also linear in the unknowns, e.g.,

Sales(t) = m ; ∑ ; i=1 x_{i}P_{i}(t)

• where x_{i} is unknown and P_{i} is known expansion function.

e.g., P_{1}(t) = 1

P_{2}(t) = (t − 1992)

P_{3}(t) = (t − 1992)^{2}, etc.

• P_{k}(t) isn’t linear in t, but Sales(t) is linear in P_{k}(t) as this is what counts!

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Example III – Nuclear Decay

• However, we consider nuclear decay

• e_{i} = number of decays observed in time interval t_{0} + (i − 1)δt to t_{0} + i δt.

ξ_{i} = c exp (−λ[t_{o} + (i − 1/2)]δt)

is a function of two unknowns, c and λ.

• ξ is linear in c but nonlinear in λ.

• Note we can analyze nuclear decays in two ways—maximum likelihood or χ^{2}. Let us compare these.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Maximum Likelihood and χ^{2} for Experimental Data Analysis

• a) Maximum Likelihood i labels events i = 1 … n_{ev}

• Minimize—ln L

L = n_{ev} ; ∏ ; i=1 \frac{1}{τ}c exp (−λt_{i})

• b) χ^{2} i labels bins (i = 1 … N_{bin}) in time

• Minimize χ^{2}

χ^{2} = N_{bin} ; ∑ ; i=1 \frac{(N_{i} − ξ_{i})^{2}}{σ_{i}^{2}}

• where σ_{i} ∼ \sqrt{N_{i}} and N_{i} is number of events in i th bin (t_{0} + (i − 1)δt ≤ t ≤ t_{0}).

• If N_{i} large, N_{i} will be Gaussian, whatever distribution of individual events.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

When Should You Use χ^{2} and When Maximum Likelihood?

• Which is best? Likelihood method is “best,” but χ^{2} method is computationally easier as N_{bin} << N_{ev}.

• Typically, as number of events in each bin must be large

N_{bin} ∼ \frac{N_{ev}}{50} .

• χ^{2} method will give an estimate with a larger standard deviation, but the difference will not be large if λδt is small, so probabilities exp(−λt) do not vary across bin. So, can choose correctly

• Small total number of events—no good method.

• Large number of events–Maximum Likelihood works and is best.

• Very large number of events—Maximum Likelihood and χ^{2} works, but χ^{2} will give a good answer with less computational complexity.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Least Squares

• Generally, we analyze

χ^{2} = N ; ∑ ; i=1 (\frac{e_{i} − ξ_{i}}{σ_{i}})^{2}
= N ; ∑ ; i=1 res_{i}^{2}  ← scaled residual

which shows why χ^{2} is also called least squares fit. As in above examples, expand ξ_{i} about \underline{x}_{0}

ξ (\underline{x}) = ξ (\underline{x}_{0}) + (\underline{x} − \underline{x}_{0})^{T} \frac{∂ξ}{∂\underline{x}}
+ Higher order in (\underline{x} − \underline{x}_{0})

• In this approximation, we now get “standard linear formulae.”

χ^{2} (\underline{x}) = χ^{2} (\underline{x}_{0}) + \frac{1}{2} (\underline{x} − \underline{x}_{0})^{T} A (\underline{x} − \underline{x}_{0})
− (\underline{x} − \underline{x}_{0})^{T} \underline{b} + corrections

• Minimize χ^{2} ignoring corrections gives

\underline{\underline{A}} (\underline{x} − \underline{x}_{0}) = \underline{b}

which is same equation we found for iterative solution of P.D.E. BUT only exact when ξ_{i} LINEAR in unknown parameters \underline{x}.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Nonlinear χ^{2} Minimization

• \underline{\underline{A}} is symmetric positive definitive

ξ_{i} = ξ_{i} (\underline{x}_{0}) + m ; ∑ ; d=1 t_{id} (x − x_{0})_{d}
where t_{id} = \frac{∂ξ_{i}}{∂x_{d}}
A_{dd′} = _{i}^{∑} t_{id} t_{id′}

• Basic idea:

• Think of initial \underline{x}_{0} and establish recursive

method

\underline{\underline{A}} (\underline{x}_{k+1} − \underline{x}_{k}) = \underline{b}

• Calculate \underline{\underline{A}} and \underline{b} using \underline{x} = \underline{x}_{k}.

\underline{x}_{k+1} = \underline{x}_{k} + A^{−1}b (NL.1)

• As A positive definite, inverse always exists.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Validity of Linearization When Does It Work?

• Why doesn’t this work?

• λ_{1} > … > λ_{n}.

• My experience at n ∼ 30, λ_{n} ∼ 10^{−8}λ_{1}. But we made approximation and so λ_{n} term swamped by corrections neglected in expanding ξ_{i} linearly. Equation (NL.2) gives big changes

δx̃_{m} = λ_{m}^{−1}b_{m} (NL.2)

Better is to use

δx̃_{m} = 0 (NL.3)

And, use of (NL.2) usually gives χ^{2} (\underline{x}_{k+1}) >> χ^{2} (\underline{k}_{n}). Not χ^{2} (x_{k+1}) < χ^{2} (\underline{x}_{k})

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Nonlinear Minimization That Works!

• Some solutions to this problem are

• 1) Conjugate Gradient, but note matrix no longer sparse and so iterative approach not so natural.

• 2) Marquardt’s method changes A in ad-hoc fashion, solving (NL.1) with the replaced value

\underline{\underline{A}} → \underline{\underline{A}} + QI

where Q is Marquardt’s Parameter

λ_{i} → λ_{i} + Q

• Suppose, Q ∼ \sqrt{λ_{1}λ_{n}}. Then shifts δx̃_{d} are unchanged for large λ_{d}, small for small λ_{d}.

• Increase Q if χ^{2} increases when you use (NL.1).

• Decrease Q if χ^{2} decreases when you use (NL.1).

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Linear Programming

• Minimize z = ∑_{j=1}^{r} c_{j}x_{j} where x_{j} ≥ 0 and there are Linear Constraints

r ; ∑ ; j=1 a_{ij} x_{j} ≤, ≥, or = b_{i} i = 1 … m

• Introduce new variables to convert all constraints to equalities.

• Suppose i = 1 is ≥. Write i = 1 constraint as

x_{r+1} = r ; ∑ ; j=1 a_{ij}x_{j} − b_{1} ≥ 0

r ; ∑ ; j=1 a_{ij}x_{j} − x_{r+1} = b_{1}
a_{1r+1} = −1
c_{r+1} = 0

• So linear programming becomes:

• Minimize z = \underline{c}\underline{x}

• Subject to \underline{\underline{A}}\underline{x} = \underline{b} as well as that all components of \underline{x} are positive.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Convex Regions and Linear Programming

• Set of \underline{x} satisfying constraints are convex region.

• Convex says that if \underline{x}_{1} and \underline{x}_{2} in a region, λ\underline{x}_{1} + (1 − λ)\underline{x}_{2} in region.

• Extreme point of convex region is one where x_{extreme} ≠ λ\underline{x}_{1} + (1 − λ)\underline{x}_{2}

where λ ≠ 0, 1.

• One can show that the optimal solution (if it exists) is an extreme point of the convex region.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Matrix Formulation of Linear Programming

• Consider \underline{\underline{A}}\underline{x} = \underline{b}. \underline{\underline{A}} is m × r matrix. Usually r >> m.

• Can show that extreme points satisfy m nonzero elements in \underline{x} n − m zero elements in \underline{x}

• Choose columns so non-zero elements are first m entries in \underline{x}

\underline{\underline{A}}^{reduced}\underline{x} = \underline{b}

where \underline{\underline{A}}^{reduced} is m × m.

• Simplex method is a strategy for moving from one solution \underline{x}^{extreme} to another \underline{x}^{extreme} by changing one column of \underline{\underline{A}}^{reduced} to another that is not being used.

• There is another class of methods—associated with Karmarkar at AT and T which are so called interior point methods. That is, you do not go from vertex to vertex of simplex but rather wander around inside simplex. Both simplex and interior point methods are important.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Computational Issues in χ^{2} and Maximum Likelihood

• Take Maximum Likelihood for definiteness

• 1. Complex nonlinear (in dependence of parameters) function to maximize

• 2. Each function value

− ln L = − N ; ∑ ; i=1 ln p (\underline{x}_{i}, \underline{α})

(where \underline{x}_{i} determines events and \underline{α} are unknown parameters) is sum of N independent calculations which is natural form of parallelism

• 3. Parallelism in optimization (in space of \underline{α}) depends on particular minimization approach adopted. Also depends on complexity of calculation

ln p (\underline{x}_{i}, \underline{α})

\frac{∂}{∂\underline{α}} ln p (\underline{x}_{i}, \underline{α})

versus complexity of manipulation of −LnL (\underline{α}) as a function of \underline{α}

• 4. χ^{2} and Likelihood are similar in aspects of analysis.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Matrices in χ^{2} Maximum Likelihood

• The arithmetic in “\underline{α} world” is that of full matrices which are n × n where n is number of parameters.

• In “linear” case (ξ_{i} linear in components of \underline{α}, χ^{2} quadratic in components), minimization requires single solution of equations

\underline{\underline{M}}\underline{α} = \underline{b} (∗∗)

• In “nonlinear” case (ξ_{i} nonlinear in α), one must iteratively solve (∗∗).

• Note that one expands ξ_{i} up to linear term which gives a quadratic approximations to χ^{2}.

• Thus, need to calculate

“value,” \frac{∂}{∂α_{k}} “value” NOT \frac{∂^{2}}{∂α_{k} ∂α_{j}} “value”

• where “value”

= ln p (\underline{x}_{i}) for Likelihood method
= ξ_{i} (\underline{x}_{i}) for χ^{2}method

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Basic Software Structure for Derivative Calculation

Often, values are formed inside loops. Inside loop over events

This code can get rather tedious. ADIFOR formally differentiates “existing” FORTRAN code and uses “chain rule” to generate derivative calculations.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Computational Complexity

• Maximum Likelihood VERY VERY expensive as

# events ; ∑ ; i=1

• Remember that the number of events could be 10^{7} at “SSC”. χ^{2} is cheaper as

# experiments or bins ; ∑ ; i=1

• Usually, number of experiments (bins) measured in 10 → 10^{3} range.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

\underline{Histograms}

• High energy physicists like to find new particles.

• For example, Nobel prize was given for the discover of ψ particle which is formed from “charm” and “anticharm” quarks.

• Mass of φ = 3.1 GeV = (p_{1} + p_{2})^{2} where ψ decays into particle 1 plus particle 2 where particle i has momentum \underline{p}_{i}.

• We measure p_{1} and p_{2} as well as lots of other things and plot

• There are an accumulation of events with (p_{1} + p_{2})^{2} near 3.1 GeV.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Histogram Formalism I

• We can see with a histogram. Let M_{0}^{2} ≤ (p_{1} + p_{2})^{2} ≤ M_{1}^{2}

be a region of interest.

• We have B Bins where each has width ΔM^{2}

ΔM^{2} = \frac{(M_{1}^{2} − M_{0}^{2})}{B}

• Let l ’th bin lie between M_{0}^{2} + (l − 1) ΔM^{2} and M_{0}^{2} + lΔM^{2}. l = 1 … B.

• Define Random Variable

η_{l} (event) = 1 event lies in l th bin
= 0 event lies outside l th bin

• Let ⟨η_{l}⟩ = μ_{l} = probability that an event will fall in bin.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Histogram Formalism II

• The means and standard deviations

⟨η_{l}⟩ = (1 − μ_{l}) · O + μ_{l} · 1 = μ_{l}
⟨η_{l}^{2}⟩ = (1 − μ_{l}) · O^{2} + μ_{l} · 1^{2} = μ_{l}
⟨η_{l}^{2}⟩ − ⟨η_{l}⟩^{2} = μ_{l} − μ_{l}^{2}  (Fulle)
= μ_{l} if μ_{l} << 1  (Approxe)

• Standard deviation σ_{l} = \sqrt{μ_{l}} = \sqrt{mean_{l}}

• Define

N_{l} = # events ; ∑ ; i=1 η_{l} (event i)
= # events in lth bin

• The central limit theorem applies to N_{l} or to

f_{l} = \frac{N_{l}}{N_{total}} ← Total # events

• This gives ⟨f_{l}⟩ = μ_{l} and f_{l} is Gaussianly distributed with error

⟨σ_{η_{l}}⟩ = \frac{⟨σ_{ηl}⟩}{\sqrt{N_{total}}}

CPS615 Opt G.C.Fox, Dec 95  55  gcf@npac.syr.edu

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Errors of Histograms – I

• Note this implies that N_{l}:

a) is GAUSSIAN

b) has Standard deviation

= \sqrt{N_{total}}⟨σ_{η_{l}}⟩ = \sqrt{N_{total} mean (η_{l})} = \sqrt{⟨N_{l}⟩}

• This is a very very important rule in all SAMPLING EXPERIMENTS. e.g., l = 1 vote for Clinton l = 2 vote for Perot, etc.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Errors of Histograms – II

• In a poll of 1,000 people, one finds that 25% are in lth bin (i.e., will vote for xyz).

N_{l} = 250
\sqrt{N_{l}} ≈ 16

• Thus, error in estimate of fraction voting for xyz is \underline{16%}—this is STATISTICAL error only.

• The systematic bias (is “1,000” chosen people representative) is a harder problem!

• As 0.25 = μ_{l} not that small, should reduces statistical error by

\sqrt{\frac{0.25 − 0.25^{2}}{0.25}} = \sqrt{\frac{3}{4}} = 0.86

using equations (fulle), not (approxe).

CPS615 Opt G.C.Fox, Dec 95  57  gcf@npac.syr.edu

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Histograms as Robust Data Exploration

• Histograms are unbiased ways of analyzing data. One does NOT need to know dependence of data on (p_{1} + p_{2})^{2}.

• General problem is data exploration. We have a set of events in a space of many dimensions

• High Energy Physics has momenta p_{1} p_{2} p_{3} p_{4} p_{5} … and

Particle Type 1 2 3 4 5 etc.

Political Polls

• (Unknown) parameters affecting a person’s choice for president.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

General Histogram

1. \underline{Select} subspace e.g., all events with > 30 particles or all events with a lepton. e.g., all people in phone book whose surname starts with Kr.

2. Choose some d functions f (event).

i = 1 … d
d = 1
f_{1} = (p_{1} + p_{2})^{2}
or f_{1} = choice for President

3. Set up set of bins—conventionally uniform in possible values for each f_{i}.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Scatterplots

• Histogram can be used if d = 1.

• If d = 2, we find a scatterplot, such as

• where we set up bins in two-dimensional space.

• e.g., f_{1} = Religion and f_{2} = # of z’s in name would be correlated (Catholic Church strong in Poland?)

• e.g., f_{2} = Maximum transverse momentum and f_{1} = (p_{1} + p_{2})^{2} would be correlated.

• as ψ’s tend to be produced in high transverse momentum (violent) collisions.

• d ≥ 3. The future with virtual reality exploration.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Possible Projects I—Maximum Likelihood

I) Produce a general χ^{2} or maximum likelihood data analysis system on a parallel machine in either explicit message passing or HPF. Apply to a simple example. You can use Marquardt’s method for minimization (see gcf records (my program “Manxcat” performed successfully 10^{6}χ^{2} fits), Numerical Recipes).

• Basic Structure:

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Possible Projects II

Separately or as follow-on to I.

II) Perform fit to specific dataset. Simplest: Generate “fake” data, e.g.,

f(x, y) = A exp [−\frac{(x − α)^{2}}{2γ^{2}} − \frac{(y − β)^{2}}{2δ^{2}}]

• Fixed Parameters are A, α, β, γ, δ

• events—choose x_{i}, y_{i} randomly.

P (\underline{x}) ∝ F (x, y)

• χ^{2} — divide two-dimensional plane into subsquares with (see figure) N_{x}N_{y} bins. Label bins with i.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Possible Projects II (cont.)

• Number of events in each bin is

N_{expected} ∝ _{∫bin}f(x, y)dxdy
∼ f (x_{central}, y_{central}) dxdy
with error \sqrt{N_{expected}}

• Assume Total number of events large.

• Generate “fake” experimental measurement as what you get from Gaussian distribution with

mean N_{expected}

error \sqrt{N_{expected}}

using chosen values A, α, β, γ, δ.

• Fit to general form and see that you find “chosen” values for these five parameters.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

More Possible Maximum Likelihood Projects

III) Find an interesting data set — High Energy Physics

— Stock Option Pricing tapes from NPAC Simplest: each event/χ^{2} bin on a separate node of a MIMD parallel machine. Complex: allow parallelism in calculation of USER1, as in parallel stock pricing model

IV) Compare χ^{2} with maximum likelihood

Performance vs. Accuracy Tradeoff V) Integrate into AVS VI) Evaluate High Performance Fortran

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

Projects—Data Exploration

I) Produce a generic histogram/scatterplot/other “robust” data analysis methods package on a parallel machine.

• This should be structured

• A Parallel Histogram is a simple database problem. Replicate in each node and combine at end.

• Scatterplots may be too big?

• Implement parallelism for histograms and parallelism for event generation.

COMPUTATIONAL SCIENCE SIMULATION TRACK SOME PRACTICAL OPTIMIZATION METHODS

More Data Exploration Projects

II) Test approach with some examples

• “Fake data” as with f(x, y) before

• Real data as before

• Theoretical calculation as in “Quantum Monte Carlo” (Cornell expertise) in Chemistry and Physics.

III) Evaluate High Performance Fortran

IV) Integrate into AVS

V) Integrate with visual data exploration—user is walking through a multi-dimensional space looking for expected or unexpected structure.

• Histograms/Scatterplots are one or two dimensions. There is a natural 3D extension which deserves consideration.

• What world to be explored has much higher dimensionality?

VI) Use Mathematica as front end to define selection of subspaces.