1 最佳化問題

目的：求函數最小值
函數乘上負號後求最大值變成求最小值問題
可應用在投資組合的配置或統計模型的預測誤差

2 Golden Search

2.1 演算法過程

從區間[a, b]開始，已知最小值在區間內
逐步縮小區間範圍，但最小值仍在區間內
當最區間上下界的差\(b' - a'\)小於設定的容忍度後停止
取區間的midpoint\((b' + a') /2\)為近似的最小值

近似值的最大誤差為\((b' - a') /2\)

2.2 決定區間上下界的方法

已知黃金比例為 \(\phi = (\sqrt{5} +1) /2\)
\(x_1 = b - (b - a) / \phi\)
\(x_2 = a + (b - a) / \phi\)
\(x_1 < x_2\)
若\(f(x_1) < f(x_2)\)，則\(b = x_2\) 為新的上界，反之則\(a = x_1\) 為新的下界
找出來的新的點的間距會保持黃金比例

f <- function(x) {
  abs(x - 3.5) + (x - 2)^2
}

curve(f, from = 1, to = 5)


# golden search function
g.search <- function (f, a, b, tol = 0.0000001)
{
  ratio <- 2 / (sqrt(5) + 1)
  x1 <- b - ratio * (b - a)
  x2 <- a + ratio * (b - a)
  f1 <- f(x1)
  f2 <- f(x2)
  while(abs(b - a) > tol) {
    if (f2 > f1) {
        b <- x2
        x2 <- x1
        f2 <- f1
        x1 <- b - ratio * (b - a)  #以新的上界計算x1
        f1 <- f(x1)  #以新的上界計算f(x1()
    } else {
        a <- x1
        x1 <- x2
        f1 <- f2
        x2 <- a + ratio * (b - a)
        f2 <- f(x2)
        }
    }
  return((a + b) / 2)
}

g.search(f, 1, 5)
## [1] 2.5

3 Newton Raphson Mathod

3.1 牛頓法使用條件

當函數的一二階導數皆為連續函數（continuous devariates）時，可用牛頓法取代Golden search

3.1.1 演算法過程

產生隨機值代入演算法
演算法將輸入值減掉一階導數和二階導數的比值，得到輸出值
將輸出值代入一階導數的泰勒近似式
重複2,3直到一階導數足夠接近0為止

3.2 牛頓法的迭代式

\(x_{n + 1} = x_{n} -\frac{f'(x_n)}{f''(x_n)}\)

3.3 一階導數的泰勒近似式

\(f'(x) \approx f'(x_0) + (x- x_0)f''(x_0)\)

3.4 迭代的停止條件

自訂tolerance \(\epsilon\)

\(|f'(x_n)| < \epsilon\)

3.5 補充

一階導數為0是極小值的必要條件，但非充分條件
二階導數大於0是極小值的充分條件
當函數有多個局部極值時，牛頓法不一定能找到最佳解
實務上若二階導數為0時，會導致無法求出泰勒近似值，解決方法是需移動一小步改變函數值

範例

f <- function(x) {
  exp(-x) + x^4
}

# 觀察函數應該在0附近有極小值
curve(f, from = 1, to = 4)


# 函數和其導數
f <- function(x) exp(-x) + x^4
fprime <- function(x) -exp(-x) + 4 * x^3
fprimeprime <- function(x) exp(-x) + 12 * x^2

# 起始值
x <- c(0.5, rep(NA, 6))
fval <- rep(NA, 7)
fprimeval <- rep(NA, 7)
fprimeprimeval <- rep(NA, 7)

# 進行6次迭代
for (i in 1:6) {
fval[i] <- f(x[i])
fprimeval[i] <- fprime(x[i])
fprimeprimeval[i] <- fprimeprime(x[i])
x[i + 1] <- x[i] - fprimeval[i] / fprimeprimeval[i]
}

# 迭代約四次後收斂，此時二階導數為正，說明其為局部極小值
# 因為此函數二階導數恆正，可確認此點是全局極小值
data.frame(x, fval, fprimeval, fprimeprimeval)
##           x      fval     fprimeval fprimeprimeval
## 1 0.5000000 0.6690307 -1.065307e-01       3.606531
## 2 0.5295383 0.6675070  5.076129e-03       3.953806
## 3 0.5282544 0.6675038  9.980020e-06       3.938266
## 4 0.5282519 0.6675038  3.881429e-11       3.938235
## 5 0.5282519 0.6675038  0.000000e+00       3.938235
## 6 0.5282519 0.6675038  0.000000e+00       3.938235
## 7 0.5282519        NA            NA             NA

3.6 Newton–Raphson

待續

4 R內建最佳化函數

4.1 一維度最佳化

f <- function (x) (x - 1/3)^2
curve(f, from = 0, to = 1)


# 使用Golden search
x.min <- optimize(f, c(0, 1), tol = 0.0001)
x.min
## $minimum
## [1] 0.3333333
## 
## $objective
## [1] 0

4.2 多維度最佳化

# 預設Nelder–Mead，可設定為Newton–Raphson或其他方法
# optim()

# par: initial value
# fn: target function
# method: optimize method

# 可處理限制式為線性不等式的函數
# constrOptim()

5 Linear Programming

線性規劃，特指函數及限制式皆為線性的最佳化問題

待續

6 Reference

“A First Course in Statistical Programming with R” [W. John Braun, Duncan J. Murdoch]

數值最佳化

hcygeorge

2019年3月18日