其實我也估不到會寫 #3 。
在 Matlab/Octave ，做功課時用過一件叫 fminunc (Find minimum of unconstrained multivariable function) 的東西。上課時講是使用 advanced optimization algorithms 時用，而課上不停講的幾隻 advanced optimization algorithms 有 BFGS, L-BFGS, DFP 之類。但別問我那些時甚麼，我仍未有時間參透。我只知道 BFGS 類，被稱為 Quasi-Newton method 。它們與 Newton-Raphson 之最大分別是，它們無需計算 Hessian Matrix ，故此無需要進行 Computational Cost 高的 Matrix Inversion 。因此，對比 Newton-Raphson 的 O(n3) ，此類 algorithms 最多也只會是 O(n2) 。
在 R 又如何使用這些 algorithms 呢？答案是用一個內建的 function 叫做 optim() 。此東西亦算方便易用，起碼無用事事自己寫，不妨作推介。使用 optim() 的重點是，一定要寫一個計算 Cost function J 的 function ，非必要的一步是要寫多一個計算 Gradient 的 function 。在寫 cost 和 gradient 的 function ，第一個 parameter 要是想 optimize 的東西，也即是 θ 。詳請請看 ?optim 。
[code]
# <-
# optim() examples
y <- mtcars$am
x <- as.matrix(cbind(mtcars$mpg, mtcars$wt))
X <- cbind(rep(1, nrow(x)), x)
scx <- as.matrix(scale(cbind(mtcars$mpg, mtcars$wt)))
scX <- cbind(rep(1, nrow(scx)), scx)
initheta <- matrix(c(0,0,0), ncol=1)
optimComputeJ <- function(theta,y,X) {
htheta <- g(X %*% theta)
J <- (1/length(y))*sum((-y*log(htheta))-((1-y)*log(1-htheta)))
return(J)
}
optimGradient <- function(theta,y,X) {
htheta <- g(X %*% theta)
gradient <- 1/(length(y)) * t(X) %*% (htheta-y)
return(gradient)
}
optimDefault <- optim(initheta, fn=optimComputeJ, y=y, X=X) #default
default.time <- system.time(optim(initheta, fn=optimComputeJ, y=y, X=X))[3]
optimLBFGSB <- optim(initheta, fn=optimComputeJ,gr=optimGradient, y=y, X=X,method="L-BFGS-B", hessian=FALSE)
LBFGSB.time <- system.time(optim(initheta, fn=optimComputeJ,gr=optimGradient, y=y, X=X,method="L-BFGS-B", hessian=FALSE))[3]
optimBFGS <- optim(initheta, fn=optimComputeJ,gr=optimGradient, y=y, X=X,method="BFGS", hessian=FALSE) # fail to converge with 100 iters
scoptimBFGS <- optim(initheta, fn=optimComputeJ,gr=optimGradient, y=y, X=scX,method="BFGS", hessian=FALSE) #OK
scBFGS.time <- system.time(optim(initheta, fn=optimComputeJ,gr=optimGradient, y=y, X=scX,method="BFGS", hessian=FALSE))[3]
barplot(c(default.time, LBFGSB.time, scBFGS.time), names.arg=c("Default (Nelder-Mead)","L-BFGS-B","Scaled BFGS"), ylab="Time",xlab="Algorithm")
[/code]
用上兩文同樣的例子。
optim() 最方便是未必需要提供 gradient 的算法，系統會自行找出來。但當然，如果有另外提供的話，可以將 convergence 時間再減少。此外， optim() 也會自行計算出 convergence ，連 iterations 的次數都無需估估下。當然，如果要設計 iterations 次數也是可以的。
從實驗可見，用 default 的 Nelder-Mead 、 L-BFGS-B 都可 Converge ，但是用 BFGS 卻不能在 100 個 iterations 之內 converge ，要用 feature scaling 才可成功。 (( optim() 回報的結果是一個 list ，其中一個 element 叫做 convergence ，可看到有否 converge 。以 optimDefault 為例，就是 optimDefault$convergence，如果是 0 的話就基本 OK 。))