Bernoulli Setting

Assume $Y_i$ follows Bernoulli distribution given the $i^{th}$ observation ${\mathbf{x}}_i$ and the parameters $\mathit{\beta }$ $$Y_i\mid {\mathbf{x}}_i\sim \text{Bernoulli}(p_i)$$ We can write the probability mass function for $y_i$ as follows $$\mathbb{P}(Y_i=y_i\mid {\mathbf{x}}_i,p_i)=p_i^{y_i}(1-p_i{)}^{1-y_i}$$ We assume that the log-odds is a linear combination of the input $$\ln \left(\frac{p_i}{1-p_i}\right)={\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }\text{i.e.}{p}_i=\frac{1}{1+\exp (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}=\frac{\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}=\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })$$

Joint Log-Likelihood

The joint likelihood is then found assuming IID-ness $$\begin{array}{rl}p(\mathbf{y}\mid \mathit{\beta })& =\prod _{i=1}^n\mathbb{P}(Y_i=y_i\mid {\mathbf{x}}_i,p_i)\\ & =\prod _{i=1}^n{p}_i^{y_i}(1-p_i{)}^{1-y_i}\end{array}$$ The log-likelihood is then $$\begin{array}{rl}\ln p(\mathbf{y}\mid \mathit{\beta })& =\ln (\prod _{i=1}^n{p}_i^{y_i}(1-p_i{)}^{1-y_i})\\ & =\sum _{i=1}^n{y}_i\ln (p_i)+(1-y_i)\ln (1-p_i)\end{array}$$ Alternatively, the likelihood can be written as follows $$\begin{array}{rl}p(\mathbf{y}\mid \mathit{\beta })& =\prod _{i=1}^n{p}_i^{y_i}(1-p_i{)}^{1-y_i}\\ & =\prod _{i=1}^n{\left(\frac{\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right)}^{y_i}{(1-\frac{\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })})}^{1-y_i}\\ & =\prod _{i=1}^n{\left(\frac{\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right)}^{y_i}{\left(\frac{1}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right)}^{1-y_i}\\ & =\prod _{i=1}^n{\left(\frac{\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right)}^{y_i}\left(\frac{1}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right){(1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))}^{y_i}\\ & =\prod _{i=1}^n\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{y}_i)\frac{1}{{(1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))}^{y_i}}\left(\frac{1}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right){(1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))}^{y_i}\\ & =\prod _{i=1}^n\frac{\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{y}_i)}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\end{array}$$

Taking the logarithm of this expression gives $$\begin{array}{rl}\ln (p(\mathbf{y}\mid \mathit{\beta }))& =\sum _{i=1}^n\ln \left(\frac{\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{y}_i)}{1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right)\\ & =\sum _{i=1}^n\ln (\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{y}_i))-\ln (1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))\\ & =\sum _{i=1}^n{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{y}_i-\ln (1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))\end{array}$$

Maximum Likelihood $\equiv$ Minimize Loss Function

Our aim is to maximize the likelihood. This is equivalent to maximizing the log likelihood, which is equivalent to minimizing the negative log likelihood. $$\underset{\mathit{\beta }}{min}-\sum _{i=1}^n{y}_i\ln (p_i)+(1-y_i)\ln (1-p_i)$$ Let’s consider the expression inside the summation $$y_i\ln (p_i)+(1-y_i)\ln (1-p_i)$$ We can notice that $$y_i\ln (p_i)+(1-y_i)\ln (1-p_i)=\{\begin{array}{lll}\ln (1-p_i)& & \text{if }{y}_i=0\\ \ln (p_i)& & \text{if }{y}_i=1\end{array}$$ Remember that $p_i=\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })$ and that $1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })=\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })$ $$y_i\ln (\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+(1-y_i)\ln (\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))=\{\begin{array}{lll}\ln (\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))& & \text{if }{y}_i=0\\ \ln (\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))& & \text{if }{y}_i=1\end{array}$$ We are thus looking for something that is $-1$ when $y_i=0$ and that it is $1$ when $y_i=1$. In particular, notice that $2y_i-1$ does the job. Thus we can write our problem as $$\begin{array}{rl}\underset{\mathit{\beta }}{min}-\sum _{i=1}^n\ln (\sigma ((2y_i-1){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))& =\underset{\mathit{\beta }}{min}\sum _{i=1}^n-\ln \left(\frac{1}{1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}\right)\\ & =\underset{\mathit{\beta }}{min}\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))\\ & =\underset{\mathit{\beta }}{min}\sum _{i=1}^{n}L({\mathbf{x}}_i^{\mathrm{\top }},y_i;\mathit{\beta })\end{array}$$ where the loss incurred using parameter $\mathit{\beta }$ when predicting the label for ${\mathbf{x}}_i^{\mathrm{\top }}$ whose true label is $y_i$ is $$L({\mathbf{x}}_i^{\mathrm{\top }},y_i;\mathit{\beta })=\{\begin{array}{lll}\ln (1+\exp ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))& & \text{if }{y}_i=0\\ \ln (1+\exp (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })& & \text{if }{y}_i=1\end{array}$$

Maximum-A-Posteriori (MAP)

Recall from Bayes Rule $$p(\mathit{\beta }\mid \mathbf{y})=\frac{p(\mathbf{y}\mid \mathit{\beta })p(\mathit{\beta })}{p(\mathbf{y})}\propto p(\mathbf{y}\mid \mathit{\beta })p(\mathit{\beta })$$ Taking the logarithm both sides and multiplying by $-1$ we obtain $$-\ln (p(\mathit{\beta }\mid \mathbf{y}))\propto -\ln (p(\mathbf{y}\mid \mathit{\beta }))-\ln (p(\mathit{\beta }))$$ We can choose to use a Gaussian Prior on the parameters $p(\mathit{\beta })=N({\mathit{\mu }}_0,{\mathbf{\Sigma }}_0)$. If $\mathit{\beta }\in {\mathbb{R}}^{p\times 1}$ then $$\ln (p(\mathit{\beta }))=-\frac{p}{2}\ln (2\pi )-\frac{1}{2}\ln (|{\mathrm{\Sigma }}_0|)-\frac{1}{2}(\mathit{\beta }-{\mathit{\mu }}_0{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_0^{-1}(\mathit{\beta }-{\mathit{\mu }}_0)$$ Plugging this in, we obtain the following. Notice how we neglect terms that do not depend on $\mathit{\beta }$ because they will not matter when we minimize this with respect to $\mathit{\beta }$. $$\begin{array}{rl}\underset{\mathit{\beta }}{min}-\ln (p(\mathit{\beta }\mid \mathbf{y}))& =\underset{\mathit{\beta }}{min}-\ln (p(\mathbf{y}\mid \mathit{\beta }))-\ln (p(\mathit{\beta }))\\ & =\underset{\mathit{\beta }}{min}-\ln (p(\mathbf{y}\mid \mathit{\beta }))+\frac{p}{2}\ln (2\pi )+\frac{1}{2}\ln (|{\mathrm{\Sigma }}_0|)+\frac{1}{2}(\mathit{\beta }-{\mathit{\mu }}_0{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_0^{-1}(\mathit{\beta }-{\mathit{\mu }}_0)\\ & =\underset{\mathit{\beta }}{min}-\ln (p(\mathbf{y}\mid \mathit{\beta }))+\frac{1}{2}(\mathit{\beta }-{\mathit{\mu }}_0{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_0^{-1}(\mathit{\beta }-{\mathit{\mu }}_0)\\ & =\underset{\mathit{\beta }}{min}\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+\frac{1}{2}(\mathit{\beta }-{\mathit{\mu }}_0{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_0^{-1}(\mathit{\beta }-{\mathit{\mu }}_0)\end{array}$$

Ridge Regularization $\equiv$ Isotripic Gaussian Prior

Now if we set ${\mathit{\mu }}_0={\mathbf{0}}_p$ and ${\mathbf{\Sigma }}_0={\sigma }_{\mathit{\beta }}^2{I}_p$ (this is equivalent to setting a univariate normal prior on each coefficient, with $p({\beta }_j)=N(0,{\sigma }_{\mathit{\beta }}^2)$) we have $$\underset{\mathit{\beta }}{min}\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+\frac{1}{2}{\mathit{\beta }}^{\mathrm{\top }}({\sigma }_{\mathit{\beta }}^2{I}_p{)}^{-1}\mathit{\beta }$$ using the fact that for an invertible matrix $A$ and a non-zero constant $c\in \mathbb{R}\mathrm{\setminus }\{0\}$ we have $(cA{)}^{-1}=\frac{1}{c}{A}^{-1}$ we obtain $$\underset{\mathit{\beta }}{min}\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+\frac{1}{2{\sigma }_{\mathit{\beta }}^2}{\mathit{\beta }}^{\mathrm{\top }}\mathit{\beta }$$ Setting $\lambda :=\frac{1}{{\sigma }_{\mathit{\beta }}^2}$ we have regularized logistic regression $$\underset{\mathit{\beta }}{min}\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+\frac{\lambda }{2}{\mathit{\beta }}^{\mathrm{\top }}\mathit{\beta }$$ It is more stable to multiply out by ${\sigma }_{\mathit{\beta }}^2$ so we get $$\underset{\mathit{\beta }}{min}{\sigma }_{\mathit{\beta }}^2\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+\frac{1}{2}{\mathit{\beta }}^{\mathrm{\top }}\mathit{\beta }$$

Ridge Regularization except on Intercept

Where $$\mathit{\beta }=\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_{p-1}\end{array}\right):=\left(\begin{array}{c}{\beta }_0\\ {\mathit{\beta }}_{1:p-1}\end{array}\right)$$ We’ve defined ${\mathit{\beta }}_{1:p-1}:=({\beta }_1,\dots ,{\beta }_p{)}^{\mathrm{\top }}$ because often we don’t really regularize the intercept. This means that we place a Multivariate Gaussian prior on ${\mathit{\beta }}_{1:p-1}$ as follows $p({\mathit{\beta }}_{1:p-1})=N({\mathbf{0}}_{p-1},{\sigma }_{{\mathit{\beta }}_{1:p-1}}^2{I}_{p-1})$ (again, this is equivalent to putting a univariate normal prior on each of ${\beta }_1,\dots ,{\beta }_{p-1}$ with $p({\beta }_j)=N(0,{\sigma }_{{\mathit{\beta }}_{1:p-1}}^2)$). Instead, on ${\beta }_0$ we could place a uniform prior, which means it’s value would not depend on ${\beta }_0$ and so could be dropped from the expression. $$\underset{\mathit{\beta }}{min}\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+\frac{1}{2{\sigma }_{{\mathit{\beta }}_{1:p-1}}^2}{\mathit{\beta }}_{1:p-1}^{\mathrm{\top }}{\mathit{\beta }}_{1:p-1}$$ It is more stable to multiply out by ${\sigma }_{{\mathit{\beta }}_{1:p-1}}^2$ therefore $$\underset{\mathit{\beta }}{min}{\sigma }_{{\mathit{\beta }}_{1:p-1}}^2\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+\frac{1}{2}{\mathit{\beta }}_{1:p-1}^{\mathrm{\top }}{\mathit{\beta }}_{1:p-1}$$ Notice that this is consistent with the implementation used by scikit-learn provided here.

Full-Bayesian (Laplace Approximation)

Full-Bayesian in intractable. Laplace Approximation approximates $p(\mathit{\beta }\mid \mathbf{y})$ with a Gaussian distribution $q(\mathit{\beta })$. To find such a distribution, we use the multivariate version of Taylor’s expansion to expand the log posterior around its mode ${\mathit{\beta }}_{\text{MAP}}$. We take a second order approximation $$\ln p(\mathit{\beta }\mid \mathbf{y})\simeq \ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})+\mathrm{\nabla }\ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}})+\frac{1}{2}(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}}{)}^{\mathrm{\top }}{\mathrm{\nabla }}^2\ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}})$$ Since ${\mathit{\beta }}_{\text{MAP}}$ is a stationary point, the gradient at this point will be zero so we have $$\ln p(\mathit{\beta }\mid \mathbf{y})\simeq \ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})+\frac{1}{2}(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}}{)}^{\mathrm{\top }}{\mathrm{\nabla }}^2\ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}})$$ We take the exponential both sides to obtain $$p(\mathit{\beta }\mid \mathbf{y})\simeq p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})\exp (\frac{1}{2}(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}}{)}^{\mathrm{\top }}{\mathrm{\nabla }}^2\ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}}))$$ We regnize this has the shape of a Multivariate Normal distribution. We therefore define our Laplace approximation to be $$q(\mathit{\beta })=(2\pi {)}^{-\frac{p}{2}}\text{det}(-{\mathrm{\nabla }}^2\ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y}))\exp (\frac{1}{2}(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}}{)}^{\mathrm{\top }}{\mathrm{\nabla }}^2\ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})(\mathit{\beta }-{\mathit{\beta }}_{\text{MAP}}))$$ That is $$q(\mathit{\beta })=N({\mathit{\beta }}_{\text{MAP}},{[-{\mathrm{\nabla }}^2\ln p({\mathit{\beta }}_{\text{MAP}}\mid \mathbf{y})]}^{-1})$$ To find ${\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathit{\beta }\mid \mathbf{y}){|}_{\mathit{\beta }={\mathit{\beta }}_{\text{MAP}}}$ we write $${\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathit{\beta }\mid \mathbf{y})={\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathbf{y}\mid \mathit{\beta })+{\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathit{\beta })$$ and we find each of the expressions on the right-hand side separately. We start with ${\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathit{\beta })$. Recall that if we have a quadratic form ${\mathbf{x}}^{\mathrm{\top }}A\mathbf{x}$ the derivative of this quadratic form with respect to $\mathbf{x}$ is given by ${\mathbf{x}}^{\mathrm{\top }}(A+A^{\mathrm{\top }})$. Applying this to our case, and using the fact that ${\mathbf{\Sigma }}_0^{-1}$ is symmetric we have $${\mathrm{\nabla }}_{\mathit{\beta }}\ln p(\mathit{\beta })=-\frac{1}{2}(\mathit{\beta }-{\mathit{\mu }}_0{)}^{\mathrm{\top }}2{\mathbf{\Sigma }}_0^{-1}=-(\mathit{\beta }-{\mathit{\mu }}_0{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_0^{-1}=-{\mathit{\beta }}^{\mathrm{\top }}{\mathbf{\Sigma }}_0^{-1}+{\mathit{\mu }}_0^{\mathrm{\top }}{\mathbf{\Sigma }}_0^{-1}$$ Taking the derivative with respect to $\mathit{\beta }$ again we get $$-{\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathit{\beta })={\mathbf{\Sigma }}_0^{-1}$$ To find ${\mathrm{\nabla }}_{\mathit{\beta }}\ln p(\mathbf{y}\mid \mathit{\beta })$ we take the derivative componentwise $$\begin{array}{rl}\frac{\partial }{\partial {\beta }_j}\ln p(\mathbf{y}\mid \mathit{\beta })& =\sum _{i=1}^n{y}_i{\partial }_{{\beta }_j}\ln \sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })+(1-y_i){\partial }_{{\beta }_j}\ln \sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })\\ & =\sum _{i=1}^n{y}_i\frac{\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))}{\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}{x}_i^{(j)}+(1-y_i)\frac{\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))}{\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })}(-x_i^{(j)})\\ & =\sum _{i=1}^n{y}_i{x}_i^{(j)}\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })+(y_i{x}_i^{(j)}-x_i^{(j)})\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })\end{array}$$ Now take the derivative with respect to ${\beta }_k$ $$\begin{array}{rl}{\partial }_{{\beta }_k}{\partial }_{{\beta }_j}\ln p(\mathbf{y}\mid \mathit{\beta })& =\sum _{i=1}^n{y}_i{x}_i^{(j)}{\partial }_{{\beta }_k}\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })+(y_i{x}_i^{(j)}-x_i^{(j)}){\partial }_{{\beta }_k}\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })\\ & =\sum _{i=1}^n{y}_i{x}_i^{(j)}\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))(-x_i^{(k)})+(y_i{x}_i^{(j)}-x_i^{(j)})\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))x_i^{(k)}\\ & =\sum _{i=1}^n-y_i{x}_i^{(j)}{x}_i^{(k)}\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+(y_i{x}_i^{(j)}{x}_i^{(k)}-x_i^{(j)}{x}_i^{(k)})\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))\\ & =-\sum _{i=1}^n{x}_i^{(j)}{x}_i^{(k)}\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))\end{array}$$ This tells us that $$[{\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathbf{y}\mid \mathit{\beta }){]}_{kj}=-\sum _{i=1}^n{x}_i^{(j)}{x}_i^{(k)}\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))$$ Note that for a vector ${\mathbf{x}}_i:=(1,x_i^{(1)},\dots ,x_i^{(p-1)}{)}^{\mathrm{\top }}$ the outer product gives the following symmetric matrix $${\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}=\left(\begin{array}{c}1\\ x_i^{(1)}\\ ⋮\\ x_i^{(p-1)}\end{array}\right)\left(\begin{array}{cccc}1& x_i^{(1)}& \cdots & x_i^{(p-1)}\end{array}\right)=\left(\begin{array}{cccc}1& x_i^{(1)}& \cdots & x_i^{(p-1)}\\ x_i^{(1)}& (x_i^{(1)}{)}^2& \cdots & x_i^{(1)}{x}_i^{(p-1)}\\ ⋮& ⋮& \ddots & ⋮\\ x_i^{(p-1)}& x_i^{(p-1)}{x}_i^{(1)}& \cdots & (x_i^{(p-1)}{)}^2\end{array}\right)$$ In particular $[{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}{]}_{kj}=x_i^{(j)}{x}_i^{(k)}$ so that $$-{\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathbf{y}\mid \mathit{\beta })=\sum _{i=1}^n\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })){\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}$$ Putting everything together we get $$-{\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathit{\beta }\mid \mathbf{y})={\mathbf{\Sigma }}_0^{-1}+\sum _{i=1}^n\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })){\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}$$

Gradient Ascent Optimization (MLE, No Regularization)

The simplest way to find $\mathit{\beta }$ that maximizes the likelihood is to do gradient ascent. Remember that when working on the Laplace approximation we found the derivative of the log-likelihood with respect to the $j^{th}$ component of $\mathit{\beta }$. We can rearrange such an expression to get a nicer form. $$\begin{array}{rl}\frac{\partial }{\partial {\beta }_j}\ln p(\mathbf{y}\mid \mathit{\beta })& =\sum _{i=1}^n{y}_i{x}_i^{(j)}\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })+(y_i{x}_i^{(j)}-x_i^{(j)})\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })\\ & =\sum _{i=1}^n{y}_i{x}_i^{(j)}(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))+(y_i{x}_i^{(j)}-x_i^{(j)})\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })\\ & =\sum _{i=1}^n{y}_i{x}_i^{(j)}-y_i{x}_i^{(j)}\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })+y_i{x}_i^{(j)}\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })-x_i^{(j)}\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })\\ & =\sum _{i=1}^n[y_i-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })]x_i^{(j)}\end{array}$$ Therefore the full gradient is given by $$\begin{array}{rl}{\mathrm{\nabla }}_{\mathit{\beta }}\ln p(\mathbf{y}\mid \mathit{\beta })& ={(\frac{\partial }{\partial {\beta }_0}\ln p(\mathbf{y}\mid \mathit{\beta }),\dots ,\frac{\partial }{\partial {\beta }_{p-1}}\ln p(\mathbf{y}\mid \mathit{\beta }))}^{\mathrm{\top }}\\ & ={(\sum _{i=1}^n(y_i-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))x_i^{(0)},\dots ,\sum _{i=1}^n(y_i-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))x_i^{(p-1)})}^{\mathrm{\top }}\\ & =\sum _{i=1}^n(y_i-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))(x_i^{(0)},\dots ,x_i^{(p-1)})\\ & =\sum _{i=1}^n(y_i-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })){\mathbf{x}}_i\end{array}$$ Gradient ascent then becomes $${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+\gamma {\mathrm{\nabla }}_{\mathit{\beta }}\ln p(\mathbf{y}\mid \mathit{\beta })={\mathit{\beta }}_k+\gamma \sum _{i=1}^n(y_i-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })){\mathbf{x}}_i$$ This can be vectorized when programming as follows $${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+\gamma X^{\mathrm{\top }}(\mathbf{y}-\sigma (X\mathit{\beta }))$$ where $X\in {\mathbb{R}}^{n\times p}$ is the design matrix. $$X=\left(\begin{array}{c}{\mathbf{x}}_1^{\mathrm{\top }}\\ ⋮\\ {\mathbf{x}}_n^{\mathrm{\top }}\end{array}\right)$$

One can change the step size at every iteration. One possible choice for ${\gamma }_k$ is as follows $${\gamma }_k=\frac{|({\mathit{\beta }}_k-{\mathit{\beta }}_{k-1}{)}^{\mathrm{\top }}[\mathrm{\nabla }\ln p(\mathbf{y}\mid {\mathit{\beta }}_k)-\mathrm{\nabla }\ln p(\mathbf{y}\mid {\mathit{\beta }}_{k-1})]|}{\parallel \mathrm{\nabla }\ln p(\mathbf{y}\mid {\mathit{\beta }}_k)-\mathrm{\nabla }\ln p(\mathbf{y}\mid {\mathit{\beta }}_{k-1}{\parallel }^2}$$

Newton’s Method (MLE, No Regularization)

Again, during the Laplace approximation section we found that the Hessian is given by $${\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathbf{y}\mid \mathit{\beta })=-\sum _{i=1}^n\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })){\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}$$ This expression can be vectorized as follows $${\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathbf{y}\mid \mathit{\beta })=-X^{\mathrm{\top }}DX$$ where $$D=\text{diag}(\sigma (X\mathit{\beta })(1-\sigma (X\mathit{\beta })))=\left(\begin{array}{cccc}\sigma ({\mathbf{x}}_1^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_1^{\mathrm{\top }}\mathit{\beta }))& 0& \dots & 0\\ 0& \sigma ({\mathbf{x}}_2^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_2^{\mathrm{\top }}\mathit{\beta }))& \dots & 0\\ ⋮& \dots & \ddots & ⋮\\ 0& 0& \dots & \sigma ({\mathbf{x}}_n^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_n^{\mathrm{\top }}\mathit{\beta }))\end{array}\right)$$ Newton’s method then updates the weights as follows (where $\alpha$ is a learning rate to control convergence) $${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+\alpha (X^{\mathrm{\top }}DX{)}^{-1}{X}^{\mathrm{\top }}(\mathbf{y}-\sigma (X{\mathit{\beta }}_k))$$

Of course in practice we never invert the matrix but rather compute the direction $\mathbf{d}$ by solving the linear system $$X^{\mathrm{\top }}DX{\mathbf{d}}_k=\alpha X^{\mathrm{\top }}(\mathbf{y}-\sigma (X{\mathit{\beta }}_k))$$ and then we find the next iterate as $${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+{\mathbf{d}}_k$$

Gradient Ascent (MAP, Ridge Regularization)

We want to maximize $$\ln p(\mathit{\beta }\mid \mathbf{y})=\ln p(\mathbf{y}\mid \mathit{\beta })+\ln p(\mathit{\beta })=-\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))-\frac{1}{2{\sigma }_{\mathit{\beta }}^2}{\mathit{\beta }}^{\mathrm{\top }}\mathit{\beta }$$ The gradient of the log posterior is given by $${\mathrm{\nabla }}_{\mathit{\beta }}\ln p(\mathit{\beta }\mid \mathbf{y})=X^{\mathrm{\top }}(\mathbf{y}-\sigma (X\mathit{\beta }))-\frac{1}{{\sigma }_{\mathit{\beta }}^2}\mathit{\beta }$$ Thus gradient descent with regularization to do MAP becomes $${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+\gamma ({\sigma }_{\mathit{\beta }}^2{X}^{\mathrm{\top }}(\mathbf{y}-\sigma (X{\mathit{\beta }}_k))-{\mathit{\beta }}_k)$$

Newton’s Method (MAP, Ridge Regularization)

Similarly, we want to maximize $\ln p(\mathit{\beta }\mid \mathbf{y})$. The Hessian is given by $${\mathrm{\nabla }}_{\mathit{\beta }}^2\ln p(\mathit{\beta }\mid \mathbf{y})=-X^{\mathrm{\top }}DX-\frac{1}{{\sigma }_{\mathit{\beta }}^2}I$$ therefore the Newton’s Method update formula becomes $${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+\alpha {[{\sigma }_{\mathit{\beta }}^2{X}^{\mathrm{\top }}DX+I]}^{-1}({\sigma }_{\mathit{\beta }}^2{X}^{\mathrm{\top }}(\mathbf{y}-\sigma (X{\mathit{\beta }}_k))-{\mathit{\beta }}_k)$$

Iteratively Reweighted Least-Squares

We can manipulate the expression in Newton’s method by defining a new variable $${\mathbf{z}}_k=X{\mathit{\beta }}_k+D^{-1}(\mathbf{y}-\sigma (X{\mathit{\beta }}_k))$$ Then the updates take the form $${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+(X^{\mathrm{\top }}{D}_{k}X{)}^{-1}{X}^{\mathrm{\top }}{D}_k{\mathbf{z}}_k$$ In practice we would follow these steps

• Evaluate ${\mathit{\eta }}_k=X{\mathit{\beta }}_k$ and $D_k$.
• Solve the system $D_k{\mathbf{r}}_k=\mathbf{y}-\sigma ({\mathit{\eta }}_k)$ for ${\mathbf{r}}_k$.
• Compute ${\mathbf{z}}_k={\mathit{\eta }}_k+{\mathbf{r}}_k$.
• Solve the system $(X^{\mathrm{\top }}{D}_{k}X){\mathbf{d}}_k=X^{\mathrm{\top }}{D}_k{\mathbf{z}}_k$ for ${\mathbf{d}}_k$.
• Compute ${\mathit{\beta }}_{k+1}\leftarrow {\mathit{\beta }}_k+{\mathbf{d}}_k$

Alternatively, noticing that $D_{ii}=\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })(1-\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))>0$ one can take the square root of its elements and rewriting the problem as $$(D_k^{1/2}X{)}^{\mathrm{\top }}(D_k^{1/2}X){\mathbf{d}}_k=(D_k^{1/2}X{)}^{\mathrm{\top }}(D_k^{1/2}{\mathbf{z}}_k)$$ which is a simple least squares problem on the new variables ${\tilde{X}}_k=D_k^{1/2}X$ and ${\tilde{\mathbf{z}}}_k=D_k^{1/2}{\mathbf{z}}_k$, and can be solve by using the QR decomposition of $\tilde{X}=QR$ by solving the following system for ${\mathbf{d}}_k$ $$R{\mathbf{d}}_k=Q^{\mathrm{\top }}{\tilde{\mathbf{z}}}_k$$

Logistic Regression $\{-1,1\}$

Notice that in the previous chapter we used $y_i\in \{0,1\}$ with $$\mathbb{P}(Y_i=y_i\mid {\mathbf{x}}_i)=\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{)}^{y_i}\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{)}^{1-y_i}=\{\begin{array}{lll}\sigma (-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })& & \text{if }{y}_i=0\\ \sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })& & \text{if }{y}_i=1\end{array}$$ In particular $$p(y_i=1)=\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })$$ This gave us the loss function $$\sum _{i=1}^n\ln (1+\exp ((1-2y_i){\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }))$$ Now the key point is to notice that $$1-2y_i=\{\begin{array}{lll}1& & \text{if }{y}_i=0\\ -1& & \text{if }{y}_i=1\end{array}$$ So that actually the mapping that makes $\{0,1\}$-Logistic Regression and $\{-1,1\}$-Logistic Regression equivalent is $0\mapsto 1$ and $1\mapsto -1$.

Now instead we have $$p(z_i=-1)=\sigma ({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })$$