# Implementing kernel logistic regression using IRWLS

I am referring to [1] for implementing Kernel Logistic Regression using IRWLS. In logistic regression, the form of the regularized negative log-likelihood we aim to minimize is the following:

$L(w) = - \sum_{i = 1}^{l} t_{i} \log \mu_{i} + (1 - t_{i}) \log(1 - \mu_{i}) + \lambda ||w||^{2},$

where $\mu_{i} = P(y_{i} = 1 \mid x_{i}) = {\exp(w^{T} x_{i})}/[1 + \exp(w^{T} x_{i})]$, and $t_{i} \in \{0, 1\}$ is the actual "label" associated to the sample (everything as in usual logistic regression; we can see the bias term as an additional 1 in $x$).

When using IRWLS (Newton-Raphson), learning $w$ should boil down to the following iterative procedure:

$$w^\text{new} = w^\text{old} - (L''(w))^{-1} L'(w) = (X^{T} W X + \lambda I)^{-1} X^{T} W z,$$

with $z = X w^\text{old} + W^{-1}(t - \mu)$, and $W$ is a diagonal matrix where $W_{ii} = \mu_{i} (1 - \mu_{i})$. Everything to here seems OK for me.

In the kernelized version, $w^{T} x$ becomes $\sum_{i = 1}^{l} \alpha_{i} K(x_{i}, x) + b$, and [1] provides the following procedure for iteratively finding $\alpha$:

$$\alpha^{(k)} = (K^{T} W K + \lambda K)^{-1} K^{T} W z$$

where $z = (K \alpha^{(k - 1)} + W^{-1} (t - \mu))$.

My questions are the following, in decreasing order of importance (to me):

1. How can I find the bias term $b$ ?
2. Which are good starting values for $W$ and $\alpha$ ?
3. Can $(K^{T} W K + \lambda K)$ turn out to be non-invertible?

[1] Zhu, J. et al. - Kernel Logistic Regression and the Import Vector Machine - NIPS'01

• If $K$ is invertible, then $K^T W K + \lambda K$ is invertible: $W$ is positive-definite, so $K^T W K$ should be also. And we assumed $K$ is positive definite, so for $\lambda > 0$ we have that $\lambda K$ is positive definite. – Dougal Jul 15 '13 at 1:17
• Thanks Dougal: I was getting some NANs from the inverse, and I was afraid there was something breaking the PSD-ness (may it depend on machine precision?). However, the gradient of the NLL in KLR is almost the same as in traditional LR: the only thing to be aware of is that the regularization term in KLR's likelihood appears as $\frac{\lambda}{2} \alpha^{T} K \alpha$ (the norm in the RKHS). – Morra Cinese Jul 15 '13 at 2:46
• Actually, I just realized that if any of the $w^T x_i$ are 0, then $W$ will be singular. Also, in practice it might become arbitrarily close to singular even if it's not actually singular. You are using a solve call rather than an explicit inverse, right? – Dougal Jul 15 '13 at 13:09
• You were right, that can be seen as a sysetem with linear equations; I was actually inverting, and using an SVD-based pseudoinverse (inspired to MATLAB's pinv.m) for the cases with troubles. For deriving the bias term $b$ in the calculations, I'm now using $[K 1]$ instead of $K$, and $\lambda [ K 0 ; 0 0 ]$ instead of $\lambda K$. Thanks for your hint :) I'm switching to solve() right now! – Morra Cinese Jul 15 '13 at 18:03
• Further realization: if $K$ is positive definite and $W$ is singular, for $\lambda > 0$ we still have $K^T W K + \lambda K$ being positive definite. But yeah, in general you very rarely want to actually compute the full inverse; solve is more stable and faster. – Dougal Jul 15 '13 at 20:53

I think this might help, although I'm still unsure myself on the relationship between kernel logistic regression and good old generalised additive models with locally weighted regression smooths. The NIPS paper suggests they're at least strongly related, if not the same.

You can think of kernel logistic regression as fitting a weighted logistic regression for each data point $x_i$, based on its neighbours $x_j$. The weights are given by $K(||x_j - x_i||)$ with $K(*) \to 0$ as the distance between $x_j$ and $x_i$ increases. From this, you can use the standard IRLS algorithm to get the solution, by applying it in turn to each $x$.

This is of course computationally inefficient, since you're applying an $O(N^2)$ algorithm to $N$ data points, making it $O(N^3)$ in all. This is the problem that the IVM is designed to solve.

• Hello, thanks a lot for your comment; the approach in the NIPS paper also may be computationally complex, since it involves a matrix inversion operation (which is $\mathcal{O}(N^{3})$). The choice of a kernel-based method was driven by the fact that such type of methods allow you to work on structured representations (e.g. trees, graphs, documents) but indeed also the Euclidean distance can be "kernelized". – Morra Cinese Jul 14 '13 at 18:46

I did some calculations by hand -- assuming that the bias term $b$ is the last value of the $\alpha$ vector, we get the following update equation:

$$\alpha^{(k + 1)} \leftarrow (\phi^{T} W \phi + \lambda R)^{-1} ((\phi^{T} W \phi + \lambda R) \alpha^{(k)} - (\phi^{T}(\mu - y) + \lambda \phi \alpha^{(k)}_{0})),$$

with $\phi = \left( \begin{array}{cc} K & 1 \end{array} \right)$, $R = \left( \begin{array}{cc} K & 0 \\ 0 & 0 \end{array} \right)$ and $\alpha_{0}$ is like $\alpha$, but with the bias term set to $0$ (so it doesn't get regularized). This follows from the Newton-Raphson update:

$$\alpha^{(k + 1)} \leftarrow \alpha^{(k)} - (\phi^{T} W \phi + \lambda R)^{-1} (\phi^{T} (\mu - y) + \lambda \phi \alpha^{(k)}_{0})$$

where the inverted term is the Hessian of the negative log-likelihood, and the other is its gradient, wrt. $\alpha$.

The negative log-likelihood has the following form:

$$NL(\alpha) = - \sum_{i = 1}^{n} (t_{i} \log \mu_{i} + (1 - t_{i}) \log (1 - \mu_{i})) + \alpha_{0}^{T} R \alpha_{0}.$$

Prediction boils down to the following (the bias term $b$ ended up in $\alpha$):

$$logit(x) = \sum_{i = 1}^{n} K(x, x_{i}) \alpha_{i} + \alpha_{n + 1};$$ $$\mu(x) = \frac{\exp(logit(x))}{1 + \exp(logit(x))}.$$

Please tell me if you spot any error. I hope my last weekend could be useful to someone else too. Thank you Hong and Dougal again for your time (the hint about using solve() is being really useful =) I would overload you in +1's if I could.)

There is a slightly better way of implementing IRWLS for kernel logistic regression that doesn't use the (weighted) normal equations and tends to have slightly better numerical properties. It basically fits the KLR model as a sequence of weighted Least-Squares Support Vector Machines (but with the logistic link function). There is a full derivation of this procedure in my paper (1) on approximate leave-one-out cross-validation for kernel logistic regression. If you use MATLAB, there is a toolbox that implements this approach.

HTH

(1) G. C. Cawley and N. L. C. Talbot, Efficient approximate leave-one-out cross-validation for kernel logistic regression, Machine Learning, vol, 71, no. 2-3, pp. 243--264, June 2008. (doi,pre-print)