C implementation

A neural network has two modes of operation: computation, and training. In computing mode, the user fixes the activity of certain input neurons, and the other neurons x_i will calculate:

   x_i  <--  1 / [ 1 + exp(-∑_j w_ij x_j - b_i) ]

Here x_i represents the activity level of neuron i, w_ij is the strength of the connection from neuron j to neuron i, and b_i is the neuron’s sensitivity. In training mode, we will update the weights and biases of each neuron using the following formula:

   w_ij  <--  w_ij + η · x_i x_j

    b_i  <--  b_i + η · x_i

The parameter η determines the learning rate, which should be large enough to train quickly, but not so big that the algorithm overshoots and destabilizes the network. This training algorithm (J. R. Movellan, Contrastive Hebbian learning in interactive networks, 1990) works for symmetric networks (w_ij = w_ji) such as we will use. As the saying goes, neurons that fire together, wire together.

These two basic operations of a neural network are encoded by simple formulas, but since those formulas will be used many times in the course of using a network they are the time-consuming step. We will therefore write that part of our program is C. Notice that our outermost function has the same form as main(), but we will use a different name so as not to conflict with Cicada’s own main() function.

NN.c

#include "NN.h"
#include <math.h>


// runNetwork(): evolves a neural network to a steady state
// Takes the params: 1 - weights; 2 - neuron activities; 3 - input; 4 - step size
// (& additionally, in training mode): 5 - target output; 6 - learning rate

ccInt runNetwork(ccInt argc, char **argv)
{
   neural_network myNN;
   double *inputs, step_size, *target_outputs, learning_rate;
   int i, numInputs, numOutputs;
   
   
   /* ----- set up data types, etc. ----- */
   
   
   for (i = 0; i < numInputs; i++)
       myNN.activity[i] = inputs[i];
   for (i = numInputs; i < myNN.numNeurons; i++)
       myNN.activity[i] = 0;
   
   if ( argc == 6 )   {     // i.e. if we're in training mode
       
       if (getSteadyState(myNN, numInputs, step_size) != 0) return 1;
       trainNetwork(myNN, -learning_rate);
       
       for (i = 0; i < numOutputs; i++)
           myNN.activity[numInputs + i] = target_outputs[i];
       
       if (getSteadyState(myNN, numInputs+numOutputs, step_size) != 0) return 1;
       trainNetwork(myNN, learning_rate);        }
   
   else if (getSteadyState(myNN, numInputs, step_size) != 0) return 1;
   
   
   /* ----- save results ----- */
   
   
   return 0;
}


// getSteadyState() evolves a network to the self-consistent state x_i = f( W_ij x_j ).

int getSteadyState(neural_network NN, int numClamped, double StepSize)
{
   const double max_mean_sq_diff = 0.001;
   const long maxIterations = 1000;
   
   double diff, sq_diff, input, newOutput;
   int iteration, i, j;
   
   if (numClamped == NN.numNeurons) return 0;
   
   
       // keep updating the network until it reaches a steady state
   
   for (iteration = 1; iteration <= maxIterations; iteration++)   {
       sq_diff = 0;
       
       for (i = numClamped; i < NN.numNeurons; i++) {
           input = 0;
           for (j = 0; j < NN.numNeurons; j++)  {
           if (i != j)  {
               input += NN.activity[j] * NN.weights[i*NN.numNeurons + j];
           }}
           newOutput = 1./(1 + exp(-input));
           
           diff = newOutput - NN.activity[i];
           sq_diff += diff*diff;
           NN.activity[i] *= 1-StepSize;
           NN.activity[i] += StepSize * newOutput;
       }
       
       if (sq_diff < max_mean_sq_diff * (NN.numNeurons - numClamped))
           return 0;
   }
   
   return 1;
}


// trainNetwork() updates the weights and biases using the Hebbian rule.

void trainNetwork(neural_network NN, double learningRate)
{
   int i, j;
   
   for (i = 0; i < NN.numNeurons; i++)   {
   for (j = 0; j < NN.numNeurons; j++)   {
   if (i != j)   {
       NN.weights[i*NN.numNeurons + j] += learningRate * NN.activity[i] * NN.activity[j];
   }}}
}

NN.h

typedef struct {
   int numNeurons;    // 'N'
   double *weights;   // N x N array of incoming synapses
   double *activity;  // length-N vector
} neural_network;

extern int runNetwork(int, char **);
extern int getSteadyState(neural_network, int, double);
extern void trainNetwork(neural_network, double);

Last update: May 8, 2024