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

Self-Organizing Networks

Self-Organizing (Kohonen) Maps with linear units

References for this Notebook:
Kohonen, T., Self-Organization and Associative Memory, Springer Verlag, Berlin, 1984
Fort, J.C., Solving a combinatorial Problem via Self-Organizing Process: An Application of the Kohonen Algorithm to the Traveling Salesman Problem, Biol. Cybern., 59, p.33-40, 1988.

Informal Problem Statement:

We will explore the dynamics of self-organization and map formation in linear units.

Network Architecture

The network consists of an array of ordered units or cells, each with m inputs.  The array is viewed as a  map. One- and two-dimensional maps are particularly useful but higher dimensional arrays can be considered as well. Each cell independently receives the same m-dimensional input vectors.

Learning Rule: Adjust only  the Winning Cell and its immediate neighbors

In the simplest version of the algorithm proposed by Kohonen, upon receiving each data sample,  the position in the array of the maximally responding cell , i.e., the cell whose weight vector most closely resembles the data vector, will be identified.  The learning process requires that winning cell and its immediate neighbors   adjust the weights in such a way as to reduce the difference between their weight vector and the data vector.  The amount of reduction applied to neighbors  decreases with distance from the winning cell.  

After repeated iterations,  the locations of the peak responses will be
 	1) correctly ordered,
  2) distributed (mapped) in a smooth fashion, 
  (the distribution will approximate the probability distribution of the data (to the extent that it was revealed by the samples drawn during the iterations)
 3)  contracted around the mean of the distribution.  
 The amount of contraction varies inversely with the variability of the input data.
 After global convergence, the process will continue to show small fluctuations since each data sample will disturb the distribution slightly. The process may be terminated whenever a "satisfactory" level of spatial organization has been reached.  What constitutes a satisfactory solution  depends of course on a particular application. In practice, best results are obtained when the amount of correction applied is decreased progressively during the iterations while at the same time the extent of the neighborhood is narrowed.
This learning rule involves a global operation,  namely, the localization of the maximally responding cell.  This may be inappropriate in some applications and can be obviated by providing lateral connections between neighboring cells with positive (excitatory) weights (for the near neighbors) that become negative (inhibitory) as distance increases.

Retrieval (Execution) Rule:

In the retrieval or execution phase, data is applied to all units in parallel.  In the absence of lateral connections, the units respond individually and independently of one another. 
"Retrieval" consists of observing the position of the cell responding maximally on the map constituted by the cell array.   This position topologically identifies the input vector parameters.

Input Data - Uniform Distribution

In this notebook the cells have two input ports, but the corresponding data vector is constrained to have unit length. The angular position theta of this vector around the unit circle will be the only independent parameter. Thus the input domain is one-dimensional and can be visualized as points on the circle.
For instance:

angle=N[Pi (2 Random[]- 1)]
datum={Cos[angle],Sin[angle]}

Network

We now establish a one-dimensional array of cells by  creating an (numData x 2)  array of  weightVectors for the two inputs. 
The initial values for the weights will  be chosen as random angles around the circle and will form vectors of unit length,  directly comparable to the data vectors.

Clear[u]
u:=N[Pi(2 Random[]- 1)]

numData =20
weightAngles=Table[u, {numData}];
Short[
weightVectors=Transpose[{Cos[weightAngles],Sin[weightAngles]}]]

According to the learning rule,  for each presentation of a datum vector, the adjustment can be pictured as a rotation of the winning cell and also, but to a lesser extent, of its neighbors,  toward  the datum vector.
The amount of "pull" must decrease with distance.
We will choose a pullProfile function of distance  that describes the relative amount of weight modification imposed on the neighbors . 
We will choose an even function with value one at distance zero and monotonically decreasing with the absolute value of distance. 
A simple exponentially decaying function will do.  The parameter spread will let us adjust the extent of interaction.  We plot the function with spread=0.5.

Clear[y, distance,spread,pullProfile]
pullProfile[y_, spread_] := Exp[ -spread Abs[y] ]

Plot[pullProfile[distance, .5],
{distance,-8,8},PlotRange->{0,1},AxesOrigin->{0,0}]

To get a discrete set of influence coefficients to apply to the cell array, we make a table of the pullProfile function, with an odd number of entries, centered on the zero distance (the winning cell). The parameter  pullExtent  sets the number of neighbors on each side of the winning cell that will receive an adjustment

pullExtent = 4;
spread=0.5;

To apply the adjustment in parallel to the whole array of cells,  we need to insert these influence coefficients into an array of length numData.  Most of the entries will be zero since only the immediate neighbors are affected.

mask = Join[ Table[
					pullProfile[distance,spread],{distance,-pullExtent,
					pullExtent}],
 					Table[0,{numData-2*pullExtent-1}] ]

Finally, to facilitate subsequent indexing, we  rotate this array to place the maximum value (unity)  in the first position.  During the calculations,  it will be shifted by the amount necessary to bring it into alignment with the "winning" cell,  i.e., the cell displaying the maximum netInput.

pullFactors =RotateLeft[mask , pullExtent]

Let us check the length of the new array

Length[pullFactors]==numData

Network Learning

We will first run one iteration of the network and analyse each step.  
First we generate one data item,  datum,  and look at the "net inputs" (the weighted sums) it produces in the cells of the array.

theta=u;
	datum= N[{Cos[theta],Sin[theta]}];
	Short[netInputs = weightVectors.datum]

Next, we find the position index of the  winning cell  .

Position[ netInputs, Max[netInputs]]

We just need the number inside the list (Flatten will do) and we should allow for the unlikely event of more than one maximum, in which case we will just take the first to show up in the array.

peakIndex=First[Flatten[%]]

We can now align the array of pullFactors with the peakIndexth cell by shifting  to the right by an amount equal to peakIndex -1.

pullFactors= RotateRight[pullFactors,peakIndex-1];

We are now in position to adjust the weights in parallel in the entire weightVectors array. Specifically, we will rotate the weightAngles toward the datum angle theta,  in proportions set by the pullFactors array and multiplied  by a gain parameter, gain,  that controls the size of the successive adjustments. At this time  gain is set to a constant value.

gain = 0.5;

Let us apply the learning rule

winningCellWeightsBefore=weightVectors[[peakIndex]];
newWeightAngles = 
	weightAngles + gain (theta-weightAngles) pullFactors;	
newWeightVectors=
Transpose[{Cos[newWeightAngles],Sin[newWeightAngles]}];

winningCellWeightsAfter=newWeightVectors[[peakIndex]];

Show[ 
	Graphics[ {
			RGBColor[0,0,1], PointSize[0.06],
			Point[datum], Line[{{0,0},datum}],
			RGBColor[0,0,1],Line[{winningCellWeightsBefore,
			winningCellWeightsAfter}],
			 PointSize[0.02],
			 Point[winningCellWeightsBefore],
			 RGBColor[0,0,1],PointSize[0.04],
			 Point[winningCellWeightsAfter]}
	 ],
Axes-> True,
AxesOrigin->{0,0},AxesLabel->{{w1},{w2}}];

We can see that the vector representing the maximally responding cell has moved closer to the datum.  We can also look at the current state of the entire array by showing the irregular polygon formed by the weight vectors of successive cells in the array

Show[ 
	Graphics[ {
			PointSize[0.02] , RGBColor[1,0,0], 
			Map[Point, weightVectors],
			Line[weightVectors],		
			RGBColor[0,0,1], PointSize[0.04],Thickness[0.01],
			Point[datum], Line[{{0,0},datum}],
			RGBColor[0,1,1],Point[winningCellWeightsAfter]
			}],
Axes->True,
AxesOrigin->{0,0},AspectRatio->1]

Spatial Self-organization

We will now reset the problem parameters and iterate the process to see how the cell weights ( represented as points around the circle)  self-organize when input data from a uniform distribution is applied repeatedly and weights are corrected according to the rules just described.
The output graphs can be animated to better illustrate the transformation.

Clear[weightAngles,newWeightAngles,weightVectors,newWeightVectors]

numData =50;
weightAngles=Table[u, {numData}];
weightVectors=Transpose[{Cos[weightAngles],Sin[weightAngles]}];

We set new parameters and recalculate pullFactors exactly as detailed earlier.

pullExtent = 5;
spread=0.6;

pullFactors =RotateLeft[Join[
	Table[ pullProfile[distance,spread],
	       {distance,-pullExtent, pullExtent}],
	Table[0,{numData-2*pullExtent-1}] ],pullExtent];

We will now iterate the network.  To avoid an excessive number of  graphics outputs,  we will look at the configuration periodically, after  iterInterval steps,  and continue the process for numSnap successive such iteration intervals (each graph will be separated by iterInterval iterations).

We need to introduce two additional heuristics.  At the beginning the network is very disorganized and large corrections are desirable.  These should be progressively reduced as spatial organization settles in.  This can be accomplished by using a decay factor for the adjustment gain gain.  To further improve convergence,  it is also desirable to start with a broad spatial spread of adjustment  (determined by the value of pullExtent).  This will require recalculating the correction table periodically and will be omitted at this time.

numSnap=10;
iterInterval=50;
gain=0.3;
decay=1/(numSnap*iterInterval);
iterNumber=0;
datum=winningCell={0.,0.};

Do[
Show[ 
	Graphics[ {
			PointSize[0.02] , RGBColor[0,0,1], 
			Map[Point, weightVectors],
			Line[weightVectors],		
			RGBColor[0,0,1], PointSize[0.06],Thickness[0.01],
			Point[datum], Line[{{0,0},datum}],
			RGBColor[0,1,1],PointSize[0.06], Point[winningCell]},
Axes->False,			
AxesOrigin->{0,0},AspectRatio->1,PlotRange->{{-1,1},{-1,1}},
PlotLabel -> StringForm["Iteration Number ``",iterNumber]]];
Do[
	theta=u;
	datum= {Cos[theta],Sin[theta]};
	netInputs = weightVectors.datum;
	peakIndex=
	Flatten[Position[ netInputs, Max[netInputs] ]][[1]];
		
	factors = gain RotateRight[pullFactors,peakIndex-1];
	newWeightAngles = 
weightAngles + (theta-weightAngles) factors;	
	newWeightVectors=
	Transpose[{Cos[newWeightAngles],Sin[newWeightAngles]}];
weightAngles = newWeightAngles; weightVectors=newWeightVectors;
winningCell=weightVectors[[peakIndex]];
gain=gain - decay gain;iterNumber++,
{iterInterval}], 

{numSnap}];