上回講過 Batch Gradient Descent ,今回講的是適用於 big data 的 Stochastic Gradient Descent 。
處理 Big data 的問題是,巨大的 data 無法一次過讀入記憶體,再一次過處理。此外,由於數據量龐大,如果次次都要將整份數據連算一次才算完成一次 iteration 更新 theta 的話,會好慢。
看看 batch gradient descent 每個 iteration 的算式,會發現有一個 summation 的計算。

因此,每次要 update theta ,就要先完成那次 summation 才可。對於一般大小的數據,此種計算還可以。但是如果數據很大就會很慢。
於是有人想出一招叫做 Stochastic gradient descent (SGD),是無需進行那個 summation 。每次處理一列數據,每處理一列更新 theta 一次。

通常將數據由第一列開始進行 SGD 一直至最後一列一圈就會 converge ,有時卻要行幾圈才可以。這個 SGD 有一個缺點,是未必可以達到 global minimum ,亦要小心處理 learning rate 的問題。
我拿來玩的數據,用現代的標準完全不算是甚麼 big data ,只有 25000 列而己。數據是由中文大學兒科捐給 ucla ,是該學系正在進行的大型流行病學研究的部份數據。今次想試試用兒童的高度估計重量,是完全無聊的,主要是求其有些數據試試 SGD algorithm 而已。此外,如此大小的數據,用 R 的 lm 進行 regression analysis 是全無問題的,是完全無需要用 SGD 。

[code]
## < -
## http://www.chainsawriot.com
stoGradDesc <- function(y,X, theta, alpha, rep.time=1, accumJ=100, jplot=FALSE) {
tempJ <- c()
Jhist <- c()
for (r in 1:rep.time) {
for (i in 1:length(y)) {
hx <- X[i,] %*% theta
J <- 0.5 * (hx-y[i])^2
newtheta <- theta - (alpha * t((hx-y[i]) %*% X[i,]))
theta <- newtheta
tempJ <- append(tempJ, J)
if (i %% accumJ == 0) {
aveJ <- mean(tempJ)
Jhist <- append(Jhist, aveJ)
tempJ <- c()
}
}
}
if (jplot) {
plot(Jhist, type="l", ylab="Cost", xlab="Batch")
}
return(theta)
}
cuhk <- read.table("http://dl.dropbox.com/u/5409929/weight.dat", header=TRUE)
cuhk.y <- cuhk$Weight
cuhk.x <- scale(cuhk$Height)
cuhk.X <- matrix(c(rep(1, length(cuhk.x)), cuhk.x), byrow=FALSE, nrow=length(cuhk.x), ncol=2)
ran.order <- sample.int(length(cuhk.y))
ran.y <- cuhk.y[ran.order]
ran.X <- cuhk.X[ran.order,]
initheta <- matrix(c(0,0), byrow=FALSE)
SGDtheta <- stoGradDesc(ran.y, ran.X, initheta, 0.001, 2, 100, TRUE) #converged

epsilon.converge <- function(y, X, theta, inialpha=0.0001, epsilon=0.0001, maxround=40) {
delta <- 1
round.no <- 1
while (delta > epsilon & round.no < maxround) {
SGDtheta <- stoGradDesc(ran.y, ran.X, theta, inialpha)
theta.new <- SGDtheta
delta <- max(abs(theta - theta.new))
theta <- theta.new
round.no <- round.no+1
}
return(list(SGDtheta=SGDtheta, round.no=round.no))
}
epsilon.converge(ran.y, ran.X, initheta, inialpha=0.0001)
[/code]

line 3-24:

SGD 的 implementation 。

line 25-31:

數據處理。注意,正常的 big data 理論上是不應一次過存入記憶體,而是應該一行一行的讀入(例如用 SQL query ,每次只 query 一行),避免記憶體存入大量數據。另外,第 29 至 31 行要為數據進行 random shuffle ,這是 SGD 何以稱為 stochastic 的原因。

line 32-33:

這是真正的 SGD 。結果與 lm(cuhk.y~cuhk.x) 差不多。從附送的 plot 可見, somehow 是 converge 了。

line 35-47:

這一 part 是我自行研發的。我發現在 SGD 選用甚麼的 learning rate (alpha) 及走幾多圈是很難預測的,故此用此一個叫 epsilon.converge 的東西找出來。這個東西是怎樣認定已經 converge 了呢?就是走完一圈之後 theta 的變化少於一個非常細小的數量 epsilon 。 line 33 如此靚的 alpha=0.001 和走兩圈,都是用這個 epsilon.converge 找出來的。