Bayesian GLM#

Bayesian generalized linear model for non-Gaussian likelihoods (binary and count data). The posterior is not conjugate; it is approximated as a Gaussian at the MAP via Laplace approximation / IRLS.

See also: NormalRegressor (exact conjugate, Gaussian likelihood).

Symbols#

Symbol

Meaning

\(p\)

Number of features

\(n\)

Number of observations in a batch

\(\mathbf{w}\)

Weight vector (MAP estimate stored as coef_)

\(\alpha\)

Prior precision scalar

\(\boldsymbol{\Lambda}\)

Posterior precision matrix (stored as cov_inv_)

\(\boldsymbol{\eta}\)

Linear predictor \(\mathbf{X}\mathbf{w}\)

\(g^{-1}\)

Inverse link function (sigmoid or exp)

\(\boldsymbol{\mu}\)

Mean response \(g^{-1}(\boldsymbol{\eta})\)

\(\mathbf{W}\)

Diagonal IRLS weight matrix

\(\mathbf{z}\)

Working response (pseudo-targets)

\(\gamma\)

Decay factor, in \((0, 1]\)

Prior#

\[\mathbf{w} \sim \mathcal{N}(\mathbf{0},\; \alpha^{-1}\mathbf{I})\]

Likelihood#

Logit link (Bernoulli):

\[p(y_i = 1 \mid \mathbf{x}_i, \mathbf{w}) = \sigma(\mathbf{x}_i^\top \mathbf{w}) = \frac{1}{1 + \exp(-\mathbf{x}_i^\top \mathbf{w})}\]

Log link (Poisson):

\[\mathbb{E}[y_i \mid \mathbf{x}_i, \mathbf{w}] = \exp(\mathbf{x}_i^\top \mathbf{w})\]

Neither likelihood is conjugate to the Gaussian prior.

Reference: Murphy (2012) Chapter 8 [1].

Laplace approximation via IRLS#

The posterior is approximated as:

\[p(\mathbf{w} \mid \mathcal{D}) \approx \mathcal{N}(\hat{\mathbf{w}},\; \boldsymbol{\Lambda}^{-1})\]

where \(\hat{\mathbf{w}}\) is the MAP estimate and \(\boldsymbol{\Lambda}\) is the Hessian of the negative log-posterior evaluated at the MAP.

The MAP is found by iteratively reweighted least squares (IRLS). Each iteration:

  1. Compute the linear predictor and link derivatives:

    \[\begin{split}\boldsymbol{\eta} &= \mathbf{X}\mathbf{w} \\ \boldsymbol{\mu} &= g^{-1}(\boldsymbol{\eta}) \\ \frac{d\boldsymbol{\mu}}{d\boldsymbol{\eta}} &= \begin{cases} \boldsymbol{\mu}(1 - \boldsymbol{\mu}) & \text{logit} \\ \boldsymbol{\mu} & \text{log} \end{cases}\end{split}\]

  2. Form IRLS weights and working response:

    \[\begin{split}\mathbf{W} &= \mathrm{diag}\!\left( \frac{d\boldsymbol{\mu}}{d\boldsymbol{\eta}} \right) \\ \mathbf{z} &= \boldsymbol{\eta} + \frac{\mathbf{y} - \boldsymbol{\mu}} {d\boldsymbol{\mu}/d\boldsymbol{\eta}}\end{split}\]

  3. Solve the weighted least squares problem:

    \[\begin{split}\boldsymbol{\Lambda} &= \gamma^n\,\boldsymbol{\Lambda}_{\text{old}} + \mathbf{X}^\top \mathbf{W} \mathbf{X} \\ \mathbf{w} &\leftarrow \boldsymbol{\Lambda}^{-1}\!\left( \gamma^n\,\boldsymbol{\Lambda}_{\text{old}}\, \mathbf{w}_{\text{old}} + \mathbf{X}^\top (\mathbf{W} \odot \mathbf{z}) \right)\end{split}\]

  4. Check convergence: \(\|\mathbf{w}_{\text{new}} - \mathbf{w}_{\text{old}}\|_\infty < \texttt{tol}\).

The default runs 5 iterations per partial_fit call (LaplaceApproximator(n_iter=5)). For fast online updates where the previous posterior is a good initialization, n_iter=1 (a single Newton step) is often sufficient.

Sampling#

Weight vectors are drawn from the Gaussian approximation:

\[\mathbf{w}_s \sim \mathcal{N}(\hat{\mathbf{w}},\; \boldsymbol{\Lambda}^{-1})\]

Predictions are transformed through the inverse link:

\[\hat{\mathbf{y}} = g^{-1}(\mathbf{X}\,\mathbf{w}_s)\]

For the logit link, samples are probabilities in \((0, 1)\). For the log link, samples are positive rates.

Decay#

\[\boldsymbol{\Lambda} \leftarrow \gamma^n\,\boldsymbol{\Lambda}\]

Same as NormalRegressor: scales the precision, widens the posterior, mean unchanged. See Choosing and Tuning a Decay Rate.

Hyperparameter semantics#

Parameter

Controls

Practical guidance

alpha

Prior precision. Sets \(\boldsymbol{\Lambda}_0 = \alpha\mathbf{I}\).

1.0 (default). Higher values regularize more.

link

Likelihood family. 'logit' for binary outcomes, 'log' for counts.

Match to your outcome type.

approximator

Posterior approximation strategy. Default is LaplaceApproximator(n_iter=5, tol=1e-4).

See the IRLS section above for n_iter guidance.

learning_rate

Decay factor \(\gamma\). See Choosing and Tuning a Decay Rate.

1.0 (default) for stationary environments.

Trade-offs vs. conjugate models#

The Laplace approximation buys flexible likelihoods (logistic, Poisson) at a cost:

  • The posterior is Gaussian at the MAP, not exact. Tail behavior and multimodality are lost.

  • IRLS iterations are more expensive than a single conjugate update (\(O(n\_iter \cdot p^2 n)\) vs. \(O(p^2 n)\)).

  • With n_iter=1, the cost matches a conjugate update, but the approximation is only good when the previous posterior is close to the new MAP.

For Gaussian outcomes, NormalRegressor is exact and cheaper. For binary or count outcomes, NormalRegressor is also viable: by Bernstein-von Mises, the linear model’s posterior concentrates correctly even under likelihood misspecification, and in practice it performs well (see the Comparing Bayesian Approaches to Disjoint Linear Bandits example). The GLM buys you a proper likelihood at the cost of IRLS iterations.

References#