Introduction to Neural Networks
U. Minn. Psy 5038

Linear discriminant

Initialize

Read in Statistical Add-in packages:

Off[General::"spell1"];
<<"MultivariateStatistics`";<<"ComputationalGeometry`"

Discriminant functions

Let's build our geometric intuitions of what a simple perceptron unit does by viewing it from a more formal point of view. Perceptron learning is an example of nonparametric statistical learning, because it doesn't require knowledge of the underlying probability distributions generating the data (such distributions are characterized by a relatively small number of "parameters", such as the mean and variance of a Gaussian distribution). Of course, how well it does will depend on the generative structure of the data. Much of the material below is covered in Duda and Hart (1978).

Linear discriminant functions: Two category case

A discriminant function, g(x) divides input space into two category regions depending on whether g(x)>0 or g(x)<0. (We've switched notation, x=f). The linear case corresponds to the simple perceptron unit we studied earlier:

g(x)=w.x+(Subscript[w, 0])

where w is the weight vector and (Subscript[w, 0]) is the (scalar) threshold (sometimes called bias, although this "bias" has nothing to do with statistical "bias"). 
Discriminant functions can be generalized, for example to quadratic decision surfaces:

g(x)=(Subscript[w, 0])+(∑_(i=1))(Subscript[w, i])(Subscript[x, i])+(∑_(i=1))(∑_(j=1))(Subscript[w, ij])(Subscript[x, i])(Subscript[x, j])

where x = {(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])...}. We've seen how g(x)=0 defines a decision surface which in the linear case is  a hyperplane. 
Suppose (Subscript[x , 1])and (Subscript[x, 2]) are vectors, with endpoints sitting on the hyperplane, then their difference is a vector lying in the hyperplane

w.(Subscript[x, 1])+(Subscript[w, 0])=w.(Subscript[x, 2])+(Subscript[w, 0])
			w.((Subscript[x, 1])-(Subscript[x, 2]))=0

so the weight vector w is normal to any vector lying in the hyperplane. Thus w determines how the plane is oriented. The normal vector w points into the region for which g(x)>0, and -w points into the region for which g(x)<0.

Let x be a point on the hyperplane. If we project x onto the normalized weight vector x.w/|w|, we have the normal distance of the hyperplane from the origin equal to:

w.x/|w|=-(Subscript[w, 0])/|w|

Thus, the threshold determines the position of the hyperplane.

One can also show that the normal distance of x to the hyperplane is given by:

g(x)/|w|

So we've seen that: 1) disriminant function divides the input space by a hyperplane decision surface; 2) The orientation of the surface is determined by the weight vector w; 3) the location is determined by the threshold (Subscript[w, 0]); 4) the discriminant function gives a measure of how far an input vector is from the hyperplane.
The figure summarizes the basic properties of the linear discriminant.

Clear[w1,w2,w,w0];
Manipulate[
(* x0={2,1};*)
w={w1,w2};wn=w/Norm[w];
g[x_]:={w1,w2}.x+w0;
gg=Plot[Tooltip[x2 /. Solve[{w1,w2}.{x1,x2}+w0==0,x2],"discriminant"],{x1,-1,3}];
ggg=Graphics[g[Dynamic[MousePosition["Graphics"]]]];
Show[{gg,Graphics[Inset["g[x]=",{1.6,2}]],Graphics[Inset[ToString[g[x0]],{2,2}]],Graphics[{Tooltip[Arrow[{{0,0},w}],"w"],Tooltip[Arrow[{{0,0},(-w0/Norm[w])*wn}],"-(Subscript[w, 0])/|w|"],Tooltip[{Arrow[{{0,0},x0}]},"x"],Tooltip[{Arrow[{x0,x0-wn*g[x0]/Norm[w]}]},"g(x)/|w|"]}]},PlotRange->{{-1,3},{-1,3}},AxesOrigin->{0,0},Axes->True,AspectRatio->1,ImageSize->Small],{{w0,-2.5},-6,3},{{w1,1},0,3},{{w2,2},0,3},{{x0,{2,1}},Locator}]

Multiple classes

Suppose there are c classes. There are a number of ways to define multiple class discriminant rules.  One way that avoids undefined regions is:

(Subscript[g, i])(x)=(Subscript[w, i]).x+(Subscript[w, i0]), i=1,...,c

Assign x to the ith class if: (Subscript[g, i])(x)>(Subscript[g, j])(x) for all j≠i.

(Adapted from Duda & Hart, 1973)

It can be shown that this classifier partitions the input space into simply connected convex regions. This means that if you connect any two feature vectors belonging to the same class by a line, all points on the line are in the same class. Thus this linear classifier won't be able to handle problems for which there are disconnected clusters of features that all belong to the same class. Also, from a probabilistic perspective, if the underlying generative probability model for a given class has multiple peaks, this linear classifier won't do a good job either.

Task-dependent Dimensionality reduction

Motivation

Later, when we consider the problem of dimensionality reduction, we will take another look at hyperplanes. But here the idea will be to find hyperplanes onto which we can project our input data, and from there divide up the hyperplane into decision regions. The idea is that the original input space may be impractically huge, but if we can find a subspace (hyperplane) that preserves the distinctions between categories as well as possible, we can make our decisions in smaller space. We will derive the Fisher linear "discriminant". 
This is closely related to the psychology idea of finding "distinctive" features. E.g. consider bird identification. If I want to discriminate cardinals from other birds in my backyard, I can make use of the fact that (males) cardinals may be the only birds that are red. So even tho' the image of a bird can have lots of dimensions, if I project the image on to the "red" axis, I can do fairly well with just one number. How about male vs. female human faces?

Fisher's linear "discriminant"

Generative model: two nearby gaussian classes

Define two bivariate base distributions

(ar = {{1, 0.99}, {0.99, 1}};
ndista = MultinormalDistribution[{0, -1}, ar];)
(br = {{1, .9}, {.9, 2}};
ndistb = MultinormalDistribution[{0, 1}, br];)

Find the expression for the probability distribution function of ndista

pdf = PDF[ndista, {x1, x2}]

Use Mean[ ] and CovarianceMatrix[ndista] to verify the population mean and the covariance matrix of ndistb

Mean[ndistb]

Covariance[ndista]

Try different covariant matrices. Should they be symmetric? Constraints on the determinant of ar, br?
Make a contour  plot of the PDF ndista

pdfa=PDF[ndista,{x1,x2}];
ContourPlot[pdfa,{x1,-3,3},{x2,-3,3},PlotPoints->64,PlotRange->All]

nsamples=500;
a=Table[Random[ndista],{nsamples}];
ga=ListPlot[a,PlotRange->{{-8,8},{-8,8}},AspectRatio->1,PlotStyle->Hue[0.2`],DisplayFunction->Identity];
b=Table[Random[ndistb],{nsamples}];
gb=ListPlot[b,PlotRange->{{-8,8},{-8,8}},AspectRatio->1,PlotStyle->Hue[0.6`],DisplayFunction->Identity];
Show[ga,gb,DisplayFunction->$DisplayFunction]

Use Mean[ ] to find the sample mean of  b. Whats is the sample covariance of b?

Mean[b]

Covariance[b]

Try out different projections of the data by varying the slope (m) of the discriminant line

m=-2/3;
wnvec={1,m}/Sqrt[1+m^2];

{{x,y}}.{n1,n2}
Map[#1*{n1,n2} &,{{x,y}}.{n1,n2}]

aproj=(#1 wnvec&)/@(a.wnvec);
gaproj=ListPlot[aproj,AspectRatio->1,PlotStyle->Hue[0.3`],DisplayFunction->Identity];
bproj=(#1 wnvec&)/@(b.wnvec);
gbproj=ListPlot[bproj,AspectRatio->1,PlotStyle->Hue[0.7`],DisplayFunction->Identity];
Show[ga,gb,gaproj,gbproj,DisplayFunction->$DisplayFunction,AspectRatio->Automatic]

By trial and error, find a value of m that separates the classes well along the projection line. Plot out the marginal distributions relative to this line.

Calculate the "signal-to-noise" ratio along the projection line.  Do this by taking the difference between the means divided by the square root of the product of the standard deviations along the line.

Mean[aproj]

Theory for simple 2-class case

(see Duda and Hart for general case)

A measure of the separation between the projections is the difference between the means:

|w.((Subscript[m, a])-(Subscript[m, b]))|
		and 
		(Subscript[m, a])=(Underoverscript[(1/N)∑, (i=1), N])x, summed over the N  x's from class a
		(Subscript[m, b])=(1/M)∑x, summed over the M  x's from class b

where w (wnvec) is the unknown unit vector along the discriminant line .

In our case above, the vector difference between the means is:

Mean[a]-Mean[b]

and the difference between the means projected onto a discriminant line is:

wnvec.(Mean[a]-Mean[b])

To improve separation, we can't just scale w, because the noise scales too.

We'd like the difference between the means to be large relative to the variation for each class. We can define a measure of the scatter for the projected samples in say class a (a==1), by:

(∑_(y∈class a))(((y-(Subscript[(OverTilde[m]), a])))^2)

where (Subscript[(OverTilde[m]), a])is the sample mean of the points from class a projected onto discriminant line and y=w.x 
Or in terms of the Mathematic example:

Apply[Plus,aproj-wnvec.Mean[a]];

The total scatter S is defined by the sum of the scatters for both classes (a and b).

S=(∑_(y∈class a))(((y-(Subscript[(OverTilde[m]), a])))^2)+(∑_(y∈class b))(((y-(Subscript[(OverTilde[m]), b])))^2)

If we divide the above number by the total number of points, we have an estimate of the variance of the combined data along the projected axis.

We now have the basic ingredients behind intuition for the Fisher linear discriminant. We'd like to find that w for which J:

J(w)=((|(Subscript[(OverTilde[m]), a])-(Subscript[(OverTilde[m]), b])((|)^2))/S)=((|w.((Subscript[m, a])-(Subscript[m, b]))((|)^2))/S)

is biggest. We want to maximize the difference between the projected class means, while minimizing the dispersion of the data on the projected line.

One can show that S = (w^T).(Subscript[S, W]).w, where

(Subscript[S, W])is measure of within-class variation called the within-class scatter matrix:

(Subscript[S, W])=(∑_(i=1)^2)(∑_(x∈Class i))(x-(Subscript[m, i]))(((x-(Subscript[m, i])))^T)

For the numerator, a measure of between class variation is the between-class scatter matrix:

(Subscript[S, B])=((Subscript[m, 1])-(Subscript[m, 2])).((((Subscript[m, 1])-(Subscript[m, 2])))^T)

and the difference between the projected means can be show to be:

|(Subscript[(OverTilde[m]), a])-(Subscript[(OverTilde[m]), b])((|)^2)=(w^T).(Subscript[S, B]).w

Find w (corresponding to slope) to maximize the criterion function

J(w)=(((w^T).(Subscript[S, B]).w)/((w^T).(Subscript[S, W]).w))

Answer:

w=((Subscript[S, W])^(-1)).((Subscript[m, a])-(Subscript[m, b]))

Demo: Finding Fisher's linear discriminant

normalize[x_]:=x/Sqrt[x.x];

ma=Mean[a];
mb=Mean[b];

Sa=Sum[Outer[Times,a[[i]]-ma,a[[i]]-ma],{i,1,nsamples}];
Sb=Sum[Outer[Times,b[[i]]-mb,b[[i]]-mb],{i,1,nsamples}];
Sw=Sa+Sb;
wldf = normalize[Inverse[Sw].(ma-mb)]

aproj=(#1 wldf&)/@(a.wldf);
gaproj=ListPlot[aproj,AspectRatio->1,PlotStyle->Hue[0.3`],DisplayFunction->Identity];
bproj=(#1 wldf&)/@(b.wldf);
gbproj=ListPlot[bproj,AspectRatio->1,PlotStyle->Hue[0.7`],DisplayFunction->Identity];
Show[ga,gb,gaproj,gbproj,DisplayFunction->$DisplayFunction]

We started off with a 2-dimensional input problem and turned it into a 1-D problem. For the n-dimensional case, see Duda and Hart.

Compare with the principal component axes

c=Join[a,b];
ListPlot[c,PlotRange->{{-8,8},{-8,8}}]

g1=ListPlot[c,PlotRange->{{-8,8},{-8,8}},AspectRatio->1,DisplayFunction->Identity];

auto=Covariance[c]
eigvalues=Eigenvalues[auto]
eigauto=Eigenvectors[auto]

gPCA = Plot[{eigauto[[1,2]]/eigauto[[1,1]] x,
	eigauto[[2,2]]/eigauto[[2,1]] x},
		{x,-4,4}, AspectRatio->1,
		DisplayFunction->Identity,
		PlotStyle->{RGBColor[.2,0,1]}];

Show[g1,gPCA,DisplayFunction->$DisplayFunction]

How does the principal component (biggest variance) compare with the Fisher discriminant line?

References

Duda, R. O., & Hart, P. E. (1973). Pattern classification and scene  analysis . New York.: John Wiley & Sons.
Duda, R. O., Hart, P. E., & Stork, D. G. (2001). Pattern classification (2nd ed.). New York: Wiley.
(Amazon.com)
Geman, S., Bienenstock, E., & Doursat, R. (1992). Neural networks and the bias/variance dilemma. Neural Computation, 4(1), 1-58.
MacKay, D. J. C. (1992). Bayesian interpolation. Neural Computation, 4(3), 415-447.
Ripley, B. D. (1996). Pattern Recognition and Neural Networks. Cambridge, UK: Cambridge University Press.
Vapnik, V. N. (1995). The nature of statistical learning. New York: Springer-Verlag.
http://neuron.eng.wayne.edu/software.html