HW 03

36-711 Fall 2002

Due October 11, 2002

1. Paper and pencil, and C. The program bisection.c that we discussed in class computes a root of f(x) = 0 by the bisection method, summarized below:

given a,b such that f(a)*f(b)<0 :
loop forever {
set x = (a+b)/2;
if f(x) is approx 0, or a-b is approx 0, break out of the loop;
if f(x)*f(a) > 0, set a = x; else set b = x;
}
return x.

(a) Let x_{n} be the x that is computed after n iterations of this loop. Assume f(x) is a continuous function with a unique root x^{∗} in the original interval with endpoints a and b. Show that

|x_{n+1} − x_{n}| ≤ k · |x_{n} − x_{n−1}|

for some constant 0 < k < 1 and all n sufficiently large. From this, show that x_{n} converges to X^{∗}, and moreover

|x_{n} − x^{∗}| ≤ k^{n} · (b − a)

Definition: A sequence x_{n} is said to converge at rate (at least) α to its limit x^{∗} if

|x_{n} − x^{∗}| ≤ k^{nα} · C

for some constants 0 < k < 1, α > 0, C > 0. Note that it is sufficient to show that

|x_{n+1} − x^{∗}| ≤ k · |x_{n} − x^{∗}|^{α}

or even

|x_{n+1} − x^{∗}| ≤ k · |x_{n} − x^{∗}|^{α}

for some constant k ≠ 0 and all n sufficiently large (k is not necessarily the same in each expression). If α = 1 we say the sequence converges linearly; if α = 2 we say it converges quadratically.

Part (a) shows that the bisection method converges linearly.

(b) If you implement bisection.c using single precision (float’s) in C, it should converge in roughly 24 iterations or less. Why? If you implement it in double precision (double’s), how many iterations should it converge in?

(c) Perhaps an attractive alternative to looking at the midpoint x = (a + b)/2 each time through the loop, is to draw a straight line from (a, f(a)) to (b, f(b)) and let x be the place where the straight line intersects the x-axis. Modify bisection.c to do this, and compare its performance with the orginal bisection.c (by running both programs). Does the new algorithm still converge linearly to the root x^{∗}? (You may need to make some assumptions about f′(x) to give a precise answer.)

2. C program. Recall from statistical theory that if

L(θ) = L(θ; y_{1}, …, y_{n}) = ∏_{i=1}^{n} p(y_{i}|θ)

is a likelihood for the θ and the iid data y_{1}, …, y_{n}, then a maximum likelihood estimate (MLE) of θ is a value θ̂ that maximizes L(θ).

In practice it is usually better to work with the log-likelihood

ℓ(θ) = log L(θ) = ∑_{i=1}^{n} log p(y_{i}|θ)

since this is less likely to underflow in floating point calculations^{1}, and since ℓ(θ) is maximized if and only if L(θ) is maximized.

We usually find an MLE by setting the derivative of ℓ(θ) equal to zero

ℓ′(θ) = ∑_{i=1}^{n} \frac{d}{dθ} log p(y_{i}|θ) ≡ 0

(this is the likelihood equation) and solving for θ; i.e. we want a root of ℓ′(θ) = 0.

Please do the following:

(a) Modify root-nr.c (or write your own program) to use Newton-Raphson to calculate θ̂

for the data 1.2 2.4 1.3 1.3 0.0 1.0 1.8 0.8 4.6 1.4 which is to be modeled as

iid observations from a Cauchy likelihood, i.e.,

p(y_{i}|θ) = \frac{1}{πs} · \frac{1}{1 + (\frac{y_{i}−θ}{s} )^{2}}

with s fixed in advance to be either s = 1 or s = 0.01.

(b) For each of the two values of s, use your program to to find the MLE, starting from a large variety of initial guesses. For each s, describe your results. What is going on?

[ please turn in your C program, results, and conclusions ] 3. Splus or C. This method explores the steepest ascent/descent method.

1It is possible that we would need to add a constant to ℓ(θ) to avoid further overflow problems, but that is an issue for another day.

(a) Verify, by differentiating f(\underline{x} + t · \underline{q}) with respect to t, that the derivative in the direction of \underline{q} is \underline{q}^{T} ∇f(\underline{x}), and therefore the direction of maximum increase of f(\underline{x}) is \underline{q} = ∇f(\underline{x}), the gradient of f(\underline{x}). [and the direction of maximum decrease is \underline{q} = −∇f(\underline{x})].

(b) Write a routine in Splus that implements the “method of steepest descent” (Monahan, p. 188), and applies it to the function in Example 8.4 (same page of Monahan). The basic algorithm is

given subroutines that calculate f(x) and df(x)=gradient of f,
and given starting value x0:
 loop forever {
  let x0=xnew
  let u minimize f(x0 + u*grad f(x0))
  let xnew=x0 + u*grad f(x0)
  if relative change from x0 to xnew is small, or xnew is near zero
      break out of the loop
}

Notes:

• In Splus you can use nlmin to do the minimization, for example, or code one of the univariate search methods by hand. In C you have these minimization routines already, but you would have to code evaluation of the vector functions by hand.

• If it seems too painful to code a search method, one could simply start with a stepsize of u = 1, and repeatedly half the stepsize until f(x_{0} + u · ∇f(x_{0})) < f(x_{0}), for example. If you try this, how does this compare with using the optimal u?

Verify (and plot) the zigzag behavior illustrated in the figure on p. 189 of Monahan.

(c) Refer to the modification of steepest descent by “taking a step in the Newton direction” at the bottom of p. 188. Do Ex. 8.19, p. 197 of Monahan, to obtain the u that minimizes the quadratic approximation described on p. 197.

(d) Use your result in part (b) to revise your routine from part (a). Check (and plot) the behavior of the new algorithm. Does it still converge to the right answer? Does it zigzag? Etc. Notes:

• Monahan states on p. 188 that the relevant step comes from equation (8.5.2) but in fact it comes from equation (8.5.4), so you will need a subroutine that can compute the Hessian as well.

• Again, an alternative to calculating u according to the recipe of Ex. 8.19 would be to start with u = 1 and repeately half it until f(x_{0} + uH^{−1}(\underline{x}_{0}) · ∇f(x_{0})) < f(x_{0}).