Problem 3: Optimization (34pts)

Write an R function returning

\[ f(x) = 1 - x + 3x^2 - x^3, \quad x \in [-0.5, 2.5]. \]

The code below defines the function in R.

f <- function(x) {
  1 - x + 3 * x^2 - x^3
}

a. Plot the function over that range and note the critical points. Check them against the truth (via calculus).

From calculus, the derivative of the function is

\[ \frac{d}{dx} f(x) = -1 + 6x - 3x^2 \]

and we can find the zeros of that function to determine the critical points. The quadratic formula says

sols <- (-6 + c(-1,1)*sqrt(6^2 - 4 *(-1)*(-3))) / (2*(-3))
sols
## [1] 1.8164966 0.1835034

The following code plots a visual of the function over the input space, agumented with vertical bars at the critical points.

x <- seq(-0.5, 2.5, length=1000)
plot(x, f(x), type="l")
abline(v=sols, col=2, lty=2)
legend("top", c("f(x)", "crits"), col=1:2, lty=1:2, bty="n")

b. Use an R library function to find those critical points numerically, and check them against the plot/truth.

The code below funds the minimum, which is in the left half of the space, by minimizing the function in an appropriate range.

optimize(f, lower=-.5, upper=1)
## $minimum
## [1] 0.1834869
## 
## $objective
## [1] 0.9113379

Next, the code below optimize (this time maximizing) the function over the right half of the input domain.

optimize(f, lower=1, upper=2.5, maximum=TRUE)
## $maximum
## [1] 1.816513
## 
## $objective
## [1] 3.088662

c. How many iterations (as counted by the number of evaluations of your R function) did it take to find each critical point?

Counting the number of iterations will require modifying the function a little bit to implement a (global) counter.

f2 <- function(x) {
  cnt <<- cnt + 1
  1 - x + 3 * x^2 - x^3
}

Now, we can initialize the counter, re-run the optimizations, and report the final count. First the minimum.

cnt <- 0
null <- optimize(f2, lower=-0.5, upper=1)
cnt
## [1] 10

Then the maximum.

cnt <- 0
optimize(f2, lower=1, upper=2.5, maximum=TRUE)
## $maximum
## [1] 1.816513
## 
## $objective
## [1] 3.088662
cnt
## [1] 10

Apparently, both take ten iterations.

Another way to do this is is by finding the roots of the derivative functions. Here is the derivative as a function in R, modified to include a counter.

df2 <- function(x) { 
    cnt <<- cnt + 1
    6*x - 1 - 3*x^2 
}

First the minimum.

cnt <- 0
uniroot(df2, interval=c(-0.5, 1))
## $root
## [1] 0.1834876
## 
## $f.root
## [1] -7.747311e-05
## 
## $iter
## [1] 6
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

Looks good. the maximum, after saving the counter.

cnt.min <- cnt
cnt <- 0
uniroot(df2, interval=c(1, 2.5))
## $root
## [1] 1.816512
## 
## $f.root
## [1] -7.747311e-05
## 
## $iter
## [1] 6
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

Also good. How does the number of iterations compare?

c(cnt.min, cnt)
## [1] 9 9