I am trying to derive the gradient and hessian of logistic regression with ridge penalty. The log-likelihood should be (correct me if I am wrong):
$$\sum_{i=0}^n\Big(\log{(P_i^{y_i}(1-P_i)^{1-y_i}- \lambda\beta^T\beta)}\Big)$$ where: $$P_i = \frac{1}{1+e^{-x_i^T \beta}}$$
Is it possible to derive the gradient and hessian analytically to apply newton's method? I made it work using the autograd library from Python but I am curious about the actual analytical derivation.
We have $y_i \mid x_i \stackrel{\text{indep}}\sim \text{Bern}(\theta_i)$ and then we model $g(\theta) = X\beta$ where $g = \text{logit}$ is the link function and $X$ is our $n\times p$ full rank feature matrix.
I'm going to first work out the log likelihood in terms of $\theta$, and then I'll use the assumed relationship between $\theta$ and $\beta$ to get the actual likelihood of interest. We have $$ \ell(\theta\mid y, x) = \log \prod_{i=1}^n \theta_i^{y_i} (1-\theta_i)^{1-y_i} \\ = \sum_{i=1}^n \left(y_i \log\theta_i + (1-y_i)\log(1-\theta_i)\right) \\ = \sum_i \left(y_i \operatorname{logit}\theta_i + \log(1 - \theta_i)\right). $$ I'll let $h = g^{-1}$ be the inverse link function, so in this case $h(x) = \frac{1}{1+e^{-x}}$.
Now from $g(\theta_i) = x_i^T\beta$ we have $$ \ell(\beta\mid y,x) = \sum_i \left(y_i \cdot x_i^T\beta + \log\left(1 - \frac{1}{1+e^{-x_i^T\beta}}\right)\right) \\ = \sum_i \left(y_i \cdot x_i^T\beta - \log\left(1 + e^{x_i^T\beta}\right)\right) . $$ We can note that $\sum_i y_i x_i^T\beta = y^TX\beta$ and $\log(1 + e^x)$ is often called the softplus function so I'll denote it via $\newcommand{\sf}{\operatorname{sf}}$$\sf(x)$. I'll treat $\sf$ as being vectorized, so e.g. $\sf((1,2,3)) = (\sf(1), \sf(2), \sf(3))$. This means that I can write the log likelihood as $\newcommand{\one}{\mathbf 1}$ $$ \ell(\beta\mid y, x) = y^TX\beta - \one^T\sf(X\beta)). $$
Adding in the ridge penalty and negating the log likelihood to make this a minimization problem leads to a full objective function of $$ f(\beta\mid y, x) = -y^TX\beta + \one^T\sf(X\beta)) + \frac{\lambda}{2} \beta^T\beta. $$ I'm dividing $\lambda$ by $2$ just to catch the factor of $2$ that comes from taking derivatives.
Now we can do some calculus. We can work out $$ \frac{\partial}{\partial \beta}\one^T\sf(X\beta) = \sum_i \frac{\partial}{\partial \beta} \log(1 + e^{x_i^T\beta}) \\ = \sum_i \frac{e^{x_i^T\beta}}{1 + e^{x_i^T\beta}} \cdot x_i \\ = \sum_i h(x_i^T\beta) \cdot x_i \\ = X^Th(X\beta) $$ (I'm also vectorizing $h$ as needed). All together this means $$ \nabla f = -X^Ty + X^Th(X\beta) + \lambda\beta \\ = -X^T(y - \hat y) + \lambda \beta $$ where I'm writing $\hat y = h(X\beta)$ since these are our predicted probabilities for this choice of $\beta$. If we were just doing minimum negative log likelihood, i.e. $\lambda = 0$, then a stationary point will correspond to $X^T(y - \hat y) = \mathbf 0$. This is really natural since this is saying that we've got a good $\hat\beta$ when our residuals are orthogonal to our predictors.
For the Hessian we have $$ \nabla^2 f= \lambda I + \frac{\partial}{\partial \beta} X^Th(X\beta) \\ = \lambda I + \sum_i x_i \frac{\partial}{\partial \beta} h(x_i^T\beta) \\ = \lambda I + \sum_i x_i h'(x_i^T\beta)x_i^T. $$ Now we can use a key feature of the inverse logit, namely that $h' = h \cdot (1-h)$, to turn this into $$ H(\beta) = \lambda I + X^T\operatorname{diag}\left[\hat y \cdot (1 - \hat y)\right]X. $$
Note that $H$ doesn't have $y$ in it anymore. This is a feature of having used the canonical link, which is that the Hessian is equal to the expected value of the Hessian (w.r.t. the distribution of $y\mid X$) so Fisher scoring and Newton-Raphson are equivalent here.
This is also an interpretable form in that $X^T\operatorname{diag}\left[\hat y \cdot (1 - \hat y)\right]X$ is like the sample covariance matrix except each case is weighted by our confidence in it, since $\hat y_i (1 - \hat y_i)$ is the variance of a Bernoulli random variable. We then add the regularization of $\lambda I$ so we're just bounding the smallest eigenvalue away from zero and guaranteeing good conditioning.
Even without the ridge penalty we can see that $X^T\operatorname{diag}\left[\hat y \cdot (1 - \hat y)\right]X$ is positive definite, so this is a convex optimization, but if there is near-perfect separation then we can end up with $\hat y_i$ approaching near zero or one and this can make the Hessian close to singular which is an issue for inverting it in a Newton-Raphson step. The ridge penalty helps guarantee that we're never close to singular.
