Matrix Multiplication

Nengo Tutorial: Matrix Multiplication

In this tutorial, you will learn how to use nengo to multiply two n x n matrices (for some fixed n) in a new network that can be easily added to any of your existing networks. The code to create the network will be in a python script.

The following techniques are covered by this tutorial:

Creating a new network
Creating nested networks
Adding nested networks to existing projects
Using python to automate ensemble creation
Connecting terminations and origins between networks
Loading data in Matlab

In addition, exercises are provided at the end.

Step 1: Create a New Network

When you run nengo the workspace is blank. You must create a new network. Select the command File|New|Network. A dialog box will appear prompting you for a name and some other options. Enter MatrixMultiplication as the name and click Ok.

Step 2: Create a Nested Network

To keep your project organised and modular it is advisable to group neural populations into nested networks. In this tutorial you will add a subnetwork to the MatrixMultiplication network: click on the icon. A yellow box should appear around the icon indicating that it is selected. Now, select File|New|Network again. The same dialog box should appear as in step 1. Enter Multiplier as the name and click Ok. If your new subnetwork is not centered in the MatrixMultiplication network viewer, you can scroll the Network Viewer window by holding down the mouse button inside of it and dragging the cursor in the direction you want to scroll. You can maximize the MatrixMultiplication network viewer by pressing the box with blunt corners in the upper right hand side of the MatrixMultiplication network viewer window.

Step 3: Create a Python Script to Automate Ensemble Creation

Now you will add ensembles to implement matrix multiplication. This will be done automatically using a python script. An ensemble will be created for each cross term multiplication and then their outputs will be summed. The creation of the cross terms will be done automatically with a python script. In order to add a python script to your project, select the command View|Open Script Editor. This should open a new window titled Script Editor. Create a new script by selecting the command File|New. Then save it with the command File|Save As. Use the filename init_multiplication.py in the Save As text box.

Copy the following script into the Script Editor. The lines starting with the symbol # are comments explaining what the script is doing. This code should also be provided at the location from which you downloaded this tutorial.

# init_multiplication.py:

# By default, the python engine does not import any of the nengo # classes, so this must be done before any of the nengo classes # can be instantiated:
from ca.nengo.model.nef.impl import *
from ca.nengo.model.neuron import *
from ca.nengo.math.impl import *

# This function will create the NEFEnsembles for matrix
# multiplication.
# network - This is the network to which we wish to add matrix
# multiplication.
# populationSize - This is the number of neurons in each
# NEFEnsemble representing the entries and cross terms.
# tau - This is the synaptic time constant of the terminations.
# n - the dimension of the matrices

def init_multiplication(network,n,populationSize,tau):

# The NEFEnsembleFactory interface lets you create

# NEFEnsembles with given properties. Here, we have asked for
# default properties.

ef = NEFEnsembleFactoryImpl();

# First, we will create NEFEnsembles to represent the
# entries of each matrix.

a = [ [None]*n for i in range(0,n) ];
b = [ [None]*n for i in range(0,n) ];
aT = [ [None]*n for i in range(0,n) ];
bT = [ [None]*n for i in range(0,n) ];
aO = [ [None]*n for i in range(0,n) ];
bO = [ [None]*n for i in range(0,n) ];

r = [ [None]*n for i in range(0,n) ];
rT = [ [None]*n for i in range(0,n) ];
rO = [ [None]*n for i in range(0,n) ];

IdentityFunction = Polynomial([0,1]);
Product = PostfixFunction(None,"x0*x1",2);

for i in range(0,n):
for j in range(0,n):

# We use ef to create an ensemble to represent the value
# of the i][j-th entry of the first matrix.
a[i][j] = ef.make("a"+str(i+1)+str(j+1), \
populationSize, 1, "a"+str(i+1)+str(j+1), 0);

# A termination is added. Other ensembles in the parent
# network will project to this terminal and set the
# i,j-th value of the first matrix.
aT[i][j] = a[i][j].addDecodedTermination \
("a"+str(i+1)+str(j+1), [[1]], tau, 0);

# An origin is added from this ensemble.
# It will project to the cross terms.
aO[i][j] = a[i][j].addDecodedOrigin \
("a"+str(i+1)+str(j+1), [IdentityFunction], \
Neuron.AXON);

# The ensemble is added to the network.
network.addNode(a[i][j]);

The terminal is exposed to the parent network.
network.exposeTermination \
(aT[i][j] , "a"+str(i+1)+str(j+1));

# This is repeated for the second matrix.
b[i][j] = ef.make("b"+str(i+1)+str(j+1), \
populationSize, 1, "b"+str(i+1)+str(j+1), 0);
bT[i][j] = b[i][j].addDecodedTermination \
("b"+str(i+1)+str(j+1), [[1]], tau, 0);
bO[i][j] = b[i][j].addDecodedOrigin \
("b"+str(istr(j+1), [IdentityFunction], \
Neuron.AXON);
network.addNode(b[i][j]);
network.exposeTermination \
(bT[i][j], "b"+str(i+1)+str(j+1));

# Now, for each row i and column j we must sum the cross
# terms to produce the i][j-th entry of the result.

for i in range(0,n):for j in range(0,n):

# This ensemble is the i][j-th entry of the result.
r[i][j] = ef.make("r"+str(i+1)+str(j+1), \
populationSize, 1, "r"+str(i+1)+str(j+1), 0);
network.addNode(r[i][j]);

for k in range (0,n):
# This is the k-th cross term between the i-th row
# of the matrix 'a' and the j-th column of the matrix
# 'b': a[i,k]b[k,j].

# Each coordinate of this two dimensional representation
# is in the range [-1,1], so the magnitude of the
# largest vector we might have to represent is
# sqrt(1+1)=sqrt(2)=1.414

c = ef.make \
("a"+str(i+1)+str(k+1)++str(k+1)+str(j+1),\
populationSize, [1.414,1.414]);

# The ensemble is added to the network.
network.addNode(c);

# The i,k-th entry of matrix 'a' will project to the
# first dimension of this cross term.
ca = c.addDecodedTermination \
"a"+str(i+1)+str(k+1)+"b"+str(k+1)+str(j+1)+"_a", \
[[1.0],[0.0]], tau, 0);
network.addProjection(aO[i][k], ca);

# The k,j-th entry of matrix 'b' will project to the
# second dimension of this cross term.
cb = c.addDecodedTermination \
"a"+str(i+1)+str(k+1)+"b"str(k+1)+str(j+1"_b", \
[[0.0],[1.0]], tau, 0);
network.addProjection(bO[k][j], cb);

# This origin will decode the product of the two
# dimensions of the cross term population.
cO = c.addDecodedOrigin \
("a"+str(i+1)+str(k+1)+"b"+str(k+1)+str(j+1), \
[Product], Neuron.AXON);

# The i,j-th entry of the result is the sum of
# these cross terms, so the cross terms must project .
rT = r[i][j].addDecodedTermination( \
"a"+str(i+1)+str(k+1)+"b"+str(k+1)+str(j+1),\,br>[[1.0]],tau,0)0;
network.addProjection(cO,rT);

# Add an origin to the i][j-th entry and expose it
# so that it can be used by parent networks.

rO = r[i][j].addDecodedOrigin("r"+str(i+1)+str(j+1), \
[IdentityFunction], Neuron.AXON);
network.exposeOrigin(rO, "r"+str(i+1)+str(j+1),);

Step 4: Run the python Script and Create Ensembles Automatically

The python script should now be in the current directory. To run it from nengo, first you must open the script console. Use the command View|Open Script Console. A text box should appear at the bottom of the screen into which you can enter python code. To run the script, type the following in the script console:

run init_multiplication.py

The script defines a new python function :

init_multiplication(network, n, PopulationSize, tau)

The first argument, network, is the network to which ensembles will be added. The second argument, n, is the dimension of the matrices you want to multiply. PopulationSize is the number of neurons that are used to represent the matrix entries and each cross term. tau is the synaptic time constant for all the terminations added in the model. To start the automatic ensemble creation, first ensure that the Multiplier network is selected. Do this by clicking on it and confirming that the icon is outlined with a yellow box. In the script console, type the following:

init_multiplication(that, 2, 100, 0.1)

In this function call, the symbol that refers to the currently selected node. If you selected the Multiplier network, then the ensembles should be added to that network. If you didn't select anything then python will likely return an error.

Double click on the Multiplier network. This will open the Multiplier network viewer. When you create ensembles automatically, they are displayed in a pile by default. To make the network look a bit nicer, ctrl+click somewhere on the black background of the Multiplier network viewer. Use the command Layout|Algorithm to sort the ensembles in various ways according to their topology, or use the command Layout|Sort|Name to sort them alphabetically.

Step 5: Specify the Input into the Multiplier Network

In order to specify the input functions to the Multiplier network, you must first Show the origins and terminations that were exposed in the automatic ensemble creation process. Close the Multiplier network viewer and ctrl+click on the Multiplier network. Use the command Origins and Terminations|Show all. The green circles are the terminations and the light blue circles are the origins. If they have no names or if the names are really small, you can zoom in to see them better. If your computer has a touchpad, you may do the touchpad scrolling technique where you roll two fingers in the direction you want to scroll (up is zooming in, and down is zooming out.)

The input to the terminations can be constant functions, time varying functions, or decoded information from other neural populations. Constant functions are the easiest place to start. You will create 3 constant functions for the values -1,0,1 and then connect them to the terminals of the Multiplier network. ctrl+click on the black background of the MatrixMultiplication network viewer. Select the command Create New|Function Input. A dialog box should appear:

Enter One into the Name box. Enter 1 into the Output Dimensions box. Click Set Functions. Another dialog box should appear asking for the details for the function you are creating. By default, the function will be a Constant Function. Verify that Function 0 is set to Constant Function and then click Set.

Another dialog box should appear asking for the value of this constant function. Enter 1 in the Value box and click Ok.

Do the same thing for Constant Functions named Minus One and Zero and make sure to set the Value box appropriately.

If you want to multiply the matrices:

Then you should drag the origins from the Function Inputs to the corresponding terminations of the Multiplier network.

Step 6: Simulate Matrix Multiplication

The origins of the Multiplier network currently project nowhere. But you can add probes to decode the function representation directly from the result ensembles. ctrl+click on the Multiplier network. Under the title Data Collection, select the command Add Probe|r11 Function of NEFEnsemble State and repeat this command for all other entries of the result from which you would like to record data.

Now you are ready to simulate 2x2 matrix multiplication of constant matrices. ctrl+click on the black background of the MatrixMultiplier network viewer and select the command Simulator|Run MatrixMultiplier Network. Enter 0.0 for the Start Time, and 0.001 for the Time Step, and 1.0 for the End Time and then click Ok.

After the simulation completes, you may plot the probed data. ctrl+click on the purple triangle (the probe) next to the origin corresponding to an entry of the resulting matrix. Use the command Plot|Plot w/ options. This will let you convolve the weighted spike train with a filter and plot a function represented by the result population. Refer to the textbook Neural Engineering for more information.

A dialog box should open. Enter 0.01 in the box for Time constant of display filter. Enter 0 in the Subsampling box and then click Ok.

A window should appear displaying a graph of the decoded matrix entry. The value might be a bit different from the expected solution. This is caused by the distortion induced by the optimal choice of weights and also by saturation of neural populations. Again, refer to the book Neural Engineering for more information.

Appendix 1: Exporting Data to matlab

To manipulate the neural data in other programs, ctrl+click on the probe from which you want to save data. Select the command Export data|Matlab. You will then be prompted for a filename. The file will be saved in the current directory unless you specify a full path. In this example, export data from the probe r11 and use the filename exported_data.

Now the data must be loaded in matlab. First, start a new matlab session. The exported data should be in a file in the directory in which you installed nengo called exported_data.mat. matlab may have a different current directory, so you might have to provide a full path when you use commands to describe and load the .mat files. A .mat file stores matlab variables. To discover what variables are stored in the file exported_data.mat, and to load them use the commands:

>>whos -file ~/simulation/nengo-1.0/exported_data.mat

      Name             Size              Bytes  Class     Attributes

      r11           1000x1                8000  double              

      r11_time         1x1000             8000  double

    >> load ~/simulation/nengo-1.0/exported_data.mat

    >>

Now, the variable r11 has the raw weighted spike train data. In order to decode a representation of the matrix entry, we must convolve with a time filter. Filter the data with an exponential decay filter using the following commands:

>> timeFilter = exp(-0.01*[1:100]);

    >> # normalise the filter:

    >> timeFilter = timeFilter/sum(timeFilter);

    >> r11_filtered = filter(timeFilter,1,r11);

Now the variable r_11 contains data that can be manipulated and plotted:

>> plot(r11_time, r11_filtered);

    >> MSE = 1/numel(r11_filtered)*(sum(r11_filtered-1))

    -0.2254

Exercises:

1. Add non-constant inputs and experiment with multiplying 2x2 time varying matrices.
2. Create two matrix multiplication networks and chain them together to compute the product of three matrices: ABC.
3. Create Matrix multiplication networks of other dimensions and explore the relationship between the number of neurons in the cross term population and the size of the matrices and the error in the representation.
4. Modify the python script so that populations are created with larger radii. Explore the relationship between the size of the radii and amount of distortion and saturation.
5. Modify the python script so that matrices of different dimensions can be multiplied.
6. Use recurrent network connections and your solution from exercise 5 to implement a simulation of the linear homogeneous differential equation dx/dt = A(t)x, where A(t) is a time varying nxn matrix, and x is an n dimensional vector.