Implementing binary logistic regression from scratch Background knowledge:
To train a logistic regression model for a classification problem with two classes (called class $0$ and class $1$), we are given a training dataset consisting of feature vectors $x_1, x_2, \ldots, x_N \in \mathbb R^d$ and corresponding target values $y_1, y_2, \ldots, y_N \in \{0,1\}$. Our goal is to find numbers $\beta_0, \beta_1, \ldots, \beta_d \in \mathbb R$ such that
$$
\tag{1} y_i \approx \sigma(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_d x_{id}) \quad \text{for } i = 1, \ldots, N.
$$
Here $x_{i1}, \ldots, x_{id}$ are the components of the feature vector $x_i$, and $\sigma:\mathbb R \to \mathbb R$ is the sigmoid function (also called "logistic function") defined by
$$
\sigma(u) = \frac{e^u}{1 + e^u} \quad \text{for all } u \in \mathbb R.
$$
The sigmoid function is useful in machine learning because it converts a real number into a probability (that is, a number between $0$ and $1$).
Equation (1) can be expressed more concisely using vector notation: we hope that
$$
y_i \approx \sigma(\hat x_i^T \beta) \quad \text{for } i = 1, \ldots N
$$
where
$$
\hat x_i = \begin{bmatrix} 1 \\ x_{i1} \\ \vdots \\ x_{id} \end{bmatrix} \quad \text{and} \quad \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_d \end{bmatrix}.
$$
(The vector $\hat x_i$ is called an "augmented" feature vector. It's the vector you get by prepending a $1$ to $x_i$.)
I think of $\sigma(\hat x_i^T \beta)$ as being a "predicted probability" which tells you how strongly our model believes that example $i$ belongs to class $1$. And $y_i$ is
a "ground truth" probability which reflects certainty about whether or not example $i$ belongs to class $1$.
We will use the binary cross-entropy loss function
$$
\ell(p,q) = -p \log(q) - (1 - p) \log(1 - q)
$$
to measure how well a predicted probability $q \in (0,1)$ agrees with a ground truth probability $p \in [0,1]$. Here $\log$ denotes the natural logarithm. The vector $\beta$ will be chosen to minimize the average cross-entropy
$$
\tag{2} L(\beta) = \frac{1}{N} \sum_{i=1}^N \ell(y_i, \sigma(\hat x_i^T \beta)).
$$
We can minimize $L(\beta)$ using an optimization algorithm such as gradient descent, stochastic gradient descent, or Newton's method. In order to use such methods, we need to know how to compute the gradient of $L$. For Newton's method, we also need to know how to compute the Hessian of $L$.
Question:
How can we compute the gradient and the Hessian of the average cross-entropy function $L(\beta)$ defined in equation (2)?
The goal is to compute the gradient and the Hessian of $L$ with finesse, making good use of vector notation.
(I'll post an answer below.)
 A: Notice that
$$
L(\beta) = \frac{1}{N} \sum_{i=1}^N h_i(\hat x_i^T \beta)
$$
where $h_i: \mathbb R \to \mathbb R$ is the function defined by
$$
h_i(u) = \ell(y_i, \sigma(u)).
$$
The sigmoid function and the cross-entropy loss function are natural companions, and they are happy to be combined together into the function $h_i$. By the chain rule,
$$
L'(\beta) = \frac{1}{N} \sum_{i=1}^N h_i'(\hat x_i^T \beta) \hat x_i^T.
$$
So, to compute $L'(\beta)$, we only need to figure out how to evaluate $h_i'(u)$. But this must be easy because $h_i$ is a function of a single variable. Evaluating $h_i'(u)$ is a routine exercise in single-variable calculus.
Before we take the derivative of $h_i$, let's first simplify $h_i(u)$ as much as possible:
\begin{align}
h_i(u) &= -y_i \log(\sigma(u)) - (1 - y_i) \log(1 - \sigma(u)) \\
&= -y_i \log(\frac{e^u}{1+e^u}) - (1 - y_i) \log(\frac{1}{1 + e^u}) \\
&= -y_i (u - \log(1 + e^u)) + (1 - y_i) \log(1 + e^u) \\
&= - y_i u + \log(1 + e^u).
\end{align}
Look how much $h_i(u)$ simplified! The sigmoid function and the cross-entropy loss function are meant to be together. This is why PyTorch offers the "binary cross-entropy with logits" loss function.
Now we are ready to take the derivative:
$$
h_i'(u) = -y_i +\frac{e^u}{1 + e^u} = \sigma(u) - y_i.
$$
This is a beautiful and delightfully simple formula. Interpretation: If the predicted probability $\sigma(u)$ agrees perfectly with the ground truth probability $y_i$, then $h_i'(u)$ is zero, suggesting that no change to the value of $u$ is necessary.
Putting the above pieces together, we see that
$$
L'(\beta) = \frac{1}{N} \sum_{i=1}^N (\sigma(\hat x_i^T \beta) - y_i) \hat x_i^T.
$$
If we use the convention that the gradient is a column vector, then
$$
\nabla L(\beta) =  L'(\beta)^T = \frac{1}{N} \sum_{i=1}^N (\sigma(\hat x_i^T \beta) - y_i) \hat x_i.
$$

Next let's compute the Hessian of $L$. The Hessian matrix $H L(\beta)$ is the derivative of the function
$$g(\beta) = \nabla L(\beta)= \frac{1}{N} \sum_{i=1}^N  \hat x_i (\sigma(\hat x_i^T \beta) - y_i).
$$
To compute the derivative of $g$, it is helpful to know that
$$
\sigma'(u) = \sigma(u) - \sigma(u)^2,
$$
which is a beautiful result that you can easily check.
Again using the chain rule, we find that the derivative of $g$ is
\begin{align}
g'(\beta) &= \frac{1}{N} \sum_{i=1}^N \hat x_i \sigma'(\hat x_i^T \beta) \hat x_i^T \\
&= \frac{1}{N} \sum_{i=1}^N (\sigma(u_i) - \sigma^2(u_i)) \hat x_i \hat x_i^T
\end{align}
where $u_i = \hat x_i^T \beta$.
So we have found our Hessian matrix $HL(\beta) = g'(\beta)$. Having computed the gradient and Hessian of $L$, it is now straightforward to train a logistic regression model from scratch using gradient descent or Newton's method in a language such as Python.
