Homework "Minimization"
Implement the Newton minimization method with numerical gradient, numerical Hessian matrix, and back-tracking linesearch.
Something like this,
def newton(φ,x,acc=1e-3) :
while True : # Newton iterations
g = gradient(φ,x)
if g.norm() < acc : break # job done
H = hessian(φ,x)
qrH = QRdecomposition(H)
dx = qrH.solve(-g)
λ = 1
while λ > 1/128 : # backtracking linesearch
if φ(x+λ*dx) < φ(x) : break # good step
λ /= 2
x=x+λ*dx
return x
You might want to limit the number of Newton's steps to, say, 1000 to stop
the infinite loop in case the convergence criterion cannot be reached.
The numerical gradient of the objective function φ(x) can be calculated like this,
def gradient(φ,x) :
φx = φ(x)
gφ = vector(len(x))
for i in range(len(x)) :
dxi = abs(x[i])*2**(-26)
x[i]+=dxi
gφ[i]=(φ(x)-φx)/dxi
x[i]-=dxi
return gφ
The Hessian matrix of the objective function φ(x) can be calculated like this,
def hessian(φ,x) :
H = matrix(len(x),len(x))
gφx = gradient(φ,x)
for j in range(len(x)) :
dxj=abs(x[j])*2**(-13)
x[j]+=dxj
dgφ=gradient(φ,x)-gφx
for i in range(len(x)) : H[i,j]=dgφ[i]/dxj
x[j]-=dxj
return H
There are several unnecessary evaluations of φ in these examples, you might want to optimize them away;
Find a minimum of the Rosenbrock's valley function,
f(x,y)=(1-x)2+100(y-x2)2.
Find a minimum of the Himmelblau's function,
f(x,y)=(x2+y-11)2+(x+y2-7)2.
Record the number of steps it takes for the algorithm to reach the minimum.
In 2012 CERN announced the discovery of a previously unknown boson with mass 125.3±0.6 GeV/c² [Wikipedia].
In one of the experiments a certain cross-section involving Higgs (something like H→γγ) was measured and the corresponding data (with background subtracted – that's why some data are negative) are given as [Physics Letters B Volume 726, pages 88-119 (2013)] [arxiv:1207.7235]
# energy E[GeV], signal σ(E) [certain units], experimental uncertainty Δσ [same units] 101 -0.25 2.0 103 -0.30 2.0 105 -0.15 1.9 107 -1.71 1.9 109 0.81 1.9 111 0.65 1.9 113 -0.91 1.9 115 0.91 1.9 117 0.96 1.6 119 -2.52 1.6 121 -1.01 1.6 123 2.01 1.6 125 4.83 1.6 127 4.58 1.6 129 1.26 1.3 131 1.01 1.3 133 -1.26 1.3 135 0.45 1.3 137 0.15 1.3 139 -0.91 1.3 141 -0.81 1.1 143 -1.41 1.1 145 1.36 1.1 147 0.50 1.1 149 -0.45 1.1 151 1.61 1.1 153 -2.21 1.1 155 -1.86 0.9 157 1.76 0.9 159 -0.50 0.9
Fit the Breit-Wigner function (where A is the scale-factor, m is the mass and Γ is the widths of the resonance),
F(E|m,Γ,A) = A/[(E-m)²+Γ²/4] ,to the data and determine the mass and the experimetnal width of the Higgs boson (the experimental width is only the upper limit on the physical width).
You should fits the data by minimizing (using your own minimizer) the deviation function
D(m,Γ,A)=Σi[(F(Ei|m,Γ,A)-σi)/Δσi]2 .
You can send the data into the standard input stream of your program like this,
mono main.exe < higgs.data.txt 1> out.txt 2> logand you can read the data table into generic lists (either our genlist or the System.Collections.Generic.List) from the standard input like this,
var energy = new genlist<double>();
var signal = new genlist<double>();
var error = new genlist<double>();
var separators = new char[] {' ','\t'};
var options = StringSplitOptions.RemoveEmptyEntries;
do{
string line=Console.In.ReadLine();
if(line==null)break;
string[] words=line.Split(separators,options);
energy.add(double.Parse(words[0]));
signal.add(double.Parse(words[1]));
error .add(double.Parse(words[2]));
}while(true);
Plot your fit together with the experimental data.
In part A, when estimating derivatives numerically, we used (for simplicity of programming) the forward difference approximations for the derivatives,
f'(x) ≈ [f(x+δx)-f(x)]/δx , f''(x) ≈ [f'(x+δx)-f'(x)]/δx .
f'(x) ≈ [f(x+δx)-f(x-δx)]/[2δx] , f''(x) ≈ [f(x+δx)-2f(x)+f(x-δx)]/δx² .