Introduction to Neural Networks

Regression, least squares, gradient descent, Widrow-Hoff & backprop

Initialization

Off[SetDelayed::write]
Off[General::spell1]

Introduction

Last time

Focused on classification: 
	Turning linear networks into classifiers with a hard threshold. 
Perceptrons, perceptron learning rule.
Linear separability

Today

Linear regression & brain-style learning -- the Widrow-Hoff learning rule
Multiple layers of weights -- backprop learning, an example of non-linear regression

Linear regression and Widrow-Hoff learning

Introduction

We have been studying a linear matrix model of memory based on the storage of connection weights that follow a particular Hebbian rule. We have studied the "psychology" of some operations of linear algebra and have seen some interesting parallels to human memory, such as interference and pattern reconstruction.

We've learned that linear networks can be configured to be supervised or unsupervised learning devices. But linear mappings are severely limited in what they can compute. With a threshold non-linearity, the linear neuron model becomes a simple perceptron unit that makes decisions.

But we are now limited by the linearity of the decision surface, i.e. by the fact that it is a hyperplane. More complex networks can be built by adding layers making multiple layer networks (“multi-layer perceptrons”), but then we have another problem: how can we learn the weights in this more complicated setting?

Overview of our strategy

First we return to strictly linear networks (no threshold) which will allow us to introduce techniques (gradient descent learning in the context of finding the weights of a linear matrix transformation). We will learn concepts that will generalize to non-linear networks with multiple layers of weights. In particular, the technique of gradient descent will lead us to the Widrow-Hoff learning rule. By treating the linear case first, we will be able to see how the Widrow-Hoff learning rule relates to the classic problem of linear regression. 

Learning, memory and generalization can be viewed from the point of view of non-linear statistical regression. The idea is to treat learning as an attempt to fit past input/output associations into a model that can be used both for recall and for future generalization.

The problem of regression in statistics is: given a set of vector inputs {xi}, and a set of corresponding vector “target” outputs {ti}, one tries to find a transformation W that will map x→t as closely as possible over the data available. For this we need a model for W, and a measure of goodness of fit. At first, we will assume a linear model for W, so for the discrete case, W is a matrix. Our measure of goodness of fit will be the sum of the squared differences between predicted output and actual output. In the linear case, this is linear least squares regression. 

Then we’ll see how the Widrow-Hoff rule provides a clue how to generalize learning to non-linear, multiple layer networks with smooth non-linear squashing functions, with the same measure of goodness of fit. This is the error back-propagation algorithm.

In later lectures,  we will develop a more powerful notion of goodness of fit -- probability of parameters given experience with past inputs and outputs, for which the sum of squared differences is a special case.

More terminology regarding supervised learning

Consider supervised learning. We have a "training set" {xi,ti} representing input vectors xi, and target outputs ti. The training set "samples" the larger space of possible input/output pairs {x, t}.   We would like to learn a general mapping: T: x -> y in such a way that T (= W) is a good fit to the training data (i.e. y is close to t), and generalizes well to novel inputs. The set of target data is the feedback for the"teacher". 

In general, feedback can vary in the degree of precision it provides for learning. It can specify whether the mapping is correct or not. Or the feedback can provide information as to how far off the map T's prediction, i.e. how far is T[xi]) from ti ? 

After training, one could require that T always maps members of the training set to exactly the target members, and generalizes appropriately for other inputs. This means that the learning should be consistent. Interpolation between data points on a graph is an example of consistent learning. A specific example of interpolation would be drawing lines connecting data points on a graph. 

Or, we may require that the  T maps the original members of the training set to outputs y, that are close to the original targets t. This is called approximation learning.

gdata=ListPlot[data=Table[{x,x+4*RandomReal[]},{x,0,7}]];
approx = Fit[data,{1,x},x];
interp= Interpolation[data]; 
Show[{Plot[{approx, interp[x]},{x,0,7}],gdata},ImageSize->Small]

Linear regression is an example of approximation learning. Approximation doesn't necessarily exactly fit the data points, but the goal is that it should come close, but not at the expense of failing to generalize to new data. 

In addition to remembering, regression generalizes, i.e. provides the basis for prediction.  So novel input values get mapped to predicted outputs based on past “experience”. We’ve already seen how linear heteroassociative learning does this. Later we will look at how the type of generalization affects “bias” and “variance”.

Let's take a close look at linear regression, a fundamental case of approximation learning.

Least squares regression - linear models

The idea behind least squares regression is given a set of N training pairs {xi, ti} , where i runs from 1 to N, we would like to find a function that given an input x, the function approximates well the output, i.e.  W.x ~ t. If the function reproduces the association between input and output that it has seen before, this is like "remembering".

An example of linear regression

A fundamental assumption that affects the modeling and performance of any learning system which needs to generalize is that there is an underlying structure to the data--the relationship between associative pairs is not arbitrary.

When trying to understand how a relationship can be learned between a set of two patterns, it helps to have some understanding of the structure of the relationship. For example, one shouldn't try to fit a straight line to data when there is reason to believe that the underlying process produces data on a circle. We will study this more when we learn about the "bias/variance dilemma". So let us assume that the data have an underlying structure that we are going to try to discover or approximate using W.

We will study a simple "toy" problem that has the following very specific generative structure. The inputs are randomly located points on a 2D plane, and the outputs are heights above these points. The outputs lie approximately on a planar surface that runs through the origin and whose orientation is characterized by two parameters, a and b.

It may seem like overkill, but we are going to estimate W for the same set of data in four different ways. 

The first two are drawn from standard linear algebra (a least squares solution using transpose and inverse, and the second using the pseudoinverse).  The third introduces the method of "gradient descent" a common numerical technique which is used in many applications. We will use it later in several contexts. The fourth method is the most relevant to neural network theory--we estimate W using a biologically plausible learning rule (the Widrow-Hoff rule).

Generative model: Synthetic training pairs

Although we will illustrate our examples with small dimensions, everything we do generalizes to higher dimensional inputs and outputs, and in fact the demonstration code below will work with higher dimensional input/output vectors.

Let x1 and x2 be the (scalar) inputs, and y be the (scalar) output:

{x1,x2} -> y, 

		and assuming the mapping is a plane through the origin and additive noise

		y = a x1 + b x2 + noise

This is a many-to-one mapping. 

We want to learn the parameters a and b from the training pairs: {{x1, x2}, y}.
(These parameters could be thought of as neural weights.)

rsurface[a_,b_]:=N[Table[{x1=1 RandomReal[],x2=1 RandomReal[],a x1+b x2+0.5 RandomReal[]-0.25},{120}],2];
data=rsurface[2,3];
Graphics3D[Point/@data,AxesLabel->{"x1","x2","y"},Axes->True,Ticks->None,BoxRatios->{1,1,1},ImageSize->Small]
Dimensions[data]

You can also use: ListPointPlot3D[data]

yy represents the 1 dimensional y values, and xx the 2-dimensional input values x = {x1,x2}:

yy = data[[All, 3]];
xx = data[[All, 1 ;; 2]];

1. Least squares regression to find W

So let's assume we want to find a matrix W that will come close to reproducing values y, given inputs x. Of  course, because we generated the data, we know the underlying structure and what the matrix W should be. The correct answer should be a 1x2 matrix = {{2,3}}. But let's assume we don't know the answer, and want to discover the weights from yy and xx.

Imagine you have knobs that you can twiddle that let you adjust every element of the matrix W. For particular settings of these knobs, you can calculate W.(Subscript[x, i]) for any input/output pair. In general you won't be so lucky as to have (Subscript[y, i])=W.(Subscript[x, i]). So you take the opportunity to adjust the knobs to produce W' so that W'. (Subscript[x, i]) is closer to (Subscript[y, i]). One measure of "how close" is (((Norm[(Subscript[y, i])-W. (Subscript[x, i])]))^2) = |(Subscript[y, i])-W. (Subscript[x, i])((|)^2)= EuclideanDistance[(Subscript[y, i]),W. (Subscript[x, i])]^2. Of course, we usually have lots of data, so we could add up all the discrepancies between what W predicts, and what the "teacher" has provided in the answers in the training pair (i.e. (Subscript[y, i]) for all of the i's).
In least squares regression, we try to find the values of the matrix that will minimize e(W).

e(W) = (∑_(i=1)^N)((\[LeftBracketingBar](Subscript[y, i])-(Subscript[Wx, i])\[RightBracketingBar])^2)

In many cases our error function will be over a very high dimensional space determined by the number of elements in W. However, it is useful to get an intuition for the problem in a low dimensional space. We can get a picture of the total error, e(W), for our synthetic data as a function of the weight parameters w1 and w2:

eW[w1_,w2_] := Sum[{yy[[i]]-{w1,w2}.xx[[i]]}.{yy[[i]]-{w1,w2}.xx[[i]]},
		{i,1,Length[xx],1}]

g=Plot3D[eW[w1,w2],{w1,0,4},{w2,0,6},ViewPoint->{1.78,-2.861,0.312},AxesLabel->{"w1","w2","e"},ImageSize->Small];

If you plot the log of the error function, you can get a better view of where the minimum is:

g=Plot3D[Log[eW[w1,w2]],{w1,0,4},{w2,0,6},ViewPoint->{1.78,-2.861,0.312},AxesLabel->{"w1","w2","e"},ImageSize->Small]

(Remember that log[] is a montonic function, so it won’t change the location of the minimum.)

Note that because our input data are 2-element vectors, and our output or target values are 1-D, the matrix W we seek is not square--it has 1 row and 2 columns. So we represent W above as a vector. 
The minimum appears to be near {2,3}. It is easier to see the location of the minimum using ContourPlot[].

gc=ContourPlot[Log[eW[w1,w2]],{w1,1,3},{w2,2,4},Contours->16,ImageSize->Small]

Note: You can get a rough idea of the coordinates near the middle by moving your computer mouse over the above graphic after initializing the cell below:

Dynamic[MousePosition["Graphics"]]

Of course, in high dimensional spaces, W is huge, so we can't visualize it, and we can't hope to manually "adjust knobs" or sliders. But we can find the exact location of the minimum by finding "the bottom of the bowl".

Calculus tells us that this should be where the gradient of the error function, e is zero. Although there are a lot of indices to worry about, the familiar analogy from introductory calculus is to find the point at which the derivative of a function f(x) is zero, i.e. the value of x where f’(x)=0. 
e(W) plays the role of f(x), but W is a matrix of variables. The solution can be written concisely in terms of vector outerproducts, inversion, and matrix multiplication:

where

.

So we see that the value of W where the gradient is zero can be calculated by first adding up two sets of outer products, taking an inverse of one of them, and then multiplying. Let's translate the last line above into Mathematica instructions:

sumX = Sum[Outer[Times,xx[[i]],xx[[i]]],{i,Length[xx]}];

sumYX = Sum[Outer[Times,{yy[[i]]}, xx[[i]]],{i,Length[xx]}];

W = sumYX.Inverse[sumX]

The values for W come close to what we would expect from the structure of our data, a plane with parameters {2,3}.

2. Pseudoinverse to find W

For linear algebra afficianodos, there is another way of solving the linear least squares regression that uses the “pseudoinverse” of matrix to find the least-squares solution. For a square matrix X, the inverse of X is chosen so that XX-1 is equal to the identity matrix, I. The pseudoinverse, X*, of a rectangular matrix X is chosen so that XX* is close to the identity matrix in the sense that the sum of the squares of all of the entries of XX* - I is least. If interested, you can look up the derivation of X*. It can be written as:

or as:

For our simple generative model, let's take all of the input vectors x, and arrange them as columns in a matrix X:

Now do the same for the y's:

And what is the matrix that maps the x's to the y's with least squared error? Again, you can look up the derivation or try to derive it yourself. The answer is X*Y:

Inverse[Transpose[xx].xx].Transpose[xx].yy

X* is the pseudoinverse  (sometimes called the generalized inverse) of the matrix X. The PseudoInverse[] function is built into Mathematica, so now that we know about it, we can calculate the least-squares solution for our problem in one simple expression:

PseudoInverse[xx].yy

Question: What are the dimensions of the above PseudoInverse? Verify that PseudoInverse[X].X is the identity matrix with our data.

Chop[PseudoInverse[xx].xx]//MatrixForm

3. Gradient descent

Let’s go back to the original global error function that we used above for standard linear least squares regression. There we found the weights that gave the minimum total error by setting the gradient of the error function to zero, and then solving for the weights. 

Gradient descent, is a more general way of finding the minimum of an error function. It can be used when the error function is much more complicated (e.g. it isn’t quadratic), and there is no linear solution. But it works best for concave error functions, where one doesn’t have to worry about multiple minima. Later we will run into this last kind of situation  when we try to learn weights via supervised learning for a non-linear network with more than one layer of weights.

The idea is to start off at some location, W0 in weight space (W0 is just an initial guess), and iteratively move towards the minimum by always taking a step downhill. (See Graphical illustration of Gradient in Appendix).  Again we take the elements of W and think of them as defining a point in weight space, i.e. a vector in weight space. We want to find a vector, say Wbest, that tells us which direction to move in weight space so that that the error measured at location W0  decreases by as much as possible. We don't want to take too big of step or too small, so we need to control the length of Wbest, with a scalar η. In general, if we are at the ith point of an iteration, we would like:

(W^(i+1)) = (W^i)+ η (Wbest at (W^i))

To derive a rule for Wbest,  we assume that the error function is a smooth function of W. The downhill direction is given by the negative gradient of the error function. 
(This can be seen by noting that ∇e . (Overscript[n, ⋀]) = |∇e| cos(θ), which is biggest when θ = 0, in other words when (Overscript[n, ⋀]) points in the direction of ∇e.)

Again, the right hand side is just short-hand for

Sometimes the change in weights can bre modeled as a dynamical system:

We’ll see an example of this later.

and then approximate the differential equation by a discrete update rule. Wbest = - ∇e evaluated at the ith weight space point in the iteration. The figure below gives a graphical view in the small dimensional case of our example problem.

From the expression for the gradient which we wrote earlier in terms of outer products, we can obtain an expression for ΔW, with η taking the place of Δt:

Clear[wg];
eWgradient[wg_] := 
			Sum[Outer[Times,{yy[[i]]}, xx[[i]]], {i,Length[xx]}]-
			wg.Sum[Outer[Times,xx[[i]],xx[[i]]], 
				{i,Length[xx]}];

Now define a function T that specifies an updated weight matrix W=wg. You may have to adjust eta ( = η) to make sure the steps are sufficiently small to get convergence, but not so small as to take long to converge.

T[wg_] := wg + eta eWgradient[wg]

i=0; eta = .005; wglist = {};

wg={{RandomReal[5.0],RandomReal[5.0]}}; 
(*wg={{0,2}};*)

T[wg_] := wg + eta eWgradient[wg];


Clear[w1List]; w1List = NestList[T, wg, 40];
gl = ListPlot[Flatten[w1List, 1], Joined -> True, 
   PlotStyle -> {Orange, Dashed, Thick}];
   
Show[gc, gl]
w1List[[40]]

We used NestList[] above to keep track of all the iteration values so we could graph them, but we could just output the final value after, say 50 steps, using Nest[]

w1 = Nest[T,wg,50]

NOTE: There are often better algorithms for finding minima depending on the problem. For example, see the conjugate gradient method.

4. "Brain-style" learning: Iterative Widrow-Hoff learning to estimate W

So far so good. But there are a several problems. First, we are interested in brain-style computation. Based on what we think we know about neurons, how could the brain compute transposes, do matrix inversion and multiplication?

Second, when we learn we don't seem to gather information on a whole set of training pairs, and then suddenly build a memory matrix. A more plausible assumption is that learning is incremental, trial by trial, rather than batch.

Another problem is purely computational. What if the dimensionality of the vectors is really big? It is computationally expensive to invert large matrices. The above gradient descent procedure avoids the problem of inverting large matrices, but it involved computing a global error term over all the training pairs. We would like a method which would learn a regression mapping without having to store all the training pairs with the accompanying computation of a global error term. Instead, we'd like to compute an error term incrementally, trial-by-trial.

Can we discover the mapping W in such a way so as to be biologically plausible, and  avoid having to invert a large matrix? The Widrow-Hoff rule provides an answer. The basic idea behind Widrow-Hoff learning is to update W iteratively with each new training pair. Let's start off with an arbitrary W, find out which direction we would have to go in weight space to reduce the discrepancy between what W tells us x should map to and what it actually is, namely y. We take our clue from gradient descent, but apply it each time a new training pair comes along. We recompute the error term each time a new training pair comes along.

Let's try out this update rule on our synthetic training pairs.

ww1 = {{0,0}}; ww1list = {}; ww2list = {};

To visualize how the learning progresses, we are using ww1list and ww2list to store the first and second weights respectively, for each learning iteration.

eta = 5;

i=0;
While[i<Length[data],
	++i;
	in = {data[[i,1]],data[[i,2]]} ; out = {data[[i,3]]};
	ww1 = ww1 + (eta/i) Outer[Times,(out - ww1.in),in];
	ww1list = Append[ww1list,ww1[[1,1]]];
	ww2list = Append[ww2list,ww1[[1,2]]];
];
ww1

Keep executing the above input cell. You may have to run through the above loop several times before reaching stable convergence. We can plot up the two weights as a function of the iteration number to see how the Widrow-Hoff rule for weight modification eventually leads to two stable weights:

Dynamic[ListPlot[{ww1list,ww2list},PlotRange->{-1,8},AxesOrigin->{0,0},Joined->True]]

Note that we introduced a small change to the learning “constant” eta, replacing it by (eta/i )--it is no longer constant, but starts of large and decreases with each iteration.

Memory recall

We've seen several ways of finding the weights of a matrix that will approximately reproduce an output, given an input it has seen before. Let's try it out.

So in order to "recall" a response, from an input xx[[6]], we run it through the "network" memory matrix w1:

w1.xx[[22]]

And we can check to see how well it recalls:

yy[[22]]

Recall that the regression model of memory should generalize--interpolate and extrapolate. We saw a 3-dimensional example of interpolation in a graphical demonstration in the last lecture. But we can also extrapolate.

For example, {11,15} wasn't in the training set, but the expected output is:

w1.{11,15}

The network has "learned" a surface, (e.g. for one particular training set, y = 1.9x1 + 3.07x2  -- the coefficients of x1 and x2 will vary with each new random set of data)  through the points specified in the training set. 

A crucial point is that the linear associator will try to fit a plane (or hyperplane) through the data. If the data are not fitted well by that model, then the memory and generalization will not be good.

An obvious generalization is to fit hypersurfaces, rather than hyperplanes. And that is the direction we will head. But first let us look at linear regression from a point of view that you may not have thought of before.

Underconstrained problems and redundancy

Learning a one-to-many mapping:
	input dimensionality < output dimensionality

Motivation from vision

This example is very similar to the preceding example, except that we are going to try to learn 2D responses from 1D stimulus inputs. At first this doesn't seem to make sense, because the mapping from 1D to 2D in general would be expected to be underconstrained by the data. But this is exactly what is needed in certain kinds of inverse problems that are said to be "ill-posed" (Poggio et al., 1985). They are made "well-posed" when there is some underlying regularity (sometimes called "smoothness") in the high-dimensional output space (e.g. Kersten et al., 1987).

In the figure below, the data reaching the eye can be represented as an intensity value at each pixel. However, you perceive a shape, and can even adjust the orientation of a “gauge figure” consisting of circle tangent to the surface at a point, and a vector perpendicular to the tangent plane (like a “thumb tack”, cf. Nefs et al., 2006). This is an example of the shape-from-shading problem in vision. The gauge figure has two degrees of freedom at each point in the image.

Image[];

Thus "shape-from-shading" can be formulated as a problem of mapping N pixel intensities to N surface normals, where each surface normal is specified by two numbers--thus the mapping goes from N to 2N. E.g. given a 15x15 grid of intensities, estimate the 15x15 grid of 2D vectors that are normal to each point of the underlying surface. Thus one might start off with 225 numbers as input, but outputs 550 numbers.

GradientFieldPlot[f_,{x_,xmin_,xmax_},{y_,ymin_,ymax_},opts:OptionsPattern[]]:=VectorPlot[Evaluate[D[f,{{x,y}}]],{x,xmin,xmax},{y,ymin,ymax},opts]

ssurface[x_,y_]:= Sin[x y];
ggss=Plot3D[ssurface[x,y],{x,-4,4},{y,-4,4},MeshFunctions->{#1&,#3&}];
gsss=Plot3D[ssurface[x,y],{x,-4,4},{y,-4,4},ViewPoint->{0,0,30},Axes->False,Lighting->{{"Directional",White,{{1,0,1},{1,1,0}}}},Mesh->15,MeshStyle->Directive[Thin,Blue]];
gf=GradientFieldPlot[ssurface[x,y], {x, -4, 4}, {y, -4, 4},Frame->False];
Show[GraphicsRow[{ggss,gsss,gf}]]

The left panel illustrates an underlying surface that generates the image intensities shown in the middle. The shape-from-shading problem tries to take the graylevel values at points (approximated by grid boxes) in the middle panel, and output a 2-D vector for that grid, as shown in the right-hand panel.

Let's try a simple version of learning a mapping from a low to high-dimensional case. We'll use low-dimensional synthetic data whose generation process we'll keep hidden for now.

Synthetic data -- try not to look closely...

Here is the generative model:

output = {x2, y} = {a x1, b x1 + noise}

We want to learn the parameters {a,b} = {w1,w2}, from training pairs of inputs and outputs: 
{x1, {{x2, y}}}

r3Dline[a_,b_]:=N[Table[{x1=1 RandomReal[],a x1,b x1+0.5 RandomReal[]-0.25},{30}],2];
data=r3Dline[2,3];
xx=data[[All,1]];
yy=data[[All,2;;3]];

So the input data is a list of 1D scalars, and the output data a list of 2D vectors. Let's apply the Widrow-Hoff algorithm to learn the relationship between the input stimuli, xx, and the output responses, yy.

Iterative Widrow-Hoff learning

w1 = {{0},{0}}; w1list = {}; w2list = {}; eta = 4.0;
Dimensions[w1];

i=0; 
While[i<Length[data],
	++i;
	out = {data[[i,2]],data[[i,3]]} ; in = {data[[i,1]]};
	w1 = w1 + (eta/i) Outer[Times,(out - w1.in),in];
	w1list = Append[w1list,w1[[1]][[1]]];
	w2list = Append[w2list,w1[[2]][[1]]];
];
w1

ListPlot[{w1list,w2list},PlotRange->{0,6},AxesOrigin->{0,0},Joined->True]

Visualizing the underlying structure

So why does it work to learn to predict a 2D value from a 1D input? The reason, of course, is that the underlying structure of the output data is highly constrained, and in fact lies close to a straight line in 3-space.

Show[Graphics3D[Point/@data],Axes->True,AxesLabel->{"x1","x2","y"},Ticks->None,ImageSize->Medium]

Here is the generative model for our synthetic data:

x1 -> {x2, y}, 
and given the form is a straightline through the origin with some additive noise, the output {x2, y} is given by:
	{x2, y} = {a x1, b x1 + noise}

An input, a scalar quantity x1 specifies a plane. We've found a line that intersects the plane and for which points in our data set our close to that line. The intersection point gives us the coordinates {x2, y} for our regression fit, or "recall".

We learned the parameters {a,b}, corresponding to the weights,  from training pairs of inputs and outputs: {x1, {{x2, y}}}

r3Dline[a_, b_] := N[Table[{x1 = 1 RandomReal[], a x1, b x1 + 0.5` RandomReal[] - 0.25`}, {30}], 2];

data = r3Dline[2,3];

It produced the coordinates of a noisy line in 3-space.

Making estimates of high-dimensional functions from lower dimensional inputs happens a lot, and is common to many so called "early vision" problems (Poggio et al., 1985). For more recent work on learning low-level vision, see Freeman et al., 2000.

Introduction to multi-layer nets

This was an earlier exercise. For linear networks, no computational power is gained by having extra layers:

y1 := W0.y0;
		 y2 := W1.y1;

is equivalent to:

y2 := W1.(W0.y0) := W1.W0.y0 := W3.y0;

where W3 is just another matrix. However, if the inner product is followed by a non-linear transformation, then concatenating layers of neural elements allows one to compute more complex transformations:

y2 :=  f[W1.f[(W0.y0)]];

where f[ ], for example, is a sigmoidal “squashing function” or logistic function:

f[x_]:=1/(1+Exp[-(x-.5)]);

NOTE: Notation alert--we are using f or f2 to represent a scalar squashing function, not an input vector.

Exercise: 2 input units, 3 hidden units, 1 output unit

Clear[y2,f2];
f2[x_]:=N[1/(1+Exp[-(x-.5)*100])];
W1={{-13.2328, 6.06398, 6.04958}};
W0={{-0.937564, -1.09841}, {-9.95589, 3.85642}, {3.79034, -10.1737}};
y2[y0_] :=  Chop[f2[W1.f2[(W0.y0)]]];

What logical function does this network compute?

{y2[{0,0}], y2[{0,1}],y2[{1,0}],y2[{1,1}]}

What logical function will the above network compute with the following weights?

W0 = {{1.9009195689645, 1.9251997561975753}, {2.177312432892121, 2.198392078532626}, 
  {-2.8297064197335953, -2.8648165030618693}};
  W1 = {{1.1241018372719358, 1.7600034188676417, -10.048562036295957}};

Clear[y2];
y2[y0_] :=  Chop[f2[W1.f2[(W0.y0)]]];

{y2[{0,0}], y2[{0,1}],y2[{1,0}],y2[{1,1}]}

The demonstration illustrates that we can compute XOR as well as  OR using a network with more than one layer of weights.

Multi-layer nets: Learning and the error back-propagation algorithm

So we've seen that the non-linear generic neuron allows increased computational power when we add an extra layer of weights.

Suppose we have a multi-layer network with 3 layers of weights. The output from the first layer is: 
y1 = f[u1] = f[W1.y0]. The output of the second layer is: y2 = f[u2] = f[W2.y1]. And so forth. To simplify a bit, we’ll assume each layer has the same number of units and weights:

The problem is: how to assign the weights? We’d like to find them through supervised learning.

For any complex system that is required to achieve a target goal, for the system to work, each component must contribute towards the goal. If the goal is not met, one has to figure out which component needs to be fixed. If the goal is met, each component contributed something towards the goal. How does one assign the credit for success or failure to a component? 

This problem is called the credit-assignment problem. In particular, for the above multi-layer network, how do we adjust the weights, which intrinsics to the whole network, in a way appropriate for learning a given input/output relation?

For a given input {(y^0)= (x^p)}, we feed forward the information to the last layer (layer L) to produce an output {y=(y^L)}. We compare the output to the target value supplied by the "teacher" {t = (t^p)}, and compute the error as the sum of squared differences:

where the sum is over all N output units. (For simplicity, we left out the superscript p in  (Subscript[t, k]^p). The subscript k in (Subscript[t, k]^p) means the component of the corresponding vector of activity ((t )^p).) λ indexes the weight layers going from λ=1 to L. Note that the y's at any point after the input depend on the u's (the weighted sum before the non-linearity), each of which in turn depends on all the (Subscript[w, ij]^λ)s before it.

The trick is to find out how to assign credit (and blame) for the error to each of the weights. Gradient descent provides the answer. Adjust the weights such that:

new (Subscript[w, ij]^λ) = previous (Subscript[w, ij]^λ) + (Subscript[Δw, ij]^λ)

where

Again, this formula means that if we calculate the gradient ∇E=(∂_(Subscript[w, ij]^λ)) E, its negative direction points in the direction of steepest descent. Thus if we update the weight vector to point in this direction, then at the next step we should in general have a lower value of E. I.e. 

	E(new (Subscript[w, ij]^λ)) < E(previous (Subscript[w, ij]^λ)).

The appendix shows how to calculate a weight adjust with just one layer of weights. This case is related to the area of Generalized Linear Regression (not to be confused with the General Linear Model!). 

We won’t go through the math here for the multiple layer case--it is really just a very complicated application of taking derivatives where the complication arises in keeping track of all the indices. Instead we assume we have the rule that tells us what the effective error term--the delta value-- and adjustment that needs to be for each weight, and see how it gets put together to converge on weights that tend to minimize the overall output errors.

Summary of backprop algorithm

1. Initialize the weights to small random values
2. Pick a pattern from the input/output collection, say the pth pattern: {xp,tp}. Run the input vector, xp , feedforward through the network. This will generate a set of values (Subscript[u, i]^L) and (Subscript[y, i]^L) in all the nodes of the network. Keep in mind that (Subscript[u, i]^L)  is the linear weighted sum of its inputs arriving from layer L-1:
(Subscript[u, i]^L) = (∑_(j=1)^3)(Subscript[w, ij]^L)(Subscript[y, j]^(L-1)). And  (Subscript[y, i]^L) =f((Subscript[u, i]^L)).

Calculate a delta term (analogous to the Widrow-Hoff rule) for the output layer L:

Note why it is important to have an expression for the derivative of the squashing function f(). The derivative has a particularly nice form when f(u) = 1/(1+(e^(-u))). (see Appendix, Logistic function - a smooth, differentiable non-linearity)

3. Propagate the errors back through the layers:

...the error back propagation or "back prop" part.

4. Calculate weight adjustments (analogous to the outerproduct part of the Widrow-Hoff) and update using:

5. Repeat steps 2 to 4 until convergence.

One can accumulate the weight adjustments for each training pair, and then update them all at once.

In practice, updating the weights after each training pair often works better than accumulating a bunch of input/output pairs, and then computing the cumulative global error. The reason is that by randomly sampling a training pair, the "descent" may actually climb the global error function defined by the entire set. As we will see later with the Boltzmann machine, occasional climbing is useful to avoid local minima.

There are a number of derivations and illustrations on the web that you might find useful. See for example, http://galaxy.agh.edu.pl/~vlsi/AI/backp_t_en/backprop.html

Also, Andrew Ng has several excellent coursera videos on error backpropagation.

Backprop simulation example: XOR

With appropriate weights, 2 weight layers with 3 hidden units can solve the XOR problem. But this is still a tough problem to learn, mainly because it requires that two very different inputs map to the same output. See the supplementary material for a Mathematica demo that learns the weights for solving the XOR problem.

Next

Relationship of neural networks to machine learning

Problems of slow learning and over-fitting in MLPs. 

Non-linear networks & “energy”: Hopfield

Probabilistic models. From energy to probability.  Boltzmann machines

Appendix

Pseudoinverse solution for the second (underconstrained redundancy) case above

r3Dline[a_,b_]:=N[Table[{x1=1 RandomReal[],a x1,b x1+0.5 RandomReal[]-0.25},{30}],2];
data=r3Dline[2,3];
xx=data[[All,1]];
yy=data[[All,2;;3]];

We can calculate the memory matrix using the PseudoInverse in this case too:

Dimensions[xx];
Dimensions[yy];

matmem = Transpose[PseudoInverse[Transpose[{xx}]].yy]

Recall

Let's check out a few values to see how well the memory matrix can recall a 2 dimensional output, given a one dimensional input.

{{"recall: ",matmem.Transpose[{xx}][[5]], matmem.Transpose[{xx}][[13]], matmem.Transpose[{xx}][[22]]},
{"true values: ", yy[[5]], yy[[13]], yy[[22]]}}//Grid

Another function for calculating the least-squares solution

Regression and learning in linear neural networks. Last time we showed 4 different ways to find the generating parameters {2,3} for the following data:

rsurface[a_,b_]:=N[Table[{x1=1 RandomReal[],x2=1 RandomReal[],a x1+b x2+0.5 RandomReal[]-0.25},{120}],2];
data=rsurface[2,3];
yy=data[[All,3]];
xx=data[[All,1;;2]];

Linear regression is so common that Mathematica has added the following function to find the least squares parameters directly:

LeastSquares[xx,yy]

Graphical illustration of Gradient

The function f(x,y) is plotted as a contour plot. The demo below calculates and plots the vector gv = ∇f. After normalized to unit length, it is plotted as red arrow. The red arrow should always point down the "hill" in the contour plot. You can manipulate the position of the point at which the gradient gets evaluated with the a and b sliders.
The θ slider controls the direction of a second vector un (in black) whose length is determined by the projection:  ∇f.un

Clear[a,b, gv,un,gm,g2,g3,x,y];
f[x_,y_]:=2 x^2 + 5 y^2 - Sin[x];
g1=ContourPlot[f[x,y],{x,-3,3},{y,-3,3},PlotPoints->50,ImageSize->Medium];
Manipulate[
gv=N[D[2 x^2 + 6 y^2 - Sin[x],{{x,y}}] /.{x->a,y->b}];
gn=-Normalize[gv];
un={Cos[θ],Sin[θ]};
gm=(un.gn) un;
g2=Graphics[{Red,Arrow[{{a,b},gn+{a,b}}]}];
g3=Graphics[Arrow[{{a,b},gm+{a,b}}]];
Show[{g1,g2,g3}],{{a,1},-3,3},{{b,2},-3,3},{θ,0,2*Pi}]

Logistic function - a smooth, differentiable non-linearity

As we saw above, it can be useful to have a non-linearity which is smooth enough to be differentiable. D[] returns an expression for the derivative, so to define a function that is the derivative of another we write:

Df1[x_] := D[f[t],t] /.t->x
Df1[x_] := Evaluate[D[f[t],t]]

The more direct way is to use "operators" or functionals that take functions as inputs and return functions as outputs. Derivative[], or f'[x]  return functions:

(*Df[x_] := Derivative[f[x],x]*)
Df[x_] := f'[x]

Note that the derivative has a particularly simple expression in terms of f[x], which you can verify by comparing Simplify[Df[x]] with Simplify[hh[x]], or Simplify[Df[x]-hh[x]]:

f[x_]:=1/(1+Exp[-x]);

hh[x_] := f[x](1-f[x]);

Here is a plot of the sigmoid and its derivative:

Plot[{Df[x],f[x]},{x,-3,3}]

Exercise: calculate the gradient for the last layer in a back-prop network

Let's see how to calculate the gradient for the last layer. Then we will generalize the result to an expression that tells us how to adjust the ijth weight for the connection in the λth layer.

Remove["Global`*"]

We start with the error term for the jth training pair. Note that (Subscript[y, j])depends on the inputs x, but these are fixed. At this point, we care about the variable weights.

(1/2)*((Subscript[t, i])-(Subscript[y, i])[(Subscript[w, ij])])^2

Take the derivative with respect to (Subscript[w, ij])

eq1=D[(1/2)*((Subscript[t, i])-(Subscript[y, i])[(Subscript[w, ij])])^2,(Subscript[w, ij])]

y depends on u which in turn depends on the w's, so using the chain rule in calculus,  we can write  (Subscript[y, i])'((Subscript[w, ij]))as:

eq2=D[(Subscript[y, i])[(Subscript[u, i])[(Subscript[w, ij])]],(Subscript[w, ij])]

But (Subscript[y, i])'((Subscript[u, i])) is the derivative of the squashing function f: (Subscript[y, i])'((Subscript[u, i]))=f'((Subscript[u, i])) So we can substitute f' for y in eq2

eq3=eq2 /. (Subscript[y, i])'[(Subscript[u, i])[(Subscript[w, ij])]]->f'[(Subscript[u, i])]

Now u is given by:

(Subscript[u, i])=(∑_(k=1)^N)(Subscript[x, k])(Subscript[w, i,k])

Let's pick an arbitrary input input, say i=5. Suppose we take the derivative with respect to (Subscript[w, j,i]), with i = 5? Then with respect to (Subscript[w, j,5]), we have:

D[(∑_(k=1)^10)(Subscript[x, k])(Subscript[w, i,k]),(Subscript[w, i,5])]

So from this we can guess that in general (i.e. any i) we should have:

(Subscript[Δw, ij])=((Subscript[t, i])-(Subscript[y, i])( x; (Subscript[w, ij]))) f' ((Subscript[u, i])) (Subscript[x, j])

where (Subscript[x, j]) is the jth "input" to the last layer. This is the delta-rule. But this delta rule only works for the output layer. We need to know how to update all the weights.

We don't have direct access to the hidden units, and the problem is how to find the error signals for the hidden layers. It turns out that the delta error terms can be propagated back in terms of a weighted sum of the delta terms at the level above. We will see how this update rule gets applied to the weights starting from the top--i..e. those closest to the outputs, and then applied recursively down towards the inputs. Hence "back-propagation". We will relabel (Subscript[x, k])to be (Subscript[y, k]^λ), where λ indicates the layer.  

So (Subscript[y, k]^0)corresponds to the set of inputs to the network (bottom of figure).
To derive the recursive rules to relate deltas at an earlier layer to the next higher layer takes some work and careful book keeping with indices. We won't go through the details here. Let's look at a summary of the algorithm that results.

References

Bishop, C. M. (1995). Neural Networks for Pattern Recognition. Oxford: Oxford Univeristy Press.
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.
William T. Freeman, Egon C. Pasztor and Owen T. Carmichael (2000). Learning Low-Level Vision. International Journal of Computer Vision 40, 25-47
Kersten, D., O'Toole, A. J., Sereno, M. E., Knill, D. C., & Anderson, J. A. (1987). Associative learning of scene parameters from  images. Appl. Opt., 26, 4999-5006.
Knill, D. C., & Kersten, D. (1990). Learning a near-optimal estimator for surface shape from shading. CVGIP, 50(1), 75-100.
Nefs, H.T., Koenderink, J.J. & Kappers, A.M.L. (2006). Shape-from-shading for matte and glossy objects. Acta Psychologica, 121 (3), 297-316 
Poggio, T., Torre, V., & Koch, C. (1985). Computational vision  and regularization theory. Nature, 317, 314-319.
Geman, S., Bienenstock, E., & Doursat, R. (1992). Neural networks and the bias/variance dilemma. Neural Computation, 4(1), 1-58.