# Kernel Logistic Classifier

## Linear logistic classifier

Consider the multi-dimensional logistic classifier
\begin{equation}
P(c\mid{\bf x},\Theta) = \frac{e^{\alpha_c+{\bf\beta}_c^{T}{\bf x}}}
{\sum_{c'=1}^{C}e^{\alpha_{c'}+{\bf\beta}_{c'}^{T}{\bf x}}}\,,
\end{equation}
for feature vector ${\bf x}\in\mathbb{R}^{F}$.
We can, if we wish, notionally consider the prior of class $c$ to be $P(c\mid\Theta)\propto e^{\alpha_c}$, and
the class density of ${\bf x}$ to be
$p({\bf x}\mid c,\Theta)\propto e^{{\bf\beta}_c^{T}{\bf x}}$, although the latter assumption poses some normalisation issues.

Now consider supervised training data comprised of $N$ known class labels ${\bf C}=[c_1,c_2,\ldots,c_N]^{T}$
and feature (or design) matrix ${\bf X}=[{\bf x}_1,{\bf x}_2,\ldots,{\bf x}_N]^{T}$.
The discriminative likelihood is then given by
\begin{eqnarray}
P({\bf C}\mid{\bf X},\Theta) & = & \prod_{d=1}^{N}P(c_d\mid{\bf x}_d,\Theta)
~=~ \prod_{d=1}^{N}\prod_{c=1}^{C}P(c\mid{\bf x}_d,\Theta)^{\delta(c_d=c)}\,,
\end{eqnarray}
where $\delta(A)=1$ and $\delta(\neg A)=0$ if proposition $A$ is true.

For notational convenience, we let $z_{cd}\doteq\delta(c_d=c)$ and $\pi_{cd}\doteq P(c\mid{\bf x}_d,\Theta)$. 
Then the discriminative log-likelihood is just
\begin{eqnarray}
L(\Theta) & = & \ln P({\bf C}\mid{\bf X},\Theta)
~=~\sum_{d=1}^{N}\sum_{c=1}^{C}z_{cd}\ln\pi_{cd}
\nonumber\\
& = & \sum_{d=1}^{N}\sum_{c=1}^{C}z_{cd}(\alpha_c+{\bf\beta}_c^{T}{\bf x}_d)
-\sum_{d=1}^{N}\ln\sum_{c=1}^{C}e^{\alpha_c+{\bf\beta}_c^{T}{\bf x}_d}\,.
\end{eqnarray}

Since discriminatively trained models are prone to overfitting, it is usual to regularise the parameters.
In the case of ridge regression, there is some dispute whether to penalise just the feature weights, proportional 
$\|{\bf\beta}_c\|^2$, or to also penalise the bias $\alpha_c$, proportional to $\|{\bf\gamma}_c\|^2$,
where ${\bf\gamma}_c\doteq [\alpha_c]\oplus{\bf\beta}_c$ is the concatenate of all parameters for class $c$.

The former case is sometimes preferred on the basis that $\alpha_c$ controls the prior on class $c$, and probably shouldn't be constrained beyond what the data suggest.
In particular, if we do not regularise the class weights, then it can be shown that
\begin{eqnarray}
\frac{1}{N}\sum_{d=1}^{N}P(c\mid{\bf x}_d,\Theta) & = & \frac{N_c}{N}\,,
\end{eqnarray}
where $N_c$ is the number of training samples of class $c$.

Alternatively, it could be noted that the observed class proportions at best only approximate the true class priors, and at worst are artificially constrained (e.g. by balancing class sizes). Hence, we might prefer the latter case of regularising all parameters.

For now we consider the more general, quadratic penalty ${\bf\gamma}_c^{T}{\bf\Lambda}_c{\bf\gamma}_c$, which allows us to not only "turn off" regularisation of some parameters, but additionally to handle differently-scaled features and correlations between features. We define the modified feature vector $\tilde{\bf x}\doteq[1]\oplus{\bf x}$, such that
\begin{eqnarray}
\pi_{cd} & = & \frac{e^{{\bf\gamma}_c^{T}\tilde{\bf x}_d}}
{\sum_{c'=1}^{C}e^{{\bf\gamma}_{c'}^{T}\tilde{\bf x}_d}}\,.
\end{eqnarray}
The ridge-regularised discriminative log-likelihood is then
\begin{eqnarray}
\tilde{L}(\Theta) & = & 
\sum_{d=1}^{N}\sum_{c=1}^{C}z_{cd}{\bf\gamma}_c^{T}\tilde{\bf x}_d
-\sum_{d=1}^{N}\ln\sum_{c=1}^{C}e^{{\bf\gamma}_c^{T}\tilde{\bf x}_d}
-\frac{1}{2}\sum_{c=1}^{C}{\bf\gamma}_c^T{\bf\Lambda}_c\mathbf{\gamma}_c
\,.
\end{eqnarray}

It can be shown that its class-specific gradient vector is given by
\begin{eqnarray}
{\bf\nabla}_{c}\tilde{L} & = & \frac{\partial\tilde{L}}{\partial{\bf\gamma}_c}
~=~\sum_{d=1}^{N}z_{cd}\tilde{\bf x}_d-\sum_{d=1}^{N}\pi_{cd}\tilde{\bf x}_d-{\bf\Lambda}_c{\bf\gamma}_c\,,
\end{eqnarray}
and the class-specific Hessian matrix is given by
\begin{eqnarray}
{\bf\nabla}_{c}^{T}{\bf\nabla}_{c}\tilde{L} & = & 
\frac{\partial^2\tilde{L}}{\partial{\bf\gamma}_c^{T}\partial{\bf\gamma}_c}
~=~-\sum_{d=1}^{N}\pi_{cd}(1-\pi_{cd})\tilde{\bf x}_d\tilde{\bf x}_d^{T}-{\bf\Lambda}_c\,.
\end{eqnarray}
Note that, for simplicity, we are going to ignore the explicit cross-class dependencies
${\bf\nabla}_{c'}^{T}{\bf\nabla}_{c}\tilde{L}$.

In matrix notation, let ${\bf z}_c=[z_{c1},\ldots,z_{cN}]^T$, ${\bf\pi}_c=[\pi_{c1},\ldots,\pi_{cN}]^T$,
and $\tilde{\bf X}={\bf 1}\oplus{\bf X}$. Then the gradient becomes
\begin{eqnarray}
{\bf\nabla}_{c}\tilde{L} & = &
\tilde{\bf X}^{T}({\bf z}_c-{\bf\pi}_c)-{\bf\Lambda}_c{\bf\gamma}_c\,.
\end{eqnarray}
Similarly, define $w_{cd}\doteq\pi_{cd}(1-\pi_{cd})$, and let ${\bf w}_c=[w_{c1},\ldots,w_{cN}]^T$
and ${\bf W}_c={\tt diag}\{{\bf w}_c\}$. Then the Hessian becomes
\begin{eqnarray}
{\bf\nabla}_{c}^{T}{\bf\nabla}_{c}\tilde{L} & = & 
-\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}-{\bf\Lambda}_c\,.
\end{eqnarray}

Now, the class-specific update for parameters ${\bf\gamma}_c$ takes the form of a single iteration of the Newton-Raphson method, namely
\begin{eqnarray}
{\bf\gamma}'_c & = & {\bf\gamma}_c - \left[{\bf\nabla}_{c}^{T}{\bf\nabla}_{c}\tilde{L}\right]^{-1}
{\bf\nabla}_{c}\tilde{L}
\nonumber\\
& = & {\bf\gamma}_c+\left[\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}+{\bf\Lambda}_c\right]^{-1}
\left[\tilde{\bf X}^{T}({\bf z}_c-{\bf\pi}_c)-{\bf\Lambda}_c{\bf\gamma}_c\right]
\nonumber\\
& = & \left[\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}+{\bf\Lambda}_c\right]^{-1}
\left[\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}{\bf\gamma}_c+{\bf\Lambda}_c{\bf\gamma}_c+
\tilde{\bf X}^{T}({\bf z}_c-{\bf\pi}_c)-{\bf\Lambda\gamma}_c\right]
\nonumber\\
& = & \left[\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}+{\bf\Lambda}_c\right]^{-1}\tilde{\bf X}^T{\bf W}_c
\left[\tilde{\bf X}{\bf\gamma}_c+{\bf W}_{c}^{-1}({\bf z}_c-{\bf\pi}_c)\right]\,.
\end{eqnarray}

## Dual optimisation

Observe from above that ${\bf z}_c-{\bf\pi}_c$ represents a vector of prediction errors for class $c$. 
Hence, if we define the weighted errors
${\bf e}_c\doteq{\bf W}_c^{-1}({\bf z}_c-{\bf\pi}_c)$
and the linear system
\begin{equation}
{\bf y}_c = \tilde{\bf X}{\bf\gamma}_c+{\bf e}_c\,,
\end{equation}
then we note that the class-specific parameter update becomes
\begin{eqnarray}
\left[\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}+{\bf\Lambda}_c\right]{\bf\gamma}'_c
& = & \tilde{\bf X}^T{\bf W}_c{\bf y}_c\,.
\end{eqnarray}
It is of interest that this update corresponds to a regularised form of the iteratively reweighted least-squares (IRLS) algorithm, except that here ${\bf y}_c$ itself varies with each iteration.

To motivate this observation, note that if we knew ${\bf y}_c$ in advance, then we could simply find the optimal
parameters ${\bf\gamma}_c$ via a weighted least-squares (WLS)
minimisation of the square error $\|{\bf W}_c{\bf e}_c\|^2$. However, we instead must obtain ${\bf y}_c$ via the following steps:
1. Choose initial parameters, ${\bf\gamma}_c$, for all class $c=1,2,\ldots,C$.
2. Compute the linear projection, $\tilde{\bf X}{\bf\gamma}_c$.
3. Compute the posterior probabilities, ${\bf\pi}_c$.
4. Compute the weighted prediction errors, ${\bf e}_c$.
5. Compute the 'observations', ${\bf y}_c$.

At this juncture, we may now find update parameters ${\bf\gamma}'_c$ that minimise the
square error $\|{\bf e}'_c\|^2=\|{\bf y}_c-\tilde{\bf X}^T{\bf\gamma}'_c\|^2$, satisfying the system
${\bf y}_c = \tilde{\bf X}{\bf\gamma}'_c+{\bf e}'_c$. This gives rise to the IRLS update above.

In conclusion, Newton-Raphson maximisation of the discriminative log-likelihood (the primal problem)
corresponds to IRLS minimisation of the (weighted) prediction error (the dual problem).

## Representer theorem

Now, a consequence of this duality is that we may apply the representer theorem. Put simply, the ridge-penalised function $f_c(\cdot)$ that minimises the square error $\sum_{d=1}^{N}\|y_{cd}-f_c(\tilde{\bf x}_d)\|^2$
satisfies $f_c(\tilde{\bf x})=\sum_{d=1}^{N}\omega_{cd}k(\tilde{\bf x}_d,\tilde{\bf x})$ for some positive-definite kernel function $k(\cdot,\cdot)$. In other words, the least-squares function interpolates over the known data points.

Hence, since the parameter update chooses ${\bf\gamma}'_c$ to minimise the square error
$\|{\bf y}_c-\tilde{\bf X}{\bf\gamma}'_c\|^2$, then we may take 
$f_c(\tilde{\bf x})\doteq{\bf\gamma}_c^{'T}\tilde{\bf x}$.

Let us now consider the scalar-product kernel $k({\bf x},{\bf y})\doteq {\bf x}^{T}{\bf y}$.
Then it follows that
\begin{equation}
{\bf\gamma}_c^{'T}\tilde{\bf x} = \sum_{d=1}^{N}\omega'_{cd}\tilde{\bf x}_d^{T}\tilde{\bf x}
= \left({\bf\omega}_c^{'T}\tilde{\bf X}\right)\tilde{\bf x}
\Rightarrow {\bf\gamma}'_c = \tilde{\bf X}^{T}{\bf\omega}'_c\,.
\end{equation}
Substituting this representation (for both the old and new parameter estimates) back into the parameter update then gives
\begin{eqnarray}
{\bf y}_c & = & \tilde{\bf X}\tilde{\bf X}^{T}{\bf\omega}_c+{\bf W}_{c}^{-1}({\bf z}_c-{\bf\pi}_c)\,,
\nonumber\\
\left[\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}+{\bf\Lambda}_c\right]\tilde{\bf X}^{T}{\bf\omega}'_c
&=&\tilde{\bf X}^T{\bf W}_c{\bf y}_c
\nonumber\\
\Rightarrow
\left[\tilde{\bf X}\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}\tilde{\bf X}^{T}
+\tilde{\bf X}{\bf\Lambda}_c\tilde{\bf X}^{T}\right]{\bf\omega}'_c
& = &\tilde{\bf X}\tilde{\bf X}^T{\bf W}_c{\bf y}_c\,.
\end{eqnarray}

In order to simplify these expressions further, we note that the general regulariser ${\bf\Lambda}_c$ now gives us a problem.
One way to avoid the problem is to reverse our generalisation, and take ${\bf\Lambda}_c=\lambda{\bf I}$ as per usual.
Hence, we obtain
\begin{eqnarray}
\left[\tilde{\bf X}\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}\tilde{\bf X}^{T}
+\lambda\tilde{\bf X}\tilde{\bf X}^{T}\right]{\bf\omega}'_c
&=&\tilde{\bf X}\tilde{\bf X}^T{\bf W}_c{\bf y}_c\,.
\end{eqnarray}
Note that this simplifying assumption now regularises both the class weights $\alpha_c$ and the feature weights
${\bf\beta}_c$. With some effort, we can avoid regularisation of $\alpha_c$ by instead choosing
${\bf\Lambda}_c\doteq\lambda\,{\tt diag}\{0,1,\ldots,1\}$, whereupon
\begin{eqnarray}
\left[\tilde{\bf X}\tilde{\bf X}^T{\bf W}_c\tilde{\bf X}\tilde{\bf X}^{T}
+\lambda{\bf X}{\bf X}^{T}\right]{\bf\omega}'_c
&=&\tilde{\bf X}\tilde{\bf X}^T{\bf W}_c{\bf y}_c\,,
\end{eqnarray}
since $\tilde{\bf X}={\bf 1}\oplus{\bf X}$.
We shall return to this point later.

## Kernel trick

Now, we define the kernel matrix $\tilde{\bf K}\doteq\tilde{\bf X}\tilde{\bf X}^{T}$, 
where $\tilde{K}_{ij}=\tilde{\bf x}_i^{T}\tilde{\bf x}_j=k(\tilde{\bf x}_i,\tilde{\bf x}_j)$.
Consequently, we obtain
\begin{eqnarray}
{\bf y}_c & = &\tilde{\bf K}{\bf\omega}_c+{\bf W}_{c}^{-1}({\bf z}_c-{\bf\pi}_c)\,,
\nonumber\\
\left[\tilde{\bf K}{\bf W}_c\tilde{\bf K}
+\lambda\tilde{\bf K}\right]{\bf\omega}'_c
& = & \tilde{\bf K}{\bf W}_c{\bf y}_c
\nonumber\\
\Rightarrow \left[{\bf W}_c\tilde{\bf K}+\lambda{\bf I}\right]{\bf\omega}'_c & = & {\bf W}_c{\bf y}_c
={\bf W}_c\tilde{\bf K}{\bf\omega}_c+({\bf z}_c-{\bf\pi}_c)\,,
\end{eqnarray}
since $\tilde{\bf K}$ is invertible, from the definition that the kernel function $k(\cdot,\cdot)$
is positive-definite.
Note that we also assumed previously that the diagonal matrix
${\bf W}_c$ is invertible, hence we could simplify further to
\begin{eqnarray}
\left[\tilde{\bf K}+\lambda{\bf W}_c^{-1}\right]{\bf\omega}'_c & = &
\tilde{\bf K}{\bf\omega}_c+{\bf W}_{c}^{-1}({\bf z}_c-{\bf\pi}_c)\,.
\end{eqnarray}
We should note at this point that we have now traded an $F\times F$ inversion problem (in the parameter space) for an $N\times N$ one (in the data space). 

We have also, as noted earlier, lost the propery of being able to *not* penalise the $\alpha_c$ bias (or class proportion) parameter.
To see what form the update would take if we did not regularise $\alpha_c$, we instead take the alternative
regularisation
${\bf\Lambda}_c\doteq\lambda\,{\tt diag}\{0,1,\ldots,1\}$. Then, from the derivation above, 
we obtain the modified version
\begin{eqnarray}
\left[{\bf W}_c\tilde{\bf K}+\lambda\tilde{\bf K}^{-1}{\bf K}\right]{\bf\omega}'_c
&=&{\bf W}_c\tilde{\bf K}{\bf \omega}_c+({\bf z}_c-\mathbb{\pi}_c)\,,
\end{eqnarray}
where $\tilde{\bf K}={\bf K}+{\bf 1}{\bf 1}^T$.
We deduce from the Sherman-Morrison formula that $\tilde{\bf K}^{-1}{\bf K}={\bf I}-{\bf v}{\bf 1}^{T}$,
where ${\bf v}={\bf u}/(1+{\bf 1}^{T}{\bf u})$ and ${\bf K}{\bf u}={\bf 1}$.
Hence, the full version of the class-specific parameter update is
\begin{eqnarray}
\left[{\bf W}_c({\bf K}+{\bf 1}{\bf 1}^T)+\lambda({\bf I}-{\bf v}{\bf 1}^T)\right]{\bf\omega}'_c
&=&{\bf W}_c({\bf K}+{\bf 1}{\bf 1}^T){\bf \omega}_c+({\bf z}_c-\mathbb{\pi}_c)\,,
\end{eqnarray}
or
\begin{eqnarray}
\left[{\bf W}_c{\bf K}+\lambda{\bf I}+({\bf w}_c-\lambda{\bf v}){\bf 1}^T\right]{\bf\omega}'_c
&=&({\bf W}_c{\bf K}+{\bf w}_c{\bf 1}^T){\bf \omega}_c+({\bf z}_c-\mathbb{\pi}_c)\,,
\end{eqnarray}
since ${\bf W}_c{\bf 1}={\bf w}_c$.

Finally, we recall that ${\bf\gamma}_c=\tilde{\bf X}^{T}{\bf\omega}_c$, such that our original parameters are recovered as
\begin{eqnarray}
\alpha_c~=~{\bf 1}^{T}{\bf\omega}_c\,, && {\bf\beta}_c~=~{\bf X}^{T}{\bf\omega}_c\,.
\end{eqnarray}

679844daf65d963ede4777262fafa50ab62864af
ipynb
Jupyter Notebook
LogisticModels/notebooks/kernel-logistic-classifier.ipynb
gaj67/gaj-data-science
aadcf6ee2cd00606563f213167c2eeeb42430c59
MIT
LogisticModels/notebooks/kernel-logistic-classifier.ipynb
gaj67/gaj-data-science
aadcf6ee2cd00606563f213167c2eeeb42430c59
MIT
LogisticModels/notebooks/kernel-logistic-classifier.ipynb
gaj67/gaj-data-science
aadcf6ee2cd00606563f213167c2eeeb42430c59
MIT
Qwen/Qwen-72B
1. YES
2. YES