Basins of Attraction & Visualizing Learning Systems

NKS 2003 Summer School Project / NKS2004 Student Presentation

Philipp Michel (philipp.michel@yale.edu)
Yale University
April (25^th),  2004

Basins of Attraction & Visualizing Learning Systems

1. Introduction & Motivation

2. Fundamentals: Basins of Attraction in Newton's Method

3. Learning using the Hopfield Model

4. Visualizing Basins of Attraction in the Hopfield Model

5. Cyclic Sequential Substitution Systems

6. Searching for Simple Classification Behavior

7. Conclusion & Future Work

Introduction and Motivation

Learning systems categorize previously unseen data as belonging to a particular class of patterns derived from familiar data

Such classification essentially identifies those training examples an unseen item is most "similar" to

Determining which basin of attraction in the space of system variables an arbitrary point lies in can be considered analogous to performing classification on the basis of similarity

NKS simple programs could potentially implement classification by exhibiting the same behavior for similar initial conditions

Project Aims

1. Consider a classic example of basins of attraction (Newton's method) by visualization & look for interesting, complex behavior

2. Investigate & visualize attractors in a basic learning system: the Hopfield model of neural networks

3. Try to find NKS-style simple programs that might exhibit elementary classification behavior: Sequential Substitution Systems

Fundamentals: Basins of Attraction in Newton's Method

A basin of attraction is the set of points in the space of system variables such that initial conditions chosen in this set dynamically evolve to a particular attractor

Newton's method uses the iterative map (Subscript[x, n + 1])=(Subscript[x, n])-((f((Subscript[x, n])))/(f'((Subscript[x, n]))))and an initial condition (Subscript[x, 0])to find the roots of f(x)

Visualization of attractor basins for polynomials using Newton's root finding algorithm

Iterate Newton's method for various polynomials taking starting approximations from the complex plane, plotting the argument of the result after convergence or 30 iterations. Use CountourPlot or Raster.

Show[
      Graphics[Raster[
          Table[(((Arg[FixedPoint[# - f[#]/(f')[#] &, i + I*j, 30]] + π)/(2 π))), {i, -5, 5, 
              0.025}, {j, -5, 5, 0.025}], ColorFunction -> Hue], 
        AspectRatio -> 1]]

Example: f[z_] := (z^7) - 1

Plotting the modulus of the result for g[z_]:=(x^7) + 3  (x^5) + 4  (x^4) + 3 (x^2) + 1

Zooming in to the center:

Example: h[z_] := (z^3) - 1

This scheme works for other iterative root finding algorithms. Here's Householder's method for g[z] above:

And some more: w[z_] := ((z^4) - 1) ((z^4) - 4)

u[z_] := ((z^4) - 1) ((z^4) - 4)

For polynomials of degree 2 or higher, the resulting map is a fractal if there are 3 or more distinct roots

Learning Using the Hopfield Model

The Hopfield model implements an associative memory

"Store a set of p patterns ((Subscript[ξ, i])^μ) in such a way that when presented with a new pattern (Subscript[ζ, i]), the network responds by producing whichever one of the stored patterns most closely resembles (Subscript[ζ, i])"

Formulation & implementation

Units are updated synchronously or asynchronously using the rule (Subscript[S, i]) :=  sgn((∑_j) (Subscript[w, ij]) (Subscript[S, j])), where the weight matrix (Subscript[w, ij])is computed from the training patterns using the Hebbian learning rule (Subscript[w, ij]) = (1/N)\(∑\+\(μ = 1\)\%p\( ξ\_i\^μ\) ξ\_j\^μ\)  (N = # of units).

tosign[x_ /; x > 0]  :=  1;
tosign[0]  :=  -1;
fromsign[x_ /; x > 0] := 1;

    
w[patterns_List]  :=  
      Module[{a  =  patterns}, 
a  = 
          Map[Flatten,  a]; 
a = 
          Map[Map[tosign, #] &, a]; 
a = 
          Map[Outer[Times, #, #] &, a]; 
(Plus @@ a)/
          N[Length[a]]
]

    
hop[n_, W_, input_]  :=  
      Module[{s = Map[tosign, Flatten[input]]}, 
l  =  
          Length[W]; 
origl  =  
          Length[input]; 
Do[
i  =  
            Random[Integer, {1, l}]; 
s[[i]] = 
            Sign[W[[i]] . 
                s], 
{n}]; 
Partition[s, 
          Length[s]/origl]
]

    
fixedhop[W_, input_]  :=  
      Module[{s = Map[tosign, Flatten[input]], 
          sprev = {}}, 
origl  =  
          Length[input]; 
While[
          sprev ≠ s, 
sprev = 
            s; 
Do[
s[[i]] = 
              Sign[W[[i]] . s], 
{i, 
              Length[s]}]
]; 
Partition[
          s, Length[s]/origl]
]

To perform classification:
1. Apply Hebb's rule to training patterns, get weight matrix.
2. Store test pattern's k bits at the k units of the network.
3. Iterate update rule until fixed point is reached. The bit pattern now stored in the units is the (or an approximation to the) most similar stored pattern.

Experimenting with simple patterns

Some simple bitmaps of numbers

\(\(four\  = \ {{0, 0, 0, 0, 0, 1, 1, 0, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 0, 1, 1, 
            1, 0, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 
            1, 1, 1, 1, 0, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, \ 1, 1, 0, 1, 
            1, 0, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 1, 1, 
            0, 0, 1, 1, 0, 0, 0}, 
	\ \ {1, 1, 0, 0, 0, 1, 
            1, 0, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {1, 1, 1, 
            1, 1, 1, 1, 1, 1, 
            1}, 
\ \ \ \ \ \ \ \ \ \ \ {1, 1, 1, 1, 1, 1, 
            1, 1, 1, 1}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, \ 0, 
            0, 0, 1, 1, 0, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 0, 0, 1, 
            1, 0, 0, 0}\ };\)\)
\(\(seven\  = \ {{0, 0, 0, 0, 0, 0, 0, 0, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 1, 1, 1, 1, 1, 
            1, 1, 1, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 1, 1, 
            1, 1, 1, 1, 1, 1, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, \ 0, 0, 0, 0, 
            1, 1, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 
            0, 0, 0, 1, 1, 0, 0}, 
	\ \ {0, 0, 0, 0, 0, 0, 
            1, 1, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 
            0, 0, 0, 1, 1, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 0, 0, 0, 
            1, 1, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, \ 0, 
            0, 0, 0, 1, 1, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 0, 0, 0, 
            1, 1, 0, 0}\ };\)\)
\(\(eight = {{0, 0, 0, 0, 1, 1, 0, 0, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 1, 1, 1, 
            1, 0, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 1, 
            1, 1, 1, 1, 1, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 1, \ 1, 0, 0, 0, 
            1, 1, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 1, 1, 
            0, 0, 0, 0, 1, 1, 0}, 
	\ \ {0, 0, 1, 1, 1, 1, 
            1, 1, 0, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 1, 1, 
            0, 0, 0, 0, 1, 1, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 1, 1, 1, 0, 0, 
            1, 1, 1, 0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, \ 1, 
            1, 1, 1, 1, 1, 0, 
            0}, 
\ \ \ \ \ \ \ \ \ \ \ {0, 0, 0, 0, 1, 1, 
            0, 0, 0, 0}\ };\)\)
ListDensityPlot[Reverse[four], ColorFunction -> Automatic,  
        Mesh -> False];

Train the Hopfield network on those patterns

W = w[{four, seven, eight}];

Create a noisy test pattern

noisy  = 
        four  +  Table[If[Random[] < 0.1, 1, 0], {10}, {10}];

ListDensityPlot[Reverse[noisy /. {1 -> 2, 2 -> 1}], 
      ColorFunction -> Automatic,  Mesh -> False]

Iterate the Network and see how the output changes as the number of iterations are increased

ListDensityPlot[
      Reverse[hop[5, W, noisy] /. {1 -> 2, 2 -> 1}], 
      ColorFunction -> Automatic,  Mesh -> False]

ListDensityPlot[
      Reverse[hop[10, W, noisy] /. {1 -> 2, 2 -> 1}], 
      ColorFunction -> Automatic,  Mesh -> False]

ListDensityPlot[
      Reverse[hop[20, W, noisy] /. {1 -> 2, 2 -> 1}], 
      ColorFunction -> Automatic,  Mesh -> False]

ListDensityPlot[
      Reverse[hop[100, W, noisy] /. {1 -> 2, 2 -> 1}], 
      ColorFunction -> Automatic,  Mesh -> False]

We get perfect reconstruction. Note however that as soon as you have >2 input patterns there are situations where the network will settle in a state that does not correspond to an input pattern!

Iterating until convergence

Some bigger bitmaps as input patterns: beer and monkey!

A noisy monkey

Thresholding the bitmaps to be black & white only and reducing size

Iterating until convergence gives the output

Visualizing Basins of Attraction in the Hopfield Model

Two dimensional plots

First, train a Hopfield network on input patterns consisting of the padded binary representation of some numbers in the range 0 - (2^16) - 1

Then construct a 2D plot by using the concatenated & padded binary representations of the x and y axes as the test pattern. After convergence, the number corresponding to the binary output pattern is plotted. Each color thus corresponds to a different attractor.

result = 
        Table[
FromDigits[
fixedhopn[
              W, 
Join[PadLeft[IntegerDigits[i, 2], 8], 
                PadLeft[IntegerDigits[j, 2], 
                  8]]
], 
2], 
{i, 0, 255}, {j, 0, 255}];

Show[
      Graphics[Raster[(res/65534), ColorFunction -> Hue], 
        AspectRatio -> 1]]

For one training pattern, this looks like the following. Note how we get two colors due to the possibility of the network settling in the reverse state (the binary complement of the most similar input pattern)

Two training patterns

Three training patterns

Ten training patterns

Pretty pictures that don't really tell us much. The structure of the basins of attraction is not really revealed here. The concatenation mapping from 2D to 1D adopted for visualization purposes is somewhat arbitrary and not as natural as the 2D decomposition of the complex numbers in the plots for Newton's method.

One dimensional plots

Simple listplots of the test pattern against the fixed point solution of the Hopfield network

We want to eliminate reverse state solutions by just replacing them by their corresponding input pattern

elim[in_, patt_]  :=  
        Module[{rev  =  
              patt  /. {0 -> 1, 
                  1 -> 0}}, 
in  /. {rev -> 
                patt}
];
elimreversed[input_, tpatterns_]  :=  
        Fold[elim, input, tpatterns];

For increasing values of k, choose k random input patterns and iterate over the numbers 0 - (2^16) - 1as the test patterns

The diagonal becomes more prevalent as the number of input patterns is increased and a symmetry about the diagonal becomes evident

For 23 patterns. Nesting?

State transition graphs

[cf. NKS pp. 962]

State transition graphs depict how many test patterns converge to a particular attractor. The following STG shows all (2^8) - 1test patterns converging to one of the two input patterns

For four input patterns

Not too helpful. For 23 input patterns

This does show how the size of the basins of attraction varies considerably. Unfortunately, it only tells us how many patterns converge to a particular attractor, but not which ones.

Cyclic Sequential Substitution Systems

Goal: Investigate attractors in sequential substitution systems, look for rules that might exhibit classification behavior when given particular initial conditions

Controlling unlimited growth in sequential substitution systems using cyclic boundary conditions

Many sequential substitution systems grow indefinitely, e.g.

To study attractors, we want systems that halt (don't evolve further) or exhibit cyclic behavior in the majority of cases. To accomplish this, proceed in a similar manner to cellular automata: limit the growth of the system at least somewhat by imposing cyclic boundary conditions

cycapply[str_, rule_] := 
        Module[{lhslen  =  Length[First[rule]],  parts, posl, pos, 
            rotated,  result}, 
parts  =  
            Partition[str, lhslen, 1, 1]; 
posl  =  
            Position[parts, First[rule], 1, 1]; 
If[
            Length[posl] == 0, 
            Return[str]]; 
pos  =  
            posl[[1, 1]]; 
rotated  =  
            RotateLeft[str, pos - 1]; 
result  =  
            RotateRight[
              Flatten[ReplacePart[
                  Partition[rotated, lhslen, lhslen, {1, 1}, {}], Last[rule], 
                  1]], pos - 1]];

cycapplylist[str_, rules_] := 
      Module[{s = str, sorig = str, riter = rules}, 
While[
          s == sorig, 
s = 
            cycapply[s, First[riter]]; 
If[
            Length[riter] > 1, riter = Drop[riter, 1], 
            Return[s]];
]; 
Return[
          s];
]

SSSEvolveListCyc[rules_, init_, t_Integer] := 
        NestList[cycapplylist[#, rules] &, init, t];

Plotting cyclic substitution systems

This uses NiceRasterGraphics started from DeBruijnSequence[{0,1},10]

CycSSSPlot[rules_, init_, t_Integer] := 
      Module[{sys = SSSEvolveListCyc[rules, init, t], 
          longest = 0}, 
longest =  
          Max[Map[Length[#] &, sys]]; 
NiceRasterGraphics[
          Map[PadRight[#, Max[Length[init], longest]] &, sys], 2, 
          1]
]

Here are some interesting cases

CycSSSPlot[{{0, 1, 0} -> {1, 1, 0}, {1, 0} ->  {0, 0}}, dbs, 
        1000] // Show

CycSSSPlot[{{0, 0} -> {0, 1, 1}, {1, 1, 0} -> {0, 1}}, dbs, 
        2200] // Show

CycSSSPlot[{{0, 1, 1} -> {1}, {1, 1, 1, 0, 1, 0, 1, 0} -> {0, 
              1}, {1, 1} -> {0}, {0, 0} -> {1, 1, 1}, {0, 1, 
              0} -> {1, 0}}, dbs, 600] // Show

CycSSSPlot[{{1, 0, 1} -> {0, 1, 0}, {0, 0, 0} -> {1, 1}, {1, 
              0} -> {1, 1}}, dbs, 5000] // Show

Using John Bacus's drawing method for the first rule above, the structure being slowly undone from the right becomes apparent as a slowly moving circular pattern:

Very dark areas indicate that the evolution has halted. Consider:

Searching for Simple Classification Behavior

Enumerating possible cyclic sequential substitution systems

To perform search, we need to be able to enumerate all possible rules for cyclic sequential substitution systems and find a unique mapping from the rules to ℕ. We need to be able to code pairs of numbers and finite lists of numbers. Use <| x, y  |  > =(2^x) (2  y + 1) and  < x, y>=(2^x) (2  y + 1) - 1.

angbar[x_, y_] :=  (2^x) (2 y  + 1);

ang[x_, y_] := (2^x) (2 y  +  1) - 1;

codelist[in_List] := 
        Fold[angbar[#2, #1] &, 0, Reverse[in]];

decodelist[in_Integer] :=  
      Module[{t = {}, s = {}, result = {}, h = 0, 
          l = 0}, 
t = 
          Reverse[IntegerDigits[in, 2]]; 
While[
          t ≠ {}, 
h = 
            First[t]; 
If[
            h == 1, 
s = 
              Append[s, 0]; 
t = 
              Drop[t, 1], 
l = 
              Length[First[Split[t]]]; 
s = 
              Append[s, l]; 
t = 
              Drop[t, l + 
                  1]
]
]; 
s
]

decangbar[in_Integer] := 
      Module[{t = {}, x = 0, y = 0}, 
t = 
          Split[Reverse[IntegerDigits[in, 2]]]; 
If[
          MemberQ[First[t], 1], 
x = 
            0; 
y = 
            FromDigits[Reverse[Drop[Flatten[t], 1]], 
              2], 
x = 
            Length[First[t]]; 
y = 
            FromDigits[Reverse[Drop[Flatten[Drop[t, 1]], 1]], 
              2]]; 
{x, y}]

decang[in_Integer] := 
      Module[{t = {}, x = 0, y = 0}, 
t = 
          Split[Reverse[IntegerDigits[in, 2]]]; 
If[
          MemberQ[First[t], 0], 
x = 
            0; 
y = 
            FromDigits[Reverse[Drop[Flatten[t], 1]], 
              2], 
x = 
            Length[First[t]]; 
y = 
            FromDigits[Reverse[Drop[Flatten[Drop[t, 1]], 1]], 
              2]]; 
{x, y}]

coderules[rules_List] := 
        codelist[
          Map[angbar[codelist[First[#]], codelist[Last[#]]] &, 
            rules]];

decoderules[in_Integer] := 
        Map[First[#] -> Last[#] &, 
          Map[decodelist[#] &, 
            Map[decangbar[#] &, decodelist[in]], {2}]];

This is decent when trying to map rules to integers. But in the converse case, we find that the mapping is sparse and many integers don't map to any sensible rule. Also, the numbers get ridiculously huge

coderules[{{0, 0} -> {1, 1}}]

coderules[{{0, 0} -> {1, 1}, {1, 1, 1} -> {0}}]

Instead, have a function that just enumerates all rules directly. This one is due to Matthew Szudzik

br[{}, rules_] := rules
\(\(br[list_, rules_] := 

        With[{len = Length[list]}, 

          Apply[Sequence, 

            Table[br[Take[list, {j + 1, len}], 

                Append[rules, 
\(Join[{a___}, 
                          Take[list, i], {___}] :> {Length[{a}] + i - 
                            1, #} &\)[
Take[list, {i + 1, j}]]]], {i, 1, 

                Min[First\  /@ \ ReplaceList[list, rules], {len}]}, {j, i, 
                len}], {0, 1}]
]\)\(
\)
      \)
\(AllSSS[k_, 
          t_] := 
\(br[#, {}] &\) /@ 

            Table[IntegerDigits[i, k, t], {i, 0, 
                k^t - 1}]\  /. \ 
\(({_, a__, _} :> {_, b_})\) -> \(({a} -> 
                b)\)\)

Searching for rules with interesting behavior

I started searching through all possible rules of total length (number of cells in the entire rule icon) 1,2,3,4,..., running each from an initial condition of DeBruijnSequence[{0,1},10] for a couple of hundred steps.

Length 1:

GraphicsArray[
        Partition[Map[CycSSSPlot[#, dbs, 500] &, AllSSS[2, 1]], 2]] // 
      Show

Length 2:

GraphicsArray[
        Partition[Map[CycSSSPlot[#, dbs, 500] &, AllSSS[2, 2]], 5]] // 
      Show

Length 3

GraphicsArray[
        Partition[Map[CycSSSPlot[#, dbs, 500] &, AllSSS[2, 3]], 6]] // 
      Show

The search space blows up fairly quickly after this. Consider

Length[AllSSS[2, 4]]

Length[AllSSS[2, 5]]

Length[AllSSS[2, 6]]

Length[AllSSS[2, 7]]

Rules that blow up or totally die out are less likely to be interesting than ones which settle to a particular configuration (although this might depend exactly on which initial condition is used). A complete search would have to consider not only all rules, but also all possible initial conditions to check for classification-like behavior where the final state or cycle reached could be considered an attractor (e.g. a "majority vote" rule which behaves differently when 1's dominate the initial condition than when 0's do).

A more targeted search for specific desired classification behavior might yield better results than just looking for rules that exhibit "interesting" behavior on a De Bruijn sequence as the initial condition.

Conclusion & Future Work

Traditional learning & classification systems are constructed explicitly to be able to exhibit complex behavior. Can simple rules do just as well? What is the simplest possible classifier system?

What about Turing Machine or CA based systems?

What is the best way of visualizing / analyzing the relationship between initial conditions and the corresponding attractors (system behavior)?

Next steps: How do you effectively structure large searches for classification behavior, particularly in an infinite rule space as in the case of Cyclic SSS? Perhaps iterate through possible initial conditions for only the very simplest "interesting" rules?