# Gradient Descent
In this problem, we use gradient descent to approximate the minima of the function,

\[
f(x) = -\frac{1}{\sqrt{2\pi}}\text{exp}\big(-\frac{1}{2}x^2\big) + \frac{1}{\sqrt{2\pi}}\text{exp}\big(-\frac{1}{2}(\frac{x-1}{0.5})^2\big)-\frac{1}{\sqrt{2\pi}}\text{exp}\big(-\frac{1}{2}(\frac{x-2}{0.4})^2\big)
\]

## Question 1 
Plot the function $f(x)$ for $x\in[-2,4]$.

## Answer 1
```{r}
f<-function(x){
  - dnorm(x,0,1) + 0.5*dnorm(x,1,0.5) - 0.4*dnorm(x,2,0.4)
}
plot(seq(-2,4,0.01),f(seq(-2,4,0.01)),pch=16,col="dodgerblue",type = "l",cex=0.5,xlab="x",ylab="f(x)")
```

## Question 2
We will approximate the derivative $f'(x)$ by replacing $f$ by a linear function within a small window. We choose some small value $\delta>0$ and compute:
\[
\hat{f}_\delta (x):=\frac{f(x+\delta)-f(x-\delta)}{|[x-\delta,x+\delta]|} = \frac{f(x+\delta)-f(x-\delta)}{2\delta}
\]
and use $\hat{f}_\delta (x)$ as our estimate of $f'(x)$. The function `f.prime()` we ask you to write below implements $\hat{f}_\delta (x)$.

Write a function, `f.prime(x)`, to calculate the numerical derivative given the location $x$. Take the approximation window to be $\delta=0.001$ (though it'd be nice to allow `f.prime()` to take this as an argument as well, with a default to 0.001). What is the output of `f.prime(-2)`?

## Answer 2
```{r}
f.prime<-function(x,delta=0.001){#default delta to 0.001
  res <- (f(x+delta)-f(x-delta))/(2*delta)
  return(res)
}
f.prime(-2)
```

## Question 3 
Write a function, `grad.des(x1)` to perform gradient descent from the starting point x1. Take the step sizes $\alpha_n=\frac{1}{n}$ and precision $\epsilon=0.000001$, and let the maximum iterations be 10,000. That is compute:
\[
x_{n+1}:=x_n-\frac{1}{n}\hat{f}_\delta(x_n)
\]
Your function should iterate through n until either the precision is satisfied or the max iteration is reached, whichever comes first. It should output a list, including the number of total iterations $N$, the minima $x^*=x_N$, the minimum value $f_{\min}=f(x^*)$ and the vectors of your search trajectory for both $x$ and $f(x)$ values: $(x_1,x_2,\cdots,x_N)$ and $(f(x_1),f(x_2),\cdots,f(x_N))$.

## Answer 3
```{r}
grad.des<-function(x, # starting value
                   precision=0.000001, # tolerance
                   max.iters=10000, # maximum number of iterations
                   verbose=FALSE){ # If TRUE, print outputs for each iteration
  x_1 <- x
  n <- 1
  close <- FALSE
  x_values <- f_values <- NULL
  x_values[1] <- x_1
  f_values[1] <- f(x_1)
  while(close == FALSE & n < max.iters){
    if(verbose == TRUE) cat("Iteration ", n,"\n")
    x <- x_1
    x_1 <- x - (1/n) * f.prime(x)
    step <- x_1 - x
    if(abs(step)<=precision){ close <-TRUE }
    n <- n+1
    x_values[n] <- x_1
    f_values[n] <- f(x_1)
  }
  if(close==FALSE & n==max.iters){ cat("Max iterations reached\n")}
  return(  list(
    "n" = n,
    "x_values" = x_values,
    "x_1" = x_1, #xstar minima here
    "f_values" = f_values,
    "f_x" = f(x_1) #f(xstar)
  ))
}
```

## Question 4
Start from $x_1=-2$, what is the minimum your function finds? Plot the points $(x_1, f(x_1)),(x_2, f(x_2)),\cdots,(x_N, f(x_N))$ in red with the curve of $f$ from Question 1.

## Answer 4
```{r}
myx <- grad.des(-2)
myx$x_1 # x*: minimum
plot(seq(-2,4,0.02),f(seq(-2,4,0.02)),pch=16,col="dodgerblue",type = "l",cex=0.5,xlab="x",ylab="f(x)")
points(myx$x_values,myx$f_values,pch=16,col="firebrick",type = "l",cex=0.5)
```

## Question 5
A gradient descent algorithm generally finds local minima, and these minima may not be global. One way to address this is to randomly start gradient descent at several different points and compare the function values at the resulting local minima. Generate 10,000 random starting points between $(-2,4)$ (after setting seed with `set.seed(2018)`), and calculate the gradient descent result for each of them. Plot the histogram of the 10,000 minima you found. What is the global minimum?

## Answer 5
```{r}
set.seed(2018)
start.val<-runif(10000,-2,4)
minima<-vector("list",length=length(start.val))
for(i in 1:length(start.val)){
  minima[[i]]<-grad.des(start.val[i],0.01,max.iters=1000)
}
mymin<-unlist(lapply(minima,function(x){return(x$x_1)}))
myfmin<-unlist(lapply(minima,function(x){return(x$f_x)}))
hist(mymin,xlab="Minima values",main="Histogram of 10000 minima values",col="dodgerblue")

mymin[which(myfmin==min(myfmin))]
```