Physics 230 -- Lab 9 (Optimization)

Fall 2015 - Transtrum

Labs developed by John S. Colton, BYU Physics & Astronomy, based in part on past labs by Branton Cambell and Steve Turley

"Optimization" refers to the process of adjusting the free parameters of a complicated system until it achieves some kind of optimal property or behavior. In this lab, we will become familiar with some of Mathematica's optimization tools. We will also explore two of the most well-known optimization problems: finding the extrema of smooth single-variable functions and fitting model functions to experimental data.  A helpful resource might be the Optimization guide page, which provides quick-reference to numerous related optimization functions and tutorials.  Also, Wikipedia has some nice information on optimization algorithms.

Clear["`*" ]
SetDirectory[NotebookDirectory[]]

Finding Extrema

Mathematica provides a number of tools for finding extrema. We’ve already seen FindMinimum and FindMaximum, used to search numerically for a local minimum or maximum beginning at a user-specified starting point. To find either type of extremum you can of course also solve for where the derivative equals zero using D (shortcut: ' ) combined with FindRoot.

New tools: NMinimize (also see NMaximize) searches numerically for a global minimum. Minimize (also see Maximize) can perform an analytical search for minima, though it defaults to NMinimize if given numeric input.  Each of these commands allows you to minimize with respect to a single variable, and also with respect to several variables at once.

Basic minimization

Evaluate the cell below to define a single-variable function, multimin[x_], which has several local minima, but only one global minimum.

multimin[x_] := 0.5x^2+Sin[4*x]^2 +0.2 x
Plot[multimin[x],{x,-3,3}, PlotRange->All]

Assignment 1. Local vs. global minima

(a) Use FindMinimum to locate a local minimum of multimin  with (Subscript[x, 0])=1.0 as a starting point.  Set a variable minim equal to the x value where the minimum is located. (Hint: use [[2]] to access the “rule” for x that is returned. Then use slash-dot to turn that rule into a regular number.)
(b) Use NMinimize to find the global minimum of of multimin.

Assignment 2. Unconstrained vs. constrained minima

What you just found with NMinimize was an unconstrained minimum. Mathematica can also perform constrained minimization, where you specify additional conditions that must be true in addition to the function having  a minimum value. Consult the help file to see how to incorporate constraints.

(a) Use NMinimize to find the global minimum of of multimin, constrained by the condition x > 0.5.
(b) Repeat for the constraint x > 5. Be able to explain your results to the TA. Hint: plot the function around this value of x.

Methods

In the Integration Lab, you saw that Mathematica has many different built-in algorithms, or “methods”, for doing numerical integrations. It’s the same for finding extrema--there are many different ways you can force Mathematica to do things if it doesn’t choose a good method automatically. These are accessed via the Method parameter. In the help file on FindMinimum under “Details and Options”, it says: “Possible settings for Method include "ConjugateGradient", "PrincipalAxis", "LevenbergMarquardt", "Newton", "QuasiNewton", "InteriorPoint", and "LinearProgramming", with the default being Automatic.”  Fortunately the default setting of Automatic nearly always works well. However, the FindMinimum help file warns that you may want to try a different method in situations where the function doesn’t have a derivative everywhere (like abs(x) at x=0).

Basic Curve Fitting

Introduction to curve fitting

Probably the most important class of optimization methods for experimental physics is “curve fitting”. That’s the process of fitting a data set with a function having one or more adjustable parameters. For example, here’s some (made-up) data for the voltage on a discharging capacitor. (Download capacitor.csv from the Physics 230 website and put it in your active directory.)

capdata = Import["capacitor.csv"] ;
plot1 = ListPlot[capdata, PlotRange->All]

From theoretical considerations, we know that the voltage of a decaying capacitor should go as: V(t)=(Subscript[V, 0])(e^(-k t)).  From the graph, it looks like the initial voltage might be around 1 V and the “time constant” might be around 0.15 s (the time constant, 1/k, is the time it takes the voltage to fall to 1/e of the original value)... but are those really the “true” or “best” values that our data is telling us? 

To get the “best values”, the FindFit function is fairly easy to use.  You define a model function with some parameters, then it tells you what the best values of the parameters are, meaning the values of the parameters that get you the least deviation from your data. To use FindFit, you specify the data, a model function, the parameters in the model that you want it to solve for, and then finally the independent variable. See this example:

modeldecay= V0 Exp[-k*t]
fitresults = FindFit[capdata,modeldecay,{V0, k},t]
bestV0 = V0 /. fitresults[[1]]
bestk = k /. fitresults[[2]]
decaytime = (1/bestk)   (* compare to guess of 0.15 s *)

We can test out how good those parameters are by plotting the modeldecay function with those parameters, on the same graph as our data. Looks great!

plot2 = Plot[modeldecay/. fitresults,{t,0,1}, PlotRange->All];
Show[plot1,plot2]

Assignment 3. A basic fit

Here’s a set of data for a damped harmonic oscillator. It’s a decaying sine wave, with a phase shift. Use FindFit to fit the data to this model function: f(x)=(ae^(-bx))sin(cx+d). Find the parameters a, b, c, and d. As I did with the capacitor decay example, make a plot of the data (as points) along with the fitted function (as a curve).

xydata = {{0.0,5.71855},{0.2,7.29998},{0.4,6.37527},{0.6,2.92445},{0.8,-0.80489},{1.0,-3.97340},{1.2,-5.20427},{1.4,-4.55806},{1.6,-2.16204},{1.8,0.65275},{2.0,2.91818},{2.2,3.67343},{2.4,3.25607},{2.6,1.43566},{2.8,-0.50996},{3.0,-1.94053},{3.2,-2.80186},{3.4,-2.19941},{3.6,-1.30169},{3.8,0.38208},{4.0,1.40752},{4.2,1.90920},{4.4,1.87848},{4.6,0.72315},{4.8,-0.42480},{5.0,-1.11011},{5.2,-1.36038},{5.4,-1.08811},{5.6,-0.66390},{5.8,0.12442},{6.0,0.78261},{6.2,1.17796},{6.4,0.83673},{6.6,0.60683},{6.8,-0.25237},{7.0,-0.63443},{7.2,-0.86788},{7.4,-0.68180},{7.6,-0.23394},{7.8,0.12773},{8.0,0.29089},{8.2,0.45893},{8.4,0.17779},{8.6,0.30576},{8.8,-0.05089},{9.0,-0.15026},{9.2,-0.35048},{9.4,-0.06596},{9.6,-0.22677},{9.8,0.15325},{10.0,0.18327}};

What is a “least squares fit”?

Back to the capacitor example. So what is Mathematica actually doing when it fits the capacitor data and calculates the best values of V0 and k? It’s the same thing we discussed in the previous section: a minimization. Typically what gets minimized when doing fits is the “squared error” between the fit and the data points. Meaning, for each data point, the “squared error” is calculated: ((((Subscript[y, data])-(Subscript[y, trial function])))^2). The trial function is the model function with a given set of parameters. The parameters are adjusted, and the error is calculated (summed over all the data points), until the computer deduces the set of parameters that give rise to the least amount of error. The result is called a “least squared error fit”, or a “least squares fit” for short.

To test this out, let’s compute the squared error for the trial function corresponding to our best fit values of V0 and k. I’ll call the modeldecay function for those particular values of V0 and k, “trial1”.

xvals =(capdata //Transpose) [[1]];  (* the column of x-values (time) *)
yvals =(capdata //Transpose )[[2]];  (* the column of y-values *)

trial1[t_] = modeldecay /. {V0->bestV0,k->bestk}   (* a specific trial function *)
yvals1 = Table[trial1[t],{t,xvals}]  (* the y-values at the given x-values, for our best fit trial function *)
squaredError1 = Total[ (yvals - yvals1)^2]  (* the squared difference between the y-values of the data and the y-values of our trial *)

Assignment 4. Capacitor again: trying other V0 and k values

Create another trial function (modeldecay, with different values of V0 and k), and compute the squared error. Do this several times, so that you become completely convinced that our “bestV0” and “bestk” values really do give the least amount of squared error.

Visualization: least squares fit

Hopefully you are now convinced that no matter what choices you make for V0 and k, you can’t get a smaller squared error than for the bestV0 and bestk choices. You can visualize this by defining a “squared error” function, a function of V0 and k, as follows:

trialx[t_,V0_,k_] = modeldecay ;
yvalstrialx[V0_,k_] := Table[trialx[t,V0,k],{t,xvals}];
squaredErrorx[V0_,k_] := Total[ (yvals - yvalstrialx[V0,k])^2];
squaredErrorx[bestV0,bestk]   (* just a test *)

I tested the error function by giving it the bestV0 and bestk values, to make sure it gave the same squared error as above.

Now we can plot the squared error, as a function of V0 and k:

Plot3D[squaredErrorx[V0,k],{V0,0,1.9},{k,6,25}]

If you look closely from a few different viewpoints, you can see that there is indeed a minimum in the squared error somewhere near the middle of the graph (a bit off to one side, actually).

Assignment 5. Find the minimum of the squared error function

Use FindMinimum to locate the minimum of that squaredErrorx function. Compare your results to the “bestV0” and “bestk” values found by the FindFit function.

Different error or “cost” definitions

Although visualizing the error as a minimum of a plot in “parameter space” (i.e., a plot with parameters as coordinate axes) is much harder for us to do when there are more than two fitting parameters, it’s really no different for Mathematica. If you think of the fitting parameters as a point (or vector) in parameter space, all the fitting algorithm has to do is move the test point around until it finds the point where the error (also called the "cost") is minimized. 

Although it may seem like our FindMinimum calculation in the last assignment just repeated what FindFit did automatically, there is at least one reason you may need to do things manually with FindMinimum rather than automatically with FindFit. That is: FindFit forces you into a particular type of cost definition, namely a least-squares fit. By contrast, with a manual minimization you can define the cost any way you want.

Assignment 6. Absolute value error

Repeat the procedure you just went through in the previous section to determine the values of V0 and k that minimize the absolute value error. Define the cost to be the absolute value of the difference between your data point and the trial function point, summed over all points. (You don’t have to make a plot unless you want to.) Compare the values of V0 and k that this method yields to the least squares method values. They should be close, but not quite the same. Hint: When I did this, I found that the FindMinimum command generated an error message. Switching to a different Method allowed it to work without errors.

Linear fits

What we have been doing so far is called “non-linear curve fitting”. Although our techniques have worked just fine in all of the examples so far, non-linear fitting can sometimes be quite challening. In a previous version of this lab Dr. Campbell wrote: “Synonyms for non-linear include difficult, painful and hard to be sure.” Many times to get good results (e.g. a global minimum rather than a local one) you need to start your parameters out with very good guesses.

Linear fitting, on the other hand, is always easy, and in fact can be solved exactly through matrix techniques which we won’t go into. A linear fit is called “linear” not because the function involved is a straight line (although it can be), but rather because the trial function is a linear combination of basis functions, those basis functions having no adjustable parameters themselves. Examples of linear fit trial functions would include:

y=m x+b   (a straight line; fit parameters m and b)
y=(Subscript[a, 0])+ (Subscript[a, 1])x+ (Subscript[a, 2])(x^2)+ (Subscript[a, 3])(x^3) (a higher order polynomial; fit parameters (Subscript[a, 0]), (Subscript[a, 1]), ...)
y=(Subscript[C, 1])sin(x) + (Subscript[C, 2])(e^(0.5x))+(Subscript[C, 3])(x^14)  (a misc. collection of functions; fit parameters (Subscript[C, 1]), (Subscript[C, 2]), (Subscript[C, 3]))

Notice how in each case the fitting parameters only exist as the coefficients multiplying the “basis functions” that make up the trial function.

Because this is computationally easier than what FindFit does, there’s a special function to determine the fit parameters for these types of situations. It’s called Fit.

Let’s generate some random data like a noisy line and test it out. Here’s what the random data looks like:

xvals = Range[0,10,0.25]; 
yvals = (10-#+RandomReal[NormalDistribution[0,.7]]&)/@xvals;
xydata = Transpose[{xvals,yvals}];
plot1 =ListPlot[xydata]

Here's how to fit the data to a straight line. The basis functions are supplied in a list as the second argument of Fit. In this case, since we want to use the function y=m x+b with fitting parameters of m and b, the basis functions are the things multiplying m and b: x and 1, respectively.

linearfit=Fit[xydata,{x,1},x]

Displaying the fit:

plot2 = Plot[linearfit,{x,0,10}];
Show[plot1,plot2]

Assignment 7. Using “Fit”

Here is some made-up data. (Download linearfitfunction.csv from the Physics 230 website and put it in your active directory.)

newxydata=Import["linearfitfunction.csv"];
ListPlot[newxydata]

(a) Use Fit to fit it to this function: y=(Subscript[C, 1])sin(x) + (Subscript[C, 2])(e^(0.1x))+(Subscript[C, 3]). That is, determine the values of (Subscript[C, 1]), (Subscript[C, 2]), (Subscript[C, 3]) that give it the least squared error. Then plot the fit along with the data.
(b) Use FindFit to do the same fit, and show that the fitted parameters have the same value.

Fourier Series

Those of you who took the majors section of Phys 123 should have been introduced to Fourier series already. Basically, the idea is that any periodic function can be written as a linear combination of sines and cosines which are multiples of the fundamental frequency. The fundamental frequency is (Subscript[ω, 0])= 2π/period  for functions in time  or  (Subscript[k, 0])=2π/(spatial period)  for functions in space.

f(x)=(∑_(n=0)^∞)(Subscript[a, n ])cos(n (Subscript[k, 0])x)+(∑_(n=0)^∞)(Subscript[b, n ])sin(n (Subscript[k, 0])x)

or, written out,

f(x)=(Subscript[a, 0])+(Subscript[a, 1])cos((Subscript[k, 0])x)+(Subscript[a, 2])cos(2 (Subscript[k, 0])x)+(Subscript[a, 3])cos(3 (Subscript[k, 0])x)+...   +(Subscript[b, 1])sin((Subscript[k, 0])x)+(Subscript[b, 2])sin(2(Subscript[k, 0])x)+(Subscript[b, 3])sin(3(Subscript[k, 0])x)+...

In other words, f(x) is a linear combination of basis functions made out of sines and cosines! Moreover, if f(x) is even, only the cosine terms (plus (Subscript[a, 0])) need to be included in the basis, and if f(x) is odd, only the sine terms need to be included.

In Phys 123 you learned and/or in Phys 318 you will learn how to calculate the a and b coefficients by doing integrals. However, when viewed in terms of today’s lab, we can calculate them numerically by considering them to be fitting parameters.

Consider this “square wave” function:

y[x_]:=If[Sin[x]>0,1,-1]
Plot[y[x],{x,-10,10}]

It’s an odd function, so let’s fit it to a set of sin(n (Subscript[k, 0])x) basis functions. The period is 2π, which means that (Subscript[k, 0])=1.

Assignment 8. Fourier series representation of a square wave

(a) Create a list of {x,y} data by “sampling” the function y(x) from -10 to 10 in steps of 0.01. Use ListPlot to view your data, to make sure it looks like the y(x) plot above.
(b) Create a list of 50 basis functions:  {sin((Subscript[k, 0])x),sin(2(Subscript[k, 0])x), sin(3(Subscript[k, 0])x), ..., sin(50 (Subscript[k, 0])x)}
(c) Do a linear fit of your sampled data over the set of basis functions. Plot the resulting function.
(d) Repeat for 5 basis functions, and for 500 basis functions. Notice as you increase the number of basis functions the fit does a better and better job of reproducing the square wave.

Incidentally, the FourierSeries package has tools for computing Fourier series expansions. If/when you get more into Fourier series, you might want to check it out.

Advanced Curve Fitting

Incorporating uncertainties

Experimental data is subject to random errors, and therefore possesses statistical uncertainty. This uncertainty affects numerical fits in at least two ways.

First, it means that there is some uncertainty in the parameters that result from a given fit. Since your data was not 100% trustworthy, the resulting parameters are not 100% trustworthy. The process of computing the statistical uncertainty in the value of fitted parameters is called “regression”. Unfortunately, Fit and FindFit can’t do regression. If you want regression, LinearModelFit replaces Fit for linear problems, and NonlinearModelFit replaces FindFit for nonlinear problems.

Second, sometimes some data points might be more reliable than others. If so, you might want to give more weight to them in the fit. From an experimental perspective, we would say that the reliable points have “smaller error bars” than the less reliable points. LinearModelFit and NonlinearModelFit are also the commands you can use to specify how much weight should be given to each data point. 

Download datawitherrorbars.csv from the Physics 230 website and put it in your active directory. Then evaluate the cell below to load in a new data set that is somewhat similar to the exponentially-decaying capacitor we started out with. This time, however, there is a third column of data which contains the experimental uncertainty for each measured y value (the “error bars” for the point). The larger the value in the third column, the less reliable the data was for that point.

xysigma = Import["datawitherrorbars.csv"]
xydata = xysigma   /. {x_,y_,s_}  -> {x,y}   (* Notice this line!! Explained below. *)
ListPlot[xydata]

What to notice: I’ve used a new trick here to rearrange and generate a new list, involving pattern recognition. In this case I’ve told it to recognize the pattern of three elements in a list (I temporarily named the elements x, y, and s), and replace that pattern with just {x,y}. You could translate that command into English as “set xydata equal to xysigma, but with the triplets of x, y, and s replaced by pairs of x and y”.  You could in similar fashion rearrange the three elements into any pattern you may want. Play around with that technique a little right now to see its power.

(* Space for you to play around. *)

To plot the data and uncertainties, I’m going to use an ErrorListPlot, which is part of the ErrorBarPlots package. It takes data in an odd format, like this: { {x,y}, ErrorBar[s] }, for each data point, so I’m going to use my new trick to rearrange the list of data again. (ErrorBar[s] uses the ErrorBar function to create an “error bar object” out of an uncertainty value.)

Needs["ErrorBarPlots`"]
oddformat = xysigma /. {x_,y_,s_}  -> {{x, y},ErrorBar[s]}  ; 
errorbarplot=ErrorListPlot[oddformat, PlotRange->All, AxesOrigin->{0,0}]

Now to fit the data. I’d like you to use a decaying exponential model, but let’s add in the possibility of a constant y-offset term since it doesn’t look like the capacitor decays all the way to 0 V.

model[x_] := a*Exp[-x/b] +c

Assignment 9. Fit with weighted uncertainties

(a) At this point let me turn it over to you. We’ll need the list of uncertainty numbers all by themselves for the next step, so use my trick to transform xysigma into just a list of s values (the third column), and call the result uncertainty.

(b) Now fit the xydata using NonlinearModelFit with Weights -> 1/uncertainty^2 as an option. The result of the NonlinearModelFit is a FittedModel object; give it the name nlm (standing for “nonlinear model”).

Behind the scenes: The Weights -> 1/uncertainty^2 option has told Mathematica to use a new cost function, somewhat like the problem above where you used the absolute value instead of the squared error. In this case, the  Weights specifier tells Mathematica to use a “weighted least squares” error function, which is computed as:  cost=(Underscript[∑ , n])(((Subscript[w, n ])((Subscript[y, n])-f((Subscript[x, n]))))^2).  Here (Subscript[w, n]) is the weighting factor for the (n^th) term. The   -> 1/uncertainty^2   part, with your definition of uncertainty from part (a) means that you want (Subscript[w, n])= 1/(Subscript[s, n]^2), i.e. you want each point weighted by the inverse square of its uncertainty. Thus we increase the importance of data points with small uncertainties and decrease the importance of data points with large uncertainties. We could have weighted things differently, but 1/uncertainty^2 turns out to be exactly the right thing to do if your (Subscript[s, n]) uncertainty values are standard deviations (due to Aitken’s extension of the Gauss-Markov thorem, if you care).

(c) Assuming that the previous step was executed correctly, you can now use nlm to study many useful properties of the fit.  Evaluate the cell below and observe that not only has it properly accounted for the "error bars" at the points, it has also produced a list of uncertainties for your fitted parameters (i.e. "ParameterErrors"). (Exactly how the parameter uncertainties are estimated is beyond the scope of this course.) See the "Scope" section of the NonlinearModelFit reference page for a complete list of the properties that can be displayed for nonlinear model fits.

nlm["BestFit"]
nlm["BestFitParameters"]
nlm["ParameterErrors"]
Show[errorbarplot, Plot[nlm["BestFit"],{x,0,32}, PlotRange->All]]

Additional Assignments

Assignment 10. Fitting Variable Star Data

The file NGC6811-35.csv, available on the Physics 230 website, is data digitized from a paper by Michael B. Rose and Eric G. Hintz of BYU entitled “A search for low-amplitude variability in six open clusters using the robust median statistic”, published in The Astronomical Journal, vol. 134, pp 2067-2078, in Nov 2007. The file contains data points representing the magnitude of the star NGC 6811-35 as a function of time. The magnitude of a star is proportional to the log of its brightness. (I.e., brightness = 10^magnitude.) The first column in the file is the time in hours at which the star was observed. The second column is its magnitude.

(a) Convert the y-axis to brightness, then ListPlot the data. (The x-axis is time.)
(b) You should see a sinusoidal variation to the brightness, albeit with a lot of noise. Fit the data to a function:  y=(Subscript[y, 0])+Acos(2πt/period - ϕ), where (Subscript[y, 0]), A, period, and ϕ are fitting parameters. Use a NonlinearModelFit so that you can also obtain the uncertainties in the parameters. Plot your data and the fit on the same graph to make sure your fit is reasonable.

Note: For reasons I couldn’t determine, Mathematica is super-picky about the functional form needed to fit this particular data. Using  + ϕ  instead of   - ϕ , for example, caused the fit to fail completely. I had similar problems when I tried Sin instead of Cos (even though they should be equivalent if given an appropriate phase shift), and when I tried to use ω instead of 2π/period. (I mentioned that nonlinear fits could be tricky!) So be sure to use the fitting function exactly as specified above.

Assignment 11, Extra Credit. Hubble’s Constant

In the Data Processing lab, you produced a graph that looks like this, illustrating Hubble’s Law. The x-axis is the distance to a galaxy (in megaparsecs); the y-axis is a measure of how fast the galaxy is traveling (redshift, in km/s). 



Recreate your x-y data from that lab and do a curve fit to deduce the slope of the best-fit line:  y = mx  (no intercept). Use the LinearModelFit command for this so that you can weight the data points in your fit according to the inverse of their distance from us (because greater distances have greater uncertainties and should therefore receive a lower weight).

The slope of the line in your fit is Hubble’s constant. Compare your value and its uncertainty to those described in the Wikipedia article on this topic. Consider that 80 years ago you could have won a Nobel prize with this analysis!

Hint: your data in the last lab may have had units associated with the distances. If so, you’ll need to first strip off the units before you can use LinearModelFit. Here’s one way to do that:

AstronomicalData["Sun", "DistanceLightYears"] 
QuantityMagnitude[%]

Wrap-up

The following functions were discussed in today's lab, or are very closely related:

Finding local minima/maxima: FindMinimum, FindMaximum
Finding global minima/maxima: NMinimize, NMaximize
Finding analytical minima/maxima: Minimize, Maximize
Nonlinear curve fit: FindFit 
Linear curve fit: Fit
Linear/nonlinear curve fits with regression: LinearModelFit, NonlinearModelFit
Result of LMF and NLMD: FittedModel (data type)
Error bars:  ErrorBarPlots (package), ErrorListPlot, ErrorBar (function to create an error bar object, as needed by ErrorListPlot)
Fourier series package: FourierSeries 


Potentially useful summary/tutorial pages: 

Guide:  Optimization