# Logistic Classification¶


This tour details the logistic classification method (for 2 classes and multi-classes).

Warning: Logisitic classification is actually called "logistic regression" in the literature, but it is in fact a classification method.

We recommend that after doing this Numerical Tours, you apply it to your own data, for instance using a dataset from LibSVM.

Disclaimer: these machine learning tours are intended to be overly-simplistic implementations and applications of baseline machine learning methods. For more advanced uses and implementations, we recommend to use a state-of-the-art library, the most well known being Scikit-Learn

In [4]:
using PyPlot
using NtToolBox


We define a few helpers.

In [5]:
#SetAR = @(ar)set(gca, 'PlotBoxAspectRatio', [1 ar 1], 'FontSize', 10);
Xm = X -> X - repeat(mean(X,1), outer=(size(X,1), 1))
Cov = X -> Xm(X)'*Xm(X)
dotp = (u,v) ->sum(u[:].*v[:]);


## Two Classes Logistic Classification¶

Logistic classification is, with support vector machine (SVM), the baseline method to perform classification. Its main advantage over SVM is that is is a smooth minimization problem, and that it also output class probabity, offering a probabilistic interpretation of the classification.

To understand the behavior of the method, we generate synthetic data distributed according to a mixture of Gaussian with an overlap governed by an offset $\omega$. Here classes indexes are set to $y_i \in \{-1,1\}$ to simplify the equations.

In [6]:
n = 1000 # number of sample
p = 2 # dimensionality
omega = [1 .5]*5 # offset
X = [randn(div(n,2),2); randn(div(n,2),2)+ones(div(n,2),1)*omega]
y = [ones(div(n,2),1);-ones(div(n,2),1)];

UndefVarError: div not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Plot the classes.

In [7]:
plot_multiclasses(X,y,disp_dim=2,ms=5);

UndefVarError: plot_multiclasses not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Logistic classification minimize a logistic loss in place of the usual $\ell^2$ loss for regression $$\umin{w} E(w) \eqdef \frac{1}{n} \sum_{i=1}^n L(\dotp{x_i}{w},y_i)$$ where the logistic loss reads $$L( s,y ) \eqdef \log( 1+\exp(-sy) )$$ This corresponds to a smooth convex minimization. If $X$ is injective, this is also strictly convex, hence it has a single global minimum.

Compare the binary (ideal) 0-1 loss, the logistic loss and the <https://en.wikipedia.org/wiki/Hinge_loss hinge loss> (the one used for SVM).

In [8]:
t = linspace(-3,3,255);
plot(t, log(1+exp(t)), c="red",lw=2, label="Logistic")
plot(t, t.>0, c="orange",lw=2, label="Hinge")
plot(t,max(t,0), c="blue", label="Binary")
axis("tight"); legend(loc="best");


This can be interpreted as a <https://en.wikipedia.org/wiki/Maximum_likelihood_estimation maximum likelihood estimator> when one models the probability of belonging to the two classes for sample $x_i$ as $$h(x_i) \eqdef (\th(x_i),1-\th(x_i)) \qwhereq \th(s) \eqdef \frac{e^{s}}{1+e^s} = (1+e^{-s})^{-1}$$

Re-writting the energy to minimize $$E(w) = \Ll(X w,y) \qwhereq \Ll(s,y)= \frac{1}{n} \sum_i L(s_i,y_i),$$ its gradient reads $$\nabla E(w) = X^\top \nabla \Ll(X w,y) \qwhereq \nabla \Ll(s,y) = \frac{y}{n} \odot \th(-y \odot s),$$ where $\odot$ is the pointwise multiplication operator, i.e. |.*| in Matlab.

Define the energies.

In [9]:
L = (s,y) -> 1/n * sum( log( 1 + exp(-s.*y) ) );
E = (w,X,y) -> L(X*w,y);


In [10]:
theta = v -> 1 ./ (1+exp(-v));
nablaL = (s,r) -> - 1/n * y.* theta(-s.*y);
nablaE = (w,X,y) -> X'*nablaL(X*w,y);


Important: in order to improve performance, it is important (especially in low dimension $p$) to add a constant bias term $w_{p+1} \in \RR$, and replace $\dotp{x_i}{w}$ by $\dotp{x_i}{w} + w_{p+1}$. This is equivalently achieved by adding an extra $(p+1)^{\text{th}}$ dimension equal to 1 to each $x_i$, which we do using a convenient macro.

In [11]:
AddBias = X -> [X ones(size(X,1),1)];


With this added bias term, once $w_{\ell=0} \in \RR^{p+1}$ initialized (for instance at $0_{p+1}$),

In [12]:
w = zeros(p+1,1);


one step of gradient descent reads $$w_{\ell+1} = w_\ell - \tau_\ell \nabla E(w_\ell).$$

In [13]:
tau = .8; # here we are using a fixed tau
w = w - tau * nablaE(w,AddBias(X),y);

UndefVarError: X not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Exercise 1

Implement a gradient descent $$w_{\ell+1} = w_\ell - \tau_\ell \nabla E(w_\ell).$$ Monitor the energy decay. Test different step size, and compare with the theory (in particular plot in log domain to illustrate the linear rate). etAR(1); etAR(1);

In [14]:
include("NtSolutions/ml_3_classification/exo1.jl");

LoadError: UndefVarError: X not defined

in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/ml_3_classification/exo1.jl:5 [inlined]
in anonymous at ./<missing>:?
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
In [15]:
# Insert your code here.


Generate a 2D grid of points.

In [16]:
q = 201; tx = linspace(minimum(X[:,1]),maximum(X[:,1]),q); ty = linspace(minimum(X[:,2]),maximum(X[:,2]),q)
B,A = meshgrid( ty,tx )
G = [A[:] B[:]];

UndefVarError: X not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Evaluate class probability associated to weight vectors on this grid.

In [17]:
Theta = theta(AddBias(G)*w);
Theta = reshape(Theta, q, q);

UndefVarError: G not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Display the data overlaid on top of the classification probability, this highlight the separating hyperplane $\enscond{x}{\dotp{w}{x}=0}$.

In [18]:
imshow(Theta'[:,end:-1:1], extent=[minimum(tx), maximum(tx), minimum(ty), maximum(ty)]);
plot_multiclasses(X,y);

UndefVarError: tx not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Exercise 2

Test the influence of the separation offset $\omega$ on the result.

In [19]:
include("NtSolutions/ml_3_classification/exo2.jl");

LoadError: UndefVarError: div not defined

in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/ml_3_classification/exo2.jl:12 [inlined]
in anonymous at ./<missing>:?
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
In [20]:
# Insert your code here.


Exercise 3

Test logistic classification on a real life dataset. You can look at the Numerical Tour on stochastic gradient descent for an example. Split the data in training and testing to evaluate the classification performance, and check the impact of regularization.

In [21]:
include("NtSolutions/ml_3_classification/exo3.jl");

In [22]:
# Insert your code here.


## Kernelized Logistic Classification¶

Logistic classification tries to separate the classes using a linear separating hyperplane $\enscond{x}{\dotp{w}{x}=0}.$

In order to generate a non-linear descision boundary, one can replace the parametric linear model by a non-linear non-parametric model, thanks to kernelization. It is non-parametric in the sense that the number of parameter grows with the number $n$ of sample (while for the basic method, the number of parameter is $p$. This allows in particular to generate decision boundary of arbitrary complexity.

The downside is that the numerical complexity of the method grows (at least) quadratically with $n$.

The good news however is that thanks to the theory of reproducing kernel Hilbert spaces (RKHS), one can still compute this non-linear decision function using (almost) the same numerical algorithm.

Given a kernel $\kappa(x,z) \in \RR$ defined for $(x,z) \in \RR^p$, the kernelized method replace the linear decision functional $f(x) = \dotp{x}{w}$ by a sum of kernel centered on the samples $$f_h(x) = \sum_{i=1}^p h_i k(x_i,x)$$ where $h \in \RR^n$ is the unknown vector of weight to find.

When using the linear kernel $\kappa(x,y)=\dotp{x}{y}$, one retrieves the previously studied linear method.

Macro to compute pairwise squared Euclidean distance matrix.

In [23]:
distmat = (X,Z) -> broadcast(+,sum(X'.*X',1)',sum(Z'.*Z',1))-2*(X*Z');


The gaussian kernel is the most well known and used kernel $$\kappa(x,y) \eqdef e^{-\frac{\norm{x-y}^2}{2\sigma^2}} .$$ The bandwidth parameter $\si>0$ is crucial and controls the locality of the model. It is typically tuned through cross validation.

In [24]:
kappa = (X,Z,sigma) -> exp( -distmat(X,Z)/(2*sigma^2) );


We generate synthetic data in 2-D which are not separable by an hyperplane.

In [25]:
n = 1000; p = 2
t = 2*pi*rand(div(n,2),1)
R = 2.5
r = R*(1 + .2*rand(div(n,2),1)); # radius
X1 = [cos(t).*r sin(t).*r]
X = [randn(div(n,2),2); X1]
y = [ones(div(n,2),1);-ones(div(n,2),1)];

UndefVarError: div not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Display the classes.

In [26]:
plot_multiclasses(X,y,disp_dim=2;disp_legend=1);
axis("off");

UndefVarError: plot_multiclasses not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Once avaluated on grid points, the kernel define a matrix $$K = (\kappa(x_i,x_j))_{i,j=1}^n \in \RR^{n \times n}.$$

In [27]:
sigma = 1
K = kappa(X,X,sigma);

UndefVarError: X not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Valid kernels are those that gives rise to positive symmetric matrices $K$. The linear and Gaussian kernel are valid kernel functions. Other popular kernels include the polynomial kernel $\dotp{x}{y}^a$ for $a \geq 1$ and the Laplacian kernel $\exp( -\norm{x-y}^2/\si )$.

The kernelized Logistic minimization reads $$\umin{h} F(h) \eqdef \Ll(K h,y).$$

In [28]:
F = (h,K,y) -> L(K*h,y);
nablaF = (h,K,y) -> K'*nablaL(K*h,y);


This minimization can be related to an infinite dimensional optimization problem where one minimizes directly over the function $f$. This is shown to be equivalent to the above finite-dimenisonal optimization problem thanks to the theory of RKHS.

Exercise 4

Implement a gradient descent to minimize $F(h)$. Monitor the energy decay. Test different step size, and compare with the theory.

In [29]:
include("NtSolutions/ml_3_classification/exo4.jl");

LoadError: UndefVarError: K not defined

in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/ml_3_classification/exo4.jl:6 [inlined]
in anonymous at ./<missing>:?
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
In [30]:
# Insert your code here.


Once this optimal $h$ has been found, class probability at a point $x$ are obtained as $$(\th(f_h(x)), 1-\th(f_h(x))$$ where $f_h$ has been defined above.

We evaluate this classification probability on a grid.

In [31]:
q = 201
tmax = 3.5
t = linspace(-tmax,tmax,q)
B,A = meshgrid( t,t )
G = [A[:] B[:]]
Theta = reshape( theta(kappa(G,X,sigma)*h) , q,q);

UndefVarError: X not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Display the classification probability.

In [32]:
imshow(Theta'[:,end:-1:1], extent=[-tmax, tmax, -tmax, tmax], cmap = get_cmap("jet"));
plot_multiclasses(X,y,disp_legend = 0);
clim(0,1)
#caxis([0 1]);
axis("off");

UndefVarError: Theta not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Exercise 5

Display evolution of the classification probability with $\sigma$

In [33]:
include("NtSolutions/ml_3_classification/exo5.jl");

LoadError: UndefVarError: X not defined

in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/ml_3_classification/exo5.jl:7 [inlined]
in anonymous at ./<missing>:?
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
In [34]:
# Insert your code here.


Exercise 6

Separate the dataset into a training set and a testing set. Evaluate the classification performance for varying $\si$. Try to introduce regularization and minmize $$\umin{h} F(h) \eqdef \Ll(K h,y) + \la R(h)$$ where for instance $R=\norm{\cdot}_2^2$ or $R=\norm{\cdot}_1$.

In [35]:
include("NtSolutions/ml_3_classification/exo6.jl");

In [36]:
# Insert your code here.


## Multi-Classes Logistic Classification¶

The logistic classification method is extended to an arbitrary number $k$ of classes by considering a familly of weight vectors $w_\ell$_{\ell=1}^k, which are conveniently stored as columns of matrix $W \in \RR^{p \times k}$.

This allows to model probabilitically the belonging of a point $x \in \RR^p$ to a the classes using an exponential model $$h(x) = \pa{ \frac{ e^{-\dotp{x}{w_\ell}} }{ \sum_m e^{-\dotp{x}{w_m}} } }_\ell$$ This vector $h(x) \in [0,1]^k$ describes the probability of $x$ belonging to the different classes, and $\sum_\ell h(x)_\ell = 1$.

The computation of $w$ is obtained by solving a maximum likelihood estimator $$\umax{w \in \RR^k} \frac{1}{n} \sum_{i=1}^n \log( h(x_i)_{y_i} )$$ where we recall that $y_i \in \{1,\ldots,k\}$ is the class index of point $x_i$.

This is conveniently rewritten as $$\umin{w} \sum_i \text{LSE}( XW )_i - \dotp{XW}{D}$$ where $D \in \{0,1\}^{n \times k}$ is the binary class index matrices $$D_{i,\ell} = \choice{ 1 \qifq y_i=\ell, \\ 0 \quad \text{otherwise}. }$$ and LSE is the log-sum-exp operator $$\text{LSE}(S) = \log\pa{ \sum_\ell \exp(S_{i,\ell}) } \in \RR^n.$$

In [37]:
LSE = S -> log( sum(exp(S), 2) );


The computation of LSE is unstable for large value of $S_{i,\ell}$ (numerical overflow, producing NaN), but this can be fixed by substracting the largest element in each row, since $\text{LSE}(S+a)=\text{LSE}(S)+a$ if $a$ is constant along rows. This is the celebrated LSE trick.

In [38]:
max2 = S -> repeat(mapslices(maximum, S,2), outer=(1, size(S,2)))
LSE_stab = S -> LSE( S-max2(S) ) + mapslices(maximum, S,2);


The gradient of the LSE operator is the <https://en.wikipedia.org/wiki/Softmax_function soft-max operator> $$\nabla \text{LSE}(S) = \text{SM}(S) \eqdef \pa{ \frac{ e^{S_{i,\ell}} }{ \sum_m e^{S_{i,m}} } }$$

In [39]:
SM = S -> exp(S) ./ repeat( sum(exp(S),2), outer=(1, size(S,2)));


Similarely to the LSE, it needs to be stabilized.

In [40]:
SM_stab = S -> SM(S-max2(S));


We load a dataset of $n$ images of size $p = 8 \times 8$, representing digits from 0 to 9 (so there are $k=10$ classes).

Load the dataset and randomly permute it. Separate the features $X$ from the data $y$ to predict information.

In [41]:
A = readdlm("NtToolBox/src/data/digits.csv", ',')
A = A[randperm(size(A,1)),:]
X = A[:,1:end-1]; y = A[:,end];


$n$ is the number of samples, $p$ is the dimensionality of the features, $k$ the number of classes.

In [42]:
n,p = size(X);
CL = unique(y); # list of classes.
k = length(CL);


Display a few samples digits

In [43]:
q = 5
for i=1:k
I = find(y.==CL[i]);
for j=1:q
f = reshape(X[I[j],:], Int(sqrt(p)), Int(sqrt(p)))';
subplot(q,k, (j-1)*k+i );
imshow(-f, extent=[0,1,0,1], cmap=get_cmap("gray")); axis("image"); axis("off")
end
end


Display in 2D.

In [44]:
plot_multiclasses(X,y,disp_dim = 2, ms = 5, disp_legend = 1);

UndefVarError: plot_multiclasses not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Display in 3D.

In [45]:
plot_multiclasses(X,y,disp_dim=3);

UndefVarError: plot_multiclasses not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Compute the $D$ matrix.

In [46]:
D = Int.(repeat(CL[:]', outer=(n,1)) .== repeat(y, outer=(1,k))); #double


Define the energy $E(W)$.

In [47]:
E = W -> 1/n*( sum(LSE(X*W)) - dotp(X*W,D)  );


Define its gradients $$\nabla E(W) = \frac{1}{n} X^\top ( \text{SM}(X W) - D ).$$

In [48]:
nablaE = W -> 1/n * X'* ( SM_stab(X*W) -  D );


Exercise 7

Implement a gradient descent $$W_{\ell+1} = W_\ell - \tau_\ell \nabla E(W_\ell).$$ Monitor the energy decay.

In [49]:
include("NtSolutions/ml_3_classification/exo7.jl");

In [50]:
W = zeros(p,k)
Elist = []
tau = .01
niter = 500
for i=1:niter
W = W - tau * nablaE(W)
append!(Elist, E(W))
end
plot(1:niter, Elist, lw=2)
axis("tight");

In [51]:
# Insert your code here.


Generate a 2D grid of points over PCA space and map it to feature space.

In [52]:
U,D,V = svd(Xm(X),thin=true)
Z = Xm(X) * V
M = maximum(abs(Z[:]))
q = 201
t = linspace(-M,M,q)
B,A = meshgrid(t,t)
G = zeros(q*q,p)
G[:,1:2] = [A[:]; B[:]]
G = G*V' + repeat(mean(X,1), outer=(q*q, 1));


Evaluate class probability associated to weight vectors on this grid.

In [53]:
Theta = SM(G*W)
Theta = reshape(Theta, q, q, k);


Display each probablity map.

In [54]:
for i=1:k
subplot(3,4,i)
imshow(Theta[:,:,i]'[:,end:-1:1], extent = [0,1,0,1], cmap=get_cmap("jet"));
title(string("Class ", i))
axis("image"); axis("off")
end
tight_layout()


Build a single color image of this map.

In [55]:
col = [ [1 0 0]; [0 1 0]; [0 0 1]; [0 0 0]; [0 1 1]; [1 0 1]; [1 1 0];
[1 .5 .5]; [.5 1 .5]; [.5 .5 1]  ]'
R = zeros(q,q,3)
for i=1:k
for a=1:3
R[:,:,a] = R[:,:,a] + Theta[:,:,i] .* col[a,i]
end
end


Display.

In [56]:
imshow(permutedims(R, [2, 1, 3]), extent = [-tmax, tmax, -tmax, tmax]);
plot_multiclasses3(X,y,disp_dim=2,ms=3,disp_legend=1);
axis("off");

UndefVarError: plot_multiclasses3 not defined

in include_string(::String, ::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
include("NtSolutions/ml_3_classification/exo8.jl");

# Insert your code here.