Import libraries:
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(animation)
Write my own NelderMead function:
NelderMead <- function(x, fun, tol = 1e-08, alpha = 1, gamma = 2, lo = 0.5, sig = 0.5,
maxiter = 50) {
n <- dim(x)[2] # The number of dimensions.
num <- 0
df <- fun(x) %>%
cbind(x, num)
names(df)[c(1, n + 2)] = c("f", "num")
area <- tol + 1
trc <- df[FALSE, ]
while (area > tol & num < maxiter) {
df[n + 2] <- num <- num + 1
# Sort the f(x).
df <- df %>%
arrange(f)
f1 <- df[1, 1]
fn <- df[n, 1]
xw <- df[n + 1, 2:(n + 1)]
x0 <- colMeans(df[1:n, 2:(n + 1)])
trc <- rbind(trc, df)
# Compute reflected point.
xr <- x0 + alpha * (x0 - xw)
fr <- fun(xr)
if (f1 <= fr & fr < fn) {
# Keep xr.
df[n + 1, 1:(n + 1)] <- cbind(fr, xr)
} else if (fr < f1) {
# Compute reflected point.
xe <- x0 + gamma * (xr - x0)
fe <- fun(xe)
if (fe < fr) {
# Keep xe.
df[n + 1, 1:(n + 1)] <- cbind(fe, xe)
} else {
# Keep xr.
df[n + 1, 1:(n + 1)] <- cbind(fr, xr)
}
} else {
# Compute contracted point.
xc <- x0 + lo * (xw - x0)
fc <- fun(xc)
if (fc < df[n + 1, 1]) {
# Keep xc.
df[n + 1, 1:(n + 1)] <- cbind(fc, xc)
} else {
# Shrink.
x1 <- t(matrix(df[1, 2:(n + 1)]))
df[, 2:(n + 1)] <- x <- x1 + sig * (df[, 2:(n + 1)] - x1)
df[, 1] <- fun(x)
}
}
area <- cbind(df[, 2:(n + 1)], 1) %>%
as.matrix %>%
det %>%
abs/2
# print(df) cat(area, end = '\n')
}
df[n + 2] <- num + 1
trc <- rbind(trc, df)
if (num == maxiter) {
cat("\nMaximum number of iterations was reached!")
}
df0 = colMeans(df)
return(list(x = df[2:(n + 1)], f = df[1], f.x = df0, trace = trc))
}
Use a triangle with vertices (1,1),(1,2),(2,2) as the starting simplex to find the minimum of the function \(f(x)=x_1^2+x_2^2\).
f <- function(x) x[1]^2 + x[2]^2
v1 <- c(1, 1)
v2 <- c(1, 2)
v3 <- c(2, 2)
v <- rbind(v1, v2, v3) %>%
data.frame
r <- NelderMead(v, f, maxiter = 80)
Display the trace of iterations in GIF on a filled contour plot.
tr <- r$trace
n <- tr["num"] %>%
max # Iteration times.
x <- seq(-2, 2, 0.2) %>%
as.data.frame
y <- seq(-2, 2, 0.2) %>%
as.data.frame
ff <- merge(x, y, by = NULL)
ff <- mutate(ff, z = ..x^2 + ..y^2)
for (i in 1:n) {
p <- ggplot(tr[(3 * i - 2):(3 * i), ], aes(X1, X2)) + geom_contour_filled(aes(x = ..x,
y = ..y, z = z), ff) + geom_polygon(fill = NA, colour = "red")
plot(p)
}