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.
Best Answer
The relationship is as follows: $l(\beta) = \sum_i L(z_i)$.
Define a logistic function as $f(z) = \frac{e^{z}}{1 + e^{z}} = \frac{1}{1+e^{-z}}$. They possess the property that $f(-z) = 1-f(z)$. Or in other words:
$$ \frac{1}{1+e^{z}} = \frac{e^{-z}}{1+e^{-z}}. $$
If you take the reciprocal of both sides, then take the log you get:
$$ \ln(1+e^{z}) = \ln(1+e^{-z}) + z. $$
Subtract $z$ from both sides and you should see this:
$$ -y_i\beta^Tx_i+ln(1+e^{y_i\beta^Tx_i}) = L(z_i). $$
Edit:
At the moment I am re-reading this answer and am confused about how I got $-y_i\beta^Tx_i+ln(1+e^{\beta^Tx_i})$ to be equal to $-y_i\beta^Tx_i+ln(1+e^{y_i\beta^Tx_i})$. Perhaps there's a typo in the original question.
Edit 2:
In the case that there wasn't a typo in the original question, @ManelMorales appears to be correct to draw attention to the fact that, when $y \in \{-1,1\}$, the probability mass function can be written as $P(Y_i=y_i) = f(y_i\beta^Tx_i)$, due to the property that $f(-z) = 1 - f(z)$. I am re-writing it differently here, because he introduces a new equivocation on the notation $z_i$. The rest follows by taking the negative log-likelihood for each $y$ coding. See his answer below for more details.