# How to correctly implement iteratively reweighted least squares algorithm for multiple logistic regression?

I'm confused about the iteratively reweighted least squares algorithm used to solve for logistic regression coefficients as described on page 121 of The Elements of Statistical Learning, 2nd Edition (Hastie, Tibshirani, Friedman 2009).

The final step of the process, after fitting a Taylor approximation to the log-likelihood of $N$ observations, is to solve the following weighted least squares problem:

$\beta^{new}\leftarrow argmin_{\beta}(\textbf{z}-\textbf{X}\beta)^T\textbf{W}(\textbf{z}-\textbf{X}\beta)$ $(1)$

by finding $\frac{\delta[(\textbf{z}-\textbf{X}\beta)^T\textbf{W}(\textbf{z}-\textbf{X}\beta)]}{\delta\beta_j}$, setting $\frac{\delta[(\textbf{z}-\textbf{X}\beta)^T\textbf{W}(\textbf{z}-\textbf{X}\beta)]}{\delta\beta_j}=0$, then solving for $\beta_j^{new}$,

where:

$\textbf{z}=\textbf{X}\beta^{old}+\textbf{W}^{-1}(\textbf{y}-\textbf{p})$,

$\textbf{W}=N\times{}N$ diagonal matrix of weights with $i$th diagonal element $p(x_i;\beta^{old})(1-p(x_i;\beta^{old}))$,

$\textbf{p}=$vector of fitted probabilities with $i$th element $p(x_i;\beta^{old})$,

$\textbf{y}=$vector of $y_i$ values,

$\textbf{X}=$matrix of $x_i$ values,

$\beta=$vector of coefficients $\beta_0,\beta_1,...,\beta_p$.

In the right hand part of expression (1), the $\beta$s are missing any superscript. Is $\beta$ presumed to be equal to $\beta^{old}$? That is, in order to solve for $\beta_j^{new}$ in (1) do we plug in the most current update of $\beta$ for all values of $\beta_{l\neq j}$ calculated in prior steps?

• Just a word of caution: While the algorithms outlined in ESL are generally ok from a descriptive standpoint, please don't rely on anything in there being numerically stable or robust. Mar 23, 2015 at 13:54

In an expression like

$$\beta^{new}\leftarrow \text{argmin}_{b}(\textbf{z}-\textbf{X}b)^T\textbf{W}(\textbf{z}-\textbf{X}b)$$

the point is that the output, $\beta^{new}$, is the result of considering all possible $b\in \mathbb{R}^p$ or whatever other space you are optimizing over. That's why there's no superscript: in the optimization problem $\beta$ is a dummy variable, just like with an integral (and I'm deliberately writing $b$ not $\beta$ to reflect $b$ being a dummy variable, not the target parameter).

The overall procedure involves getting a $\beta^{(t)}$, computing the "response" for the WLS, and then solving the WLS problem for $\beta^{(t+1)}$; as you know, we can use derivatives to get a nice closed-form solution for the optimal $\hat \beta$ for this problem. Thus $\beta^{old}$, which is fixed, appears in the vector $\textbf{z}$ in the WLS computation and then leads to $\beta^{new}$. That's the "iteration" part, that we use our current solution to create a new response vector; the WLS part then is solving for the new $\hat \beta$ vector. We keep doing this until there's no "significant" change.

Remember that the WLS procedure doesn't know that it is being used iteratively; as far as it is concerned, it is presented with an $X$, $y$, and $W$ and then outputs

$$\hat{\beta} = (X^T W X)^{-1} X^T W y$$ like it would do in any other instance. We are being clever with our choice of $y$ and $W$ and iterating.

Update: We can derive the solution to the WLS problem without using any component-wise derivatives. Note that if $Y \sim \mathcal N(X\beta, I)$ then $W^{1/2}Y \sim \mathcal N(W^{1/2}X\beta, W)$ from which we have that $$\frac{\text d}{\text d\beta}\|W^{1/2}Y - W^{1/2}X\beta\|^2 = -2X^TWY + 2X^TWX\beta.$$

Setting the derivative equal to 0 and solving we obtain

$$\hat{\beta} = (X^TWX)^{-1} X^TWY.$$

Thus for any inputs $W$, $X$, and $Y$ (provided W is positive definite and $X$ is full column rank) we get our optimal $\hat{\beta}$. It doesn't matter what these inputs are. So what we do is we use our $\beta^{old}$ to create our $Y$ vector and then we plug that in to this formula which outputs the optimal $\hat \beta$ for the given inputs. The whole point of the WLS procedure is to solve for $\hat \beta$. It in and of itself doesn't require plugging in a $\hat \beta$.

• -Thank you for the explanation. I understand that the $\beta$s are dummy variables in the sense that $\beta_{l\neq j}$s are treated as constants when finding the derivative $\frac{\delta[(\textbf{z}-\textbf{X}\beta)^T\textbf{W}(\textbf{z}-\textbf{X} \beta)]}{\delta\beta_j}$. However the solution to the derivative is going to be something like $2(\textbf{z}-\textbf{X}\beta)\textbf{W}\textbf{x}_j^T$. We're still left with a $\beta$ vector in the expression when solving for $\beta_j$ in the current iteration. Do we just plug in $\beta=\beta^{old}$? Mar 23, 2015 at 15:03
• Ok, I think I understand. When solving for the logistic regression coefficients, we substite $\textbf{z}$ for $\textbf{y}$ giving us: $\hat{\beta} = (X^TWX)^{-1} X^TW \textbf{z}$ where $\textbf{z}=\textbf{X}\beta^{old}+\textbf{W}^{-1}(\textbf{y}-\textbf{p})$. The $argmin$ expression is really more of a shorthand, is that correct? Mar 23, 2015 at 17:17
• Just to be clear, with what you're doing, $y$ is your actual response data where $z$ is the "response" as used in the WLS formulas (sometimes $z$ is called the working response) which is what the $Y$ in my answers refers to. The argmin comes from the fact that regression is a minimization problem: we find $\beta$ to minimize the residual sum of squares. That's why we can set derivatives equal to 0, because we're looking for a minimum. The argmin notation is a very convenient way to express this.
• That's correct, when applying IRLS to a logistic regression, I'm assuming $\textbf{y}$ is the response vector of 0s and 1s. Mar 23, 2015 at 17:58