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.

  • $\begingroup$ Post your code and a sample dataset.... $\endgroup$
    – Hong Ooi
    Mar 23, 2015 at 4:39
  • 1
    $\begingroup$ I deliberately didn't post my code because I didn't want to make this about code review. I just want to know if the math and overall procedure that I'm using is correct. $\endgroup$
    – jld
    Mar 23, 2015 at 13:27

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.