(Background info taken from my blog) In logistic regression, the hypothesis function, which models the relationshiop between the dependent variable $P(y = 1)$ and the independent variable $X$, is : \begin{align*} H_i = h(X_i) &= \frac{1}{1 + e^{-X_i \cdot \beta}} \end{align*} where $X_i$ is the $i$th row of the design matrix $\underset{m \times n}{X}$, or in matrix form: \begin{align*} H &= \frac{1}{1 + e^{-X \beta}} \end{align*}
H is a $m\times 1$ matrix. Except for $X\beta$ all operations are element-wise. The cost function $J$ is a measure of deviance of the modeled dependent variable from the observed $y$
\begin{align*} J &= (1/m)\sum_{i = 1}^m [-y_i\log H_i - (1-y_i)\log (1-H_i)] \\ &= (1/m)\sum_{i = 1}^m \left[ -y_i \log \frac{1}{1+e^{- X_i \cdot \beta}} - (1 - y_i) \log \left( 1 - \frac{1}{1+e^{- X_i \cdot \beta}} \right) \right] \\ &= (1/m)\sum_{i = 1}^m \left[ y_i \log (1+e^{-X_i \cdot \beta}) + (1 - y_i) \log \left( 1+e^{X_i \cdot \beta} \right) \right] \\ \end{align*}
\begin{align*} \frac{\partial J}{\partial \beta_j} &= \dfrac {1} {m} \sum_{i=1}^m \left[ y_{i}H_{i}e^{-X_{i}\cdot \beta }\left( -X_{ij}\right ) + \left( 1-y_{i}\right) \dfrac {1} {1+e^{X_{i}\cdot\beta }}e^{X_{i}\cdot\beta }X_{ij} \right] \\ &= \dfrac {1} {m} \sum_{i=1}^m \left[ y_{i}H_{i}e^{-X_{i}\cdot \beta }\left( -X_{ij}\right ) + \left( 1-y_{i}\right) H_i X_{ij} \right] \\ &= \sum_{i=0}^m H_{i}X_{ij}\left( -y_{i}e^{-X_{i}\cdot \beta }+1-y_{i}\right) \\ &= \sum_{i=0}^m H_{i}X_{ij}\left( 1-y_{i}\left( 1+e^{X_i\cdot\beta }\right) \right) \\ &= \sum_{i=0}^m H_{i}X_{ij}\left( 1-y_{i} / H_i\right) \\ &= \sum_{i=0}^m \left( H_{i}-y_{i}\right) X_{ij} \\ &= (H - y) \cdot X_j \\\\ \frac{\partial J}{\partial \beta} &= X^T (H-y) \end{align*}
Let $f_i' = \left( H_{i}-y_{i}\right) X_{ij}$, then according to this video:
batch gradient descent can be described as:
Until convergence:
for all $j$:
$$\theta_j := \theta_j - \alpha \sum f_i'$$
and stochastic gradient descent can be described as:
Shuffle the rows of data, and until convergence:
for all $i$ in $1\cdots m$:
for all $j$ in $0\cdots n$: \begin{align*} \theta_j := \theta_j - \alpha f_i' \end{align*}
This looks straight-forward, but when I implement stochastic gradient descent in R, it's unable to converge anywhere close to the optimum, here is the code:
logreg = function(y, x) {
alpha = 1.15
x = as.matrix(x)
x = cbind(1, x)
m = nrow(x)
m1 = sample(m)
n = ncol(x)
b = matrix(rep(1, n))
newb = b + .1
h = 1 / (1 + exp(-x %*% b))
J = -(t(y) %*% log(h) + t(1-y) %*% log(1 -h))
newJ = J+.5
while(1) {
cat("outer while...\n")
for(i in m1) {
Vi = exp(-as.numeric(x[i, ]%*%b))
Hi = 1 / (1 + Vi)
Ei = (Hi - y[i])
sDerivJ = matrix(Ei * x[i, ])
newb = b - alpha * sDerivJ
}
h = 1 / (1 + exp(-x %*% newb))
newJ = -(t(y) %*% log(h) + t(1-y) %*% log(1 -h))
if((newJ - J)/J > .15) {
alpha = alpha/2
next
}
print(b)
print(newb)
b = newb
J = newJ
if(max(abs(b - newb)) < 0.001)
{
break
}
}
b
}
nr = 5000
nc = 20
set.seed(17)
x = matrix(rnorm(nr*nc), nr)
y = matrix(sample(0:1, nr, repl=T), nr)
testglm = function() {
res = summary(glm(y~x, family=binomial))
print(res)
}
testlogreg = function() {
res = logreg(y, x)
print(res)
}
print(system.time(testlogreg()))
print(system.time(testglm()))
I am wondering what went wrong.
