from  
  Exploring Neural Networks with Mathematica™ 
  © Jacques J. Vidal,  1990

Multi-layered Feedforward Networks - Backpropagation

We will look at the backpropagation algorithm in the context of feedforward, layered networks (Figure 1), which have found broad applications in the literature. These networks have one layer of input ports, the points of application of the input vectors.  From each input port, connections radiate toward the next level,  the first "hidden" layer of  intermediate units. Units are continuous valued,  typically with a superSigmoid transfer function. There can be  one or several hidden layers.
Finally, the output layer contains the units that are considered outputs.   To apply the backpropagation algorithm,  which is a supervised learning rule, a training set of inputs must be given along with the corresponding  "desired outputs". The algorithm is a generalization of the Widrow-Hoff procedure. Observed errors serve to determine the direction of successive corrections to the weights.  However, since weights affect hidden units,  error measures must be propagated backward to find the correcting terms for the hidden units.

Informal Problem Statement:

Specifically, we will explore the simplest two-layered networks of continuous units with a sigmoid output function. 
We will test the Backpropagation Algorithm on the EXOR logical function. The EXOR and the EQUIVALENCE functions partition  input data  into two classes that are not linearly separable.

Processing Unit - Formulating the output function

We will use the superSigmoid function, a scaled-up version of the classical sigmoid function :

Clear[superSigmoid]

superSigmoid[x_] =N[-m((1 - E^(a x))/(1 + E^(a x)))];

dsuperSigmoid[y_]=superSigmoid'[y]

With these values of a and m  the  function will cross the values -1 and  +1 at  the the two points where the second derivative is extremal.

Network

A 2-layer three node (2-1) feedforward network will be used. (Figure 2) . Note that there is also (not shown below) a bias input into the top three units.

Training Set

We first create  the training data, namely,  the input and output lists for the logical function selected (the EXOR). Recall that the first component of the input vector is always constant and provides us with a bias.

xVectors =  List[{-1,-1,-1},{-1,-1,1},{-1,1,-1},{-1,1,1}];

To complete the training data,  outputs must be associated
with each xVector. For the EXOR function:

desiredOutputs = List[ -1, 1, 1, -1 ];

We choose the three sets of initial weights randomly, in appropriate list structures

u:=2*Random[]-1;

topWeights =2.0 * {u,u,u};
hiddenWeights =2.0 * {{u,u,u},{u,u,u}};

n=2;
numData=2^n (* 2^n is the number of distinct data entries *)

v:=Random[Integer,{1,numData}];

Forward Pass

We will generate one input-output pair by selecting one of them at random in the function truth table.  Then we will calculate (forward) and print out the net inputs and outputs with the initially selected weight values.

index = v;
 Off[General::spell1];
	Print[ 
		xVector = xVectors[[index]],"  ",
		desiredOutput = desiredOutputs[[index]],"  ",
		
		bottomNetInputs = hiddenWeights.xVector,"  ",
  bottomOutputs = 
     Prepend[superSigmoid[bottomNetInputs], -1],"  ",
		
		topNetInput = topWeights.bottomOutputs," ",
		topOutput = superSigmoid[topNetInput]
		];
		On[General::spell1];

At this point we can calculate the cost function for the initial set of weights:

startingCost = Apply[Plus,  
0.5 (
	desiredOutputs -
		superSigmoid[topWeights.Transpose[Map[Prepend[#,-1]&,
		       Transpose[superSigmoid[hiddenWeights.Transpose[xVectors]]]]]]
    )^2]

This cost function calculation subsumes the feedforward "execution" phase of the network.

Backward Error Propagation

Now we will update the weights by backpropagating the errors,  starting with the top layer (here a single node):

cost :=  Apply[Plus,  
0.5 (
	desiredOutputs - superSigmoid[topWeights.Transpose[
		Map[Prepend[#,-1]&,Transpose[superSigmoid[hiddenWeights.
		                           Transpose[xVectors]]]]]])^2];

initTopWeights :={u,u,u};
initHiddenWeights :={{u,u,u},{u,u,u}};

eta = 0.1;

We will repeat the procedure, and produce numPlots plots of the global cost, as weights are being updated, to illustrate the dependency of convergence on iniital conditions.  All runs will converge but at different rate.

numPlots=16;

Do[
topWeights=initTopWeights;
hiddenWeights=initHiddenWeights;
iterMaximumCount = 300;
iterCount = 1;
topWeightSequence=
hiddenWeightSequence=costSequence={};
While[(Chop[cost,0.001]>=0 && iterCount<iterMaximumCount),
Do[
index = v;
	xVector = xVectors[[index]];
		desiredOutput = desiredOutputs[[index]];
		bottomNetInputs = hiddenWeights.xVector;
		bottomOutputs = Prepend[ superSigmoid[bottomNetInputs], -1];
		topNetInput = topWeights.bottomOutputs;
		topOutput = superSigmoid[topNetInput];
		 deltaErrorTop =
		   ((topOutput - desiredOutput)*
				 dsuperSigmoid[topNetInput]);

updateTopWeights = eta* deltaErrorTop * bottomOutputs;
	topWeights = topWeights - updateTopWeights;
updateHiddenWeightsLeft =
			  (eta*dsuperSigmoid[bottomNetInputs[[1]]])*
			  deltaErrorTop*topWeights[[2]]*xVector;		  
		
			updateHiddenWeightsRight =
			  (eta*dsuperSigmoid[bottomNetInputs[[2]]])*
			  deltaErrorTop*topWeights[[3]]*xVector;		
	
	hiddenWeights[[1]] = hiddenWeights[[1]] -
	  updateHiddenWeightsLeft;
	hiddenWeights[[2]] = hiddenWeights[[2]] -
	  updateHiddenWeightsRight;
				 (*Now calculate cost with the new weights:*)
costSequence=Append[costSequence,cost];
topWeightSequence=Append[topWeightSequence,topWeights];
hiddenWeightSequence=Append[hiddenWeightSequence,
hiddenWeights]; Increment[iterCount]]];
costPlot = ListPlot[costSequence,PlotRange->{0,5}],
{numPlots}];

The next plot illustrates the trajectories of the three weight vectors in one typical run.

Show [ Graphics3D [{Thickness[0.002], RGBColor[1,0,0],
Line[topWeightSequence],
RGBColor[0,1,0],
Line[Transpose[hiddenWeightSequence][[1]]],RGBColor[0,0,1],
Line[Transpose[hiddenWeightSequence][[2]]]}],
 PlotRange ->Automatic, Axes->Automatic]

Empirical evidence indicates that the cost function defined above has some undesirable but unavoidable properties.  It is of course non convex and its global minimum may not be unique. It can have local minima, although experience indicates that they do not constitute a  major problem: they are indeed few since they require that all partial derivatives vanish simultaneously and these minima can be by-passed by various techniques.  More serious is the fact that ravines and plateaus form easily.
In ravines the derivatives are large in one direction and small in another.  On plateaus,  cost is high but gradients are low so that each step brings only minute improvements. Plateaus are a consequence of the saturation characteristics of the superSigmoid function. Significant learning  occurs only when the units are in the central,  near linear, region. With large weights,  large NetInputs will occur,  bringing the units into the saturation zone.  One possible solution is to force the weights to systematically decrease with the iterations, but this technique has the potential drawback of bringing the weights to zero,  i.e., to a guaranteed plateau existing around the origin, before reaching convergence. 
The consequence of these problems is that,  although the procedure almost always converges,  the number of iterations for convergence varies widely with the initial choice of weights.  There is a zone (probably fractal in configuration) from which convergence is rapid.  But other areas will require orders of magnitude more iterations to reach a solution.
It must also be realized that there will be an infinity of possible solutions (weights range) generally forming nonconnected regions in the weight space.