Minimization

Minimization (D=1)

Here is the golden ratio minimizer, Implementation 12.1 - it takes a function, f, to minimize, the initial bracketing Xbrack (three values, with F(xl) > F(xm) and F(xr) > F(xm)), and the tolerance eps (describing the minimum width of the bracket):

GSMin[f_,Xbrack_,eps_] := Module[{retval,x,xl,xm,xr},
xl = Xbrack[[1]];
xm = Xbrack[[2]];
xr = Xbrack[[3]];
If[xm - xl < eps, 
Return[xm];
];
If[(xr - xm) > (xm - xl),
x =N[ xm+(1/2)(-1 + Sqrt[5])   (xm- xl)];
If[f[x] > f[xm],
retval= GSMin[f,{xl,xm,x},eps];
,
retval = GSMin[f,{xm, x,xr},eps];
];
,
x = N[xm - (1/2)(3 - Sqrt[5]) (xm-xl)];
If[f[x] < f[xm],
retval = GSMin[f,{xl,x,xm},eps];
,
retval = GSMin[f,{x,xm,xr},eps];
];
];
Return[retval];
]

and the initial bracketing routine, Implementation 12.2:

PreGSMin[f_,xl_,step_] := Module[{xr,xm,done},
done = False;
xm = xl + step;
While[done == False,
xr = N[(1/2)(3 + Sqrt[5])] (xm-xl)+ xl;
If[f[xl] > f[xm] && f[xr] > f[xm], 
done = True;
,
xm = xm + step;
];
];
Return[{xl,xm,xr}];
]

We’ll test on a function with a known minimum:

f[x_] := (x - 2)(x-3)

Plot[f[x],{x,0,4}]

First, get the initial bracketing:

X0 = PreGSMin[f,0.2,.1]

and then minimize, using a tolerance of eps = 10^(-5):

x = GSMin[f,X0,10^(-5)]

Variational ground state of hydrogen

Using our minimizer, we can compute a bound on the ground state energy of hydrogen - we start with a function of the appropriate form, and numerically integrate it using Simpson’s method from Chapter 6 (Implementation 6.1):

Simpson[f_,a_,b_,NN_] := Module[{Ia,j,dx},
dx = (b-a)/NN;
Ia = 0.0;
For[j = 0, j ≤ NN-2,j=j+2,
Ia = Ia + (1/3) dx(f[a + j dx] + 4 f[a + (j+1)dx] + f[a + (j+2) dx]);
];
Return[Ia];
]

We define our one-parameter wavefunction and associated energy integrand  from (12.25) - for limits, infinity is moved in to the dimensionless spatial value 20 (>>1), and we use 1000 steps for the Simpson integration:

Eofb[b_] := Module[{igrand,retval},
igrand[q_] := -4 (b^3) (E^(-2 b q)) q (2-2 b+(b^2) q);
retval = Simpson[igrand,0.0,20.0,1000];
Return[retval];
]

A plot of Eofb, for reference:

Plot[Eofb[b],{b,0,2}]

here’s the minimizing value of b, computed using the golden ratio minimizer (with initial bracketing performed automatically):

minval = GSMin[Eofb,PreGSMin[Eofb,0.0,.01],.0001]

giving a ground state energy of

Eofb[minval]

which is correct in these dimensionless quantities (in which the bound-state spectrum of hydrogen is -(1/(n^2))).

Levenberg-Marquardt

This is an implementation of Algorithm 12.3 - we must provide the data, d, the function to minimize, u (the norm, in this case), and the assumed functional form F, together with its derivative in terms of the parameters x, called J.  The starting values for the parameters are given in x0.

LM[d_, u_,F_, J_, x0_, lambda0_,lambdax_,eps_] := Module[{lambda, y, z, w,Jt,B,A,RHS,dx,Fc},
y = x0;
z =  Table[0.0,{i,1,Length[y]}];
lambda = lambda0;
While[Norm[y - z] ≥ eps,
Jt = Transpose[J[y]];
B =Jt.Transpose[Jt];
Fc = F[y];
RHS = Jt.(d - Fc);
w = y;
While[u[F,d,y] ≤ u[F,d,w],
A = B + lambda  DiagonalMatrix[Table[B[[i,i]],{i,1,Length[B]}]];
dx = LinearSolve[A,RHS];
w = y + dx;
If[u[F,d,w] < u[F,d,y], 
lambda = lambda/lambdax;
,
lambda = lambda lambdax;
];
];
z=y;
y=w;
];
Return[y];
]

Example 1

Below, we define the data to be a half-Gaussian, centered at zero.  We’ll use 10 data points for the fit:

NN = 10;
data = Table[2 Exp[-4 t^2],{t,0,1.0,1.0/(NN-1)}];

Our assumed functional form is A (e^(-B (t^2))) where A = (Subscript[x, 1]) and B = (Subscript[x, 2]) are the parameters we want to find

F[x_] := Table[x[[1]] Exp[-x[[2]] t^2],{t,0,1.0,1.0/(NN-1)}]

the derivative of F with respect to these two parameters defines J:

J[x_] := Table[{Exp[-x[[2]] t^2], -t^2 x[[1]]Exp[-x[[2]] t^2]},{t,0,1.0,1.0/(NN-1)}]

Finally, given a particular set of parameters x, we want to minimize u defined by:

u[F_,d_,x_] := (d - F[x]).(d - F[x])

starting off far from the solution, LM converges to the correct values

x0  = {-5,29};
{x1, x2} = LM[data, u, F,J,x0,.001,10,.01]

Example 2

This oscillating, decaying, exponential example is described in Section 12.4.

NN = 10;
data = Table[8 Exp[-Sqrt[2] t] Sin[2 Pi t],{t,0,1.0,1.0/(NN-1)}];

A plot of the 10 points of the data we’ll use to determine the fit parameters in our model:

G0 = ListPlot[data]

F and J are the functional form of the model and its derivative (with respect to the three parameters we set up)

F[x_] := Table[x[[1]] Exp[x[[2]] t] Sin[x[[3]] t],{t,0,1.0,1.0/(NN-1)}]

J[x_] := Table[{ Exp[x[[2]] t] Sin[x[[3]] t], x[[1]] Exp[x[[2]] t] t Sin[x[[3]] t], x[[1]] Exp[x[[2]] t] t Cos[x[[3]] t]},{t,0,1.0,1.0/(NN-1)}]

again, we are minimizing the length (squared) of the residual:

u[F_,d_,x_] := (d - F[x]).(d - F[x])

x0  = {1,1,1};
{x1, x2,x3} = LM[data, u, F,J,x0,.001,10,.01]

we recover good approximations to the correct coefficients (namely 8, -Sqrt[2], and 2 Pi).  The resulting curve matches the data nicely:

data = Table[{t,8 Exp[-Sqrt[2] t] Sin[2 Pi t]},{t,0,1.0,1.0/(NN-1)}];
Show[ListPlot[data],Plot[x1 Exp[x2 t] Sin[x3 t],{t,0,1}]]

Steepest Descent Line Minimization - Implementation 12.3

This version of steepest descent uses golden ratio minimization to determine the optimal step size in the one-dimensional gradient direction.

Descent[u_, du_, x_,minstep_,eps_] := Module[{y, z,F,etamin,grad,bracket},
y = x;
z = Table[0.0,{i,1,Length[y]}];
While[Norm[y - z] ≥ eps, 
grad = du[y];
 F[eta_] := u[y - eta grad];
bracket = PreGSMin[F,0.0,minstep];
etamin = GSMin[F,bracket, eps];
z = y;
y = y -etamin grad; 
];
Return[y];
]

Example 1

We need to provide the function u, and its gradient, du:

u[X_] :=(X[[1]] - 2)^2 + (X[[2]]-3)^2

Plot3D[u[{x,y}],{x,0,4},{y,0,4}]

du[X_] :=N[ {2 (-2+X[[1]]),2 (-3+X[[2]])}]

then descent finds a local minimum, given some minimum step size (for cutoff) and a tolerance eps:

Descent[u,du,{1,1},.001,10^(-6)]

Example 2

u[X_] := (X[[1]]-2)^2 (X[[2]] - 3)^2  + X[[1]]^2

Plot3D[u[{x,y}],{x,-10,10},{y,-10,10}, PlotRange->All]

du[X_] := {2 X[[1]]+2 (-2+X[[1]]) (((-3+X[[2]]))^2),2 (((-2+X[[1]]))^2) (-3+X[[2]])}

Descent[u,du,{1,1},.001,10^(-6)]

Action Example - Figure 12.2

We have a harmonic potential,

V[x_] := (1/2) 2 x^2

leading to the usual linear force,

F[x_] := -2 x

Here is the approximate action associated with the potential V, from (12.49)

S[m_,V_, X_,dt_] := Module[{vecs,retS},
vecs = Table[(1/2) m (( X[[j+1]] - X[[j-1]])/(2 dt) )^2 - V[X[[j]]],{j,2,Length[X]-1}];
retS = Sum[vecs[[j]],{j,1,Length[vecs]}];
Return[retS];
]

and its derivative, from (12.50)

gradS[m_,V_,F_, X_,dt_] := Module[{gradout,curS,dX, k},
gradout = Table[0.0,{i,1,Length[X]}];
For[k= 3,k ≤  Length[X]-2,k = k+1,
gradout[[k]] = m/(2 dt)^2 (2.0 X[[k]] - X[[k-2]] - X[[k+2]]) + F[ X[[k]]];
];
Return[gradout];
]

The descent method is modified slightly to accomodate the particulars of this problem:

ActionDescent[m_,V_, Fo_, x_,dt_, minstep_,eps_] := Module[{y, z,F,etamin,grad,bracket,count},
y = x;
z = Table[0.0,{i,1,Length[y]}];
count = 1;
While[Norm[y - z] ≥ eps, 
grad = gradS[m,V,Fo,y,dt];
 F[eta_] := S[m,V,y - eta grad,dt];
bracket = PreGSMin[F,0.0,minstep];
etamin = GSMin[F,bracket, eps];
z = y;
y = y -etamin grad; 
count = count+1;
If[Mod[count,100] == 0,Print[S[m,V,y,dt]]];
];
Return[y];
]

Set the physical parameters, m and T (the “final” time) - and the numerical Nn (number of steps) and dt.

m = 1.0;
T = 1.0;
Nn = 100;
dt = T/(Nn);

Our starting point is the linear function

x0 = Table[t,{t,0,T,dt}];

The (approximate) action is printed every hundred steps:

xf = ActionDescent[m,V,F,x0,dt,.0001,.0001];

S[m,V,xf,dt]

Here’s the exact solution to the problem:

xx[t_] := Csc[(Sqrt[2])] Sin[(Sqrt[2]) t]

and below, we plot the exact solution (solid black line), the approximation (dotted gray line), and linear starting point (gray line).

Ga = ListLinePlot[Table[{(j-1) dt,x0[[j]]},{j,1,Length[x0]}],PlotStyle->GrayLevel[0.5]];
G0 = Plot[xx[t],{t,0,1}, PlotStyle-> GrayLevel[0]];
G1 = ListPlot[Table[{(j-1) dt,xf[[j]]},{j,1,Length[xf]}],PlotStyle->GrayLevel[0.75]];
Show[{Ga,G1,G0}]

Quantum Mechanics and Monte Carlo

UpdateQM is an implementation of Algorithm 12.4 - it represents a single pass through the starting psiin, and updates all entries (outside of the Nzero indices) with random values drawn from a Gaussian with mean zero and standard deviation mag (an input parameter).  The user must also provide a function that calculates the (approximate) energy of the state psi, and that will be problem-specific.

UpdateQM[EQM_,psiin_,mag_, Nzero_] := Module[{psiout,psiup, index,Ecur,Eup},
psiout = psiin;
Ecur = EQM[psiout];
For[index = Nzero + 1, index ≤ Length[psiin] - Nzero, index=index+1,
psiup = psiout;
psiup[[index]] ={psiin[[index,1]],psiin[[index,2]] + RandomVariate[NormalDistribution[0,mag]]};
Eup = EQM[psiup];
If[Ecur > Eup,
psiout = psiup;
Ecur = Eup;
];
];
Return[psiout];
]

Next we define the approximate energy calculator - it is just (12.54), the approximate energy for a harmonic oscillator potential.

EQM[psibar_] := (1/Sqrt[2])Sum[(-psibar[[j,2]](psibar[[j+1,2]] - 2 psibar[[j,2]] + psibar[[j-1,2]])/(psibar[[j+1,1]] - psibar[[j,1]])^2 + (1/2) psibar[[j,1]]^2 psibar[[j,2]]^2) (psibar[[j+1,1]] - psibar[[j,1]]),{j,2,Length[psibar]-1}]/Sum[psibar[[j,2]]^2 (psibar[[j+1,1]] - psibar[[j,1]]),{j,2,Length[psibar]-1}]

We set up the parameters in the problem, the grid extends from -5 -> 5 (that’s set by qf, below), and we use a step size of dq = .2.  The first and last 10 entries of the wave function are zero (NZ  = 10) to enforce boundary conditions (that the wave function go to zero near the ends).  The initial wave function is generated randomly, and is shown below

SeedRandom[1111];
dq= .2;
qf = 5.0;
mag = 2.;
Nn = 2 qf/dq + 1;
NZ = 10;
psi = Table[If[q ≤ NZ|| q ≥ Nn - NZ,{-qf + (q-1) dq,0.0},{-qf + (q-1) dq,RandomReal[{-mag,mag}]}],{q,1,Nn}];
ListLinePlot[psi]

Below, we run through 2000 passes, printing the energy every 100 steps:

For[index = 1, index ≤ 2000, index=index+1,
psi =UpdateQM[EQM, psi,mag,NZ];
If[Mod[index,100] == 0, Print[EQM[psi]];
];
];

ListLinePlot[psi]

The harmonic oscillator ground state is supposed to be a Gaussian, and the above is a reasonably good approximation to it.