Fitting GLM with Quasi-Newton method

Question

I'm trying to code my own quasi-Newton algorithm for fitting GLMs in R. My results do not match up with glm and I've been over my code many, many times so I'm wondering if I have the math wrong.

My understanding is that the $\beta$ vector is updated in the following way: $$ \beta^{(t+1)} = \beta^{(t)} + \left[-S'(\beta^{(t)})\right]^{-1} \times S(\beta^{(t)}) $$

where $S$ is the score function. Because I'm doing quasi-Newton I will replace the $\left[-S(\beta^{(t)})\right]^{-1}$ with an estimate denoted by $V_t$, which is a function of the score function evaluated at $\beta^{(t)}$ and $\beta^{(t-1)}$.

Therefore my actual procedure is $$ \beta^{(t+1)} = \beta^{(t)} + V_t \times S(\beta^{(t)}). $$

If $Y_1, \dots, Y_n \sim \ iid \ f(y | \theta, \phi)$ where $f$ is an exponential family, $\theta$ is the canonical parameter, $\phi$ is the dispersion parameter, and $l$ is the log-likelihood, then we have

$$ \frac{\partial l}{\partial \beta_j}(\beta_j|\vec{y}) = \sum\limits_{i=1}^n \frac{y_i - \mu_i}{a_i \nu(\mu_i)} \times \frac{d\mu_i}{d\eta_i} \times x_{i,j} $$

where $a_i = a_i(\phi) / \phi$, $\eta_i = x_{(i)}^T\beta \ $ [$x_{(i)}$ is the $i$th row of $X$], and $g(\mu_i) = \eta_i$ for link function $g$. From this we have that $\frac{d\mu_i}{d\eta_i} = \frac{1}{g'(\eta_i)}$.

source: page 26.

I wrote a function to compute a vector score according to the above formula, taking as inputs all the things typically passed to glm such as x, y, g, g.prime, g.inv, a.vec, v.func, and initial values beta1, beta2.

I iterate the following until convergence:

1. compute $S(\beta^{(t)})$ and $S(\beta^{(t-1)})$
2. update $V_t$ using those two vectors
3. compute $\beta^{(t+1)}$ according to the above formula

The problem is that my $V_t$ and $\beta$ blow up leading to divisions by 0 and lots of NaNs. What am I missing? Is there a problem with the math or do I just have a bug that I missed? Thanks very much for any comments/suggestions.

Answer 1

In honor of it having been exactly two years since I asked this, I figured I'll post an answer now.

For some differentiable function $f : \mathbb R^p \to \mathbb R$, Newton's method seeks a stationary point by iterating $$ x^{(i+1)} = x^{(i)} - H^{-1}(x^{(i)}) \nabla f(x^{(i)}) $$ which is just an application of Newton's method for root finding to the derivative of $f$.

Computing $H^{-1}$ at every point can be a major pain for any sizable $p$, so we can use a Quasi-Newton method whereby we approximate it with a matrix $V^{(i)}$ that is a function of $x^{(i)}$ and $x^{(i-1)}$.

It turns out the formula that I was using to compute my particular $V^{(i)}$ was just wrong. The approach I describe in the question is indeed correct.