多層パーセプトロンのC言語による実装を示します。
この実装は、下記のWebページにあった実装をもとにしています。
http://www.ccodechamp.com/c-program-of-multilayer-perceptron-net-using-backpropagation/
主な変更点は以下のとおりです。
- ReLUを追加。
- 最適化アルゴリズムをAdaGradに変更（An overview of gradient descent optimization algorithms）。
- 乱数をSFMTに変更（SIMD-oriented Fast Mersenne Twister (SFMT)）。
- いくつかの固定長の配列を動的配列に変更。
- loadNet関数を削除。
# include <stdio.h>
# include <math.h>
# include "SFMT-src-1.4.1/SFMT.h"
# define DTYPE float
# define M_PI 3.14159265358979323846
//================================================================
// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/
sfmt_t *sfmt;
unsigned int initRand(unsigned int nSeed)
{
sfmt = (sfmt_t *) malloc(sizeof(sfmt_t));
if (sfmt == NULL) return 1;
sfmt_init_gen_rand(sfmt, nSeed);
return 0;
}
DTYPE genrand_real1() // [0,1]
{ return (DTYPE) sfmt_genrand_real1(sfmt); }
DTYPE genrand_real2() // [0,1)
{ return (DTYPE) sfmt_genrand_real2(sfmt); }
DTYPE genrand_normal() // Box-Muller
{ return (DTYPE) sqrt(- 2.0 * log(1.0 - sfmt_genrand_real2(sfmt)))
* cos(2.0 * M_PI * (1.0 - sfmt_genrand_real2(sfmt))); }
//================================================================
struct neuron{
DTYPE *output;
char axonFamily;
DTYPE *weights;
DTYPE *weightGrad;
DTYPE *weightHistoricalGrad;
DTYPE *bias;
DTYPE *biasGrad;
DTYPE *biasHistoricalGrad;
int noOfInputs;
DTYPE *inputs;
DTYPE *error;
};
struct layer {
int noOfNeurons;
struct neuron *neurons;
};
struct neuralNet {
int noOfInputs;
DTYPE *inputs;
DTYPE *outputs;
int noOfLayers;
struct layer *layers;
int noOfBatchChanges;
int noOfIterations;
} theNet;
DTYPE weightInit() { return 0.3 * genrand_normal(0); }
DTYPE biasInit() { return 0.1; }
DTYPE getRand() { return 2.0 * M_PI * (2.0 * genrand_real1(0) - 1.0); }
int totalNoOfNeurons;
int totalNoOfWeights;
struct neuron *netNeurons;
DTYPE *netInputs;
DTYPE *netNeuronOutputs;
DTYPE *netErrors;
struct layer *netLayers;
DTYPE *netWeights;
DTYPE *netWeightGrad;
DTYPE *netWeightHistoricalGrad;
DTYPE *netBiases;
DTYPE *netBiasGrad;
DTYPE *netBiasHistoricalGrad;
void createNet(int noOfLayers, int *noOfNeurons, int *noOfInputs, char *axonFamilies, int initWeights) {
int i, j, counter, counter2, counter3, counter4;
theNet.noOfLayers = noOfLayers;
netLayers = (struct layer *) malloc(sizeof(struct layer) * theNet.noOfLayers);
theNet.layers = netLayers;
theNet.noOfInputs = noOfInputs[0];
netInputs = (DTYPE *) malloc(sizeof(DTYPE) * theNet.noOfInputs);
theNet.inputs = netInputs;
totalNoOfNeurons = 0;
for(i = 0; i < theNet.noOfLayers; i++) {
totalNoOfNeurons += noOfNeurons[i];
}
netNeurons = (struct neuron *) malloc(sizeof(struct neuron) * totalNoOfNeurons);
netNeuronOutputs = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
netErrors = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
netBiases = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
netBiasGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
netBiasHistoricalGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
for(i = 0; i < totalNoOfNeurons; i++) { netNeuronOutputs[i] = 0.0; }
totalNoOfWeights = 0;
for(i = 0; i < theNet.noOfLayers; i++) {
totalNoOfWeights += noOfInputs[i] * noOfNeurons[i];
}
netWeights = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfWeights);
netWeightGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfWeights);
netWeightHistoricalGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfWeights);
counter = 0;
counter2 = 0;
counter3 = 0;
counter4 = 0;
for(i = 0; i < theNet.noOfLayers; i++) {
for(j = 0; j < noOfNeurons[i]; j++) {
if(i == theNet.noOfLayers-1 && j == 0) { // beginning of the output layer
theNet.outputs = &netNeuronOutputs[counter];
}
netNeurons[counter].output = &netNeuronOutputs[counter];
netNeurons[counter].axonFamily = axonFamilies[i];
netNeurons[counter].weights = &netWeights[counter2];
netNeurons[counter].weightGrad = &netWeightGrad[counter2];
netNeurons[counter].weightHistoricalGrad = &netWeightHistoricalGrad[counter2];
netNeurons[counter].bias = &netBiases[counter];
netNeurons[counter].biasGrad = &netBiasGrad[counter];
netNeurons[counter].biasHistoricalGrad = &netBiasHistoricalGrad[counter];
netNeurons[counter].noOfInputs = noOfInputs[i];
if (i == 0) {
netNeurons[counter].inputs = netInputs;
}
else {
netNeurons[counter].inputs = &netNeuronOutputs[counter3];
}
netNeurons[counter].error = &netErrors[counter];
counter2 += noOfInputs[i];
counter++;
}
netLayers[i].noOfNeurons = noOfNeurons[i];
netLayers[i].neurons = &netNeurons[counter4];
if (i > 0) {
counter3 += noOfNeurons[i-1];
}
counter4 += noOfNeurons[i];
}
if (initWeights == 1) {
for(i = 0; i < totalNoOfNeurons; i++) { netBiases[i] = biasInit(); }
for(i = 0; i < totalNoOfNeurons; i++) { netBiasHistoricalGrad[i] = 0.0; }
for(i = 0; i < totalNoOfWeights; i++) { netWeights[i] = weightInit(); }
for(i = 0; i < totalNoOfWeights; i++) { netWeightHistoricalGrad[i] = 0.0; }
}
for(i = 0; i < totalNoOfNeurons; i++) { netBiasGrad[i] = 0.0; }
for(i = 0; i < totalNoOfWeights; i++) { netWeightGrad[i] = 0.0; }
theNet.noOfBatchChanges = 0;
theNet.noOfIterations = 0;
}
void feedNetInputs(DTYPE *inputs) {
int i;
for ( i = 0; i < theNet.noOfInputs; i++ ) {
netInputs[i] = inputs[i];
}
}
void updateNeuronOutput(struct neuron *myNeuron) {
DTYPE activation = 0.0;
int i;
for (i = 0; i < myNeuron->noOfInputs; i++) {
activation += myNeuron->inputs[i] * myNeuron->weights[i];
}
activation += *(myNeuron->bias);
DTYPE temp;
switch (myNeuron->axonFamily) {
case 'g': // sigmoid
temp = - activation;
if (temp > 45.0) { *(myNeuron->output) = 0.0; }
else if (temp < -45.0) { *(myNeuron->output) = 1.0; }
else { *(myNeuron->output) = 1.0 / (1.0 + exp(temp)); }
break;
case 't': // tanh
temp = - activation;
if (temp > 45.0) { *(myNeuron->output) = -1.0; }
else if (temp < -45.0) { *(myNeuron->output) = 1.0; }
else { *(myNeuron->output) = (2.0 / (1.0 + exp(temp))) - 1.0; }
break;
case 'r': // ReLU
if (activation > 0.0) { *(myNeuron->output) = activation; }
else { *(myNeuron->output) = 0.0; }
break;
default:
break;
}
}
void updateNetOutput() {
int i, j;
for(i = 0; i < theNet.noOfLayers; i++) {
for(j = 0; j < theNet.layers[i].noOfNeurons; j++) {
updateNeuronOutput(&(theNet.layers[i].neurons[j]));
}
}
}
DTYPE derivative(struct neuron *myNeuron) {
DTYPE temp;
switch (myNeuron->axonFamily) {
case 'g': // sigmoid
temp = *(myNeuron->output) * (1.0 - *(myNeuron->output)); break;
case 't': // tanh
temp = 1.0 - 0.5 * *(myNeuron->output) * *(myNeuron->output); break;
case 'r': // ReLU
if (*(myNeuron->output) > 0.0) { temp = 1.0; }
else { temp = 0.0; }
break;
default:
temp = 0; break;
}
return temp;
}
void trainNet(DTYPE *outputTargets) {
int i,j,k;
DTYPE temp;
struct layer *currLayer, *nextLayer;
// calculate errors
for(i = theNet.noOfLayers - 1; i >= 0; i--) {
currLayer = &theNet.layers[i];
if (i == theNet.noOfLayers - 1 ) { // output layer
for (j = 0; j < currLayer->noOfNeurons; j++ ) {
*(currLayer->neurons[j].error)
= derivative(&currLayer->neurons[j]) * (outputTargets[j] - *(currLayer->neurons[j].output));
}
}
else { // other layers
nextLayer = &theNet.layers[i+1];
for (j = 0; j < currLayer->noOfNeurons; j++ ) {
temp = 0.0;
for (k = 0; k < nextLayer->noOfNeurons; k++ ) {
temp += *(nextLayer->neurons[k].error) * nextLayer->neurons[k].weights[j];
}
*(currLayer->neurons[j].error) = derivative(&currLayer->neurons[j]) * temp;
}
}
}
// update weights n biases
for(i = theNet.noOfLayers - 1; i >= 0; i--) {
currLayer = &theNet.layers[i];
for (j = 0; j < currLayer->noOfNeurons; j++) {
*(currLayer->neurons[j].biasGrad) += *(currLayer->neurons[j].error);
// weights
for(k = 0; k < currLayer->neurons[j].noOfInputs; k++) {
currLayer->neurons[j].weightGrad[k]
+= *(currLayer->neurons[j].error) * currLayer->neurons[j].inputs[k];
}
}
}
theNet.noOfBatchChanges++;
}
void setGradient(DTYPE alpha) {
int i,j,k;
struct layer *currLayer;
for(i = theNet.noOfLayers - 1; i >= 0; i--) {
currLayer = &theNet.layers[i];
for (j = 0; j < currLayer->noOfNeurons; j++) {
// biases
*(currLayer->neurons[j].biasGrad) /= theNet.noOfBatchChanges;
if (theNet.noOfIterations == 0) {
*(currLayer->neurons[j].biasHistoricalGrad)
= *(currLayer->neurons[j].biasGrad) * *(currLayer->neurons[j].biasGrad);
}
else {
*(currLayer->neurons[j].biasHistoricalGrad)
= alpha * *(currLayer->neurons[j].biasHistoricalGrad)
+ (1.0 - alpha) * *(currLayer->neurons[j].biasGrad) * *(currLayer->neurons[j].biasGrad);
}
// weights
for (k = 0; k < currLayer->neurons[j].noOfInputs; k++) {
currLayer->neurons[j].weightGrad[k] /= theNet.noOfBatchChanges;
if (theNet.noOfIterations == 0) {
currLayer->neurons[j].weightHistoricalGrad[k]
= currLayer->neurons[j].weightGrad[k] * currLayer->neurons[j].weightGrad[k];
}
else {
currLayer->neurons[j].weightHistoricalGrad[k]
= alpha * currLayer->neurons[j].weightHistoricalGrad[k]
+ (1.0 - alpha) * currLayer->neurons[j].weightGrad[k] * currLayer->neurons[j].weightGrad[k];
}
}
}
}
}
void update(DTYPE stepsize) {
int i,j,k;
struct layer *currLayer;
DTYPE fudgeFactor = 1e-8;
for(i = theNet.noOfLayers - 1; i >= 0; i--) {
currLayer = &theNet.layers[i];
for (j = 0; j < currLayer->noOfNeurons; j++) {
// biases
*(currLayer->neurons[j].bias)
+= stepsize * (*(currLayer->neurons[j].biasGrad)
/ sqrt(fudgeFactor + *(currLayer->neurons[j].biasHistoricalGrad)));
*(currLayer->neurons[j].biasGrad) = 0.0;
// weights
for (k = 0; k < currLayer->neurons[j].noOfInputs; k++) {
currLayer->neurons[j].weights[k]
+= stepsize * (currLayer->neurons[j].weightGrad[k]
/ (sqrt(fudgeFactor + currLayer->neurons[j].weightHistoricalGrad[k])));
currLayer->neurons[j].weightGrad[k] = 0.0;
}
}
}
theNet.noOfBatchChanges = 0;
theNet.noOfIterations++;
}
DTYPE *getOutputs() {
return theNet.outputs;
}
int main() {
if (initRand(0)) return 1;
/* determine layer paramaters */
int noOfLayers = 4; // input layer excluded
int noOfNeurons[] = {20,20,20,2};
int noOfInputs[] = {2,20,20,20};
char axonFamilies[] = {'r','r','r','t'};
DTYPE *inputs;
DTYPE *outputTargets;
inputs = (DTYPE *) malloc(sizeof(DTYPE) * noOfInputs[0]);
outputTargets = (DTYPE *) malloc(sizeof(DTYPE) * noOfNeurons[noOfLayers - 1]);
createNet(noOfLayers, noOfNeurons, noOfInputs, axonFamilies, 1);
/* train it using batch method */
int i;
DTYPE tempTotal;
int counter = 0;
for(i = 0; i < 10000000; i++) {
inputs[0] = getRand();
inputs[1] = getRand();
tempTotal = inputs[0] + inputs[1];
feedNetInputs(inputs);
updateNetOutput();
outputTargets[0] = (DTYPE) sin(tempTotal);
outputTargets[1] = (DTYPE) cos(tempTotal);
trainNet(outputTargets);
counter++;
if (counter == 50) { setGradient(0.9); update(0.001); counter = 0; }
}
/* test it */
DTYPE *outputs;
printf("Sin Target \t Output \t Cos Target \t Output\n");
printf("---------- \t -------- \t ---------- \t --------\n");
for(i = 0; i < 50; i++) {
inputs[0] = getRand();
inputs[1] = getRand();
tempTotal = inputs[0] + inputs[1];
feedNetInputs(inputs);
updateNetOutput();
outputs = getOutputs();
printf( "%f \t %f \t %f \t %f \n", sin(tempTotal), outputs[0], cos(tempTotal), outputs[1]);
}
return 0;
}
実行結果は以下のような感じになると思います。
$ ./mlp
Sin Target Output Cos Target Output
---------- -------- ---------- --------
0.640082 0.746855 -0.768307 -0.812494
0.069101 0.079599 -0.997610 -0.988040
-0.985831 -0.984480 -0.167741 -0.157381
0.521136 0.532716 0.853474 0.864914
-0.097435 -0.109023 0.995242 0.994873
-0.114213 -0.112507 0.993456 0.992926
-0.917489 -0.925177 -0.397762 -0.380658
-0.046450 -0.051414 0.998921 0.996574
-0.692018 -0.679392 -0.721880 -0.727636
0.409927 0.424983 -0.912118 -0.924427
-0.916550 -0.913544 0.399921 0.423543
-0.824584 -0.721065 0.565740 0.638374
0.023286 0.022497 0.999729 0.967530
-0.650853 -0.650613 0.759204 0.782016
0.972549 0.942901 -0.232696 -0.201588
0.014787 -0.000262 -0.999891 -0.962320
-0.969422 -0.968556 0.245399 0.293858
0.338011 0.355989 -0.941142 -0.948886
-0.981046 -0.981085 -0.193777 -0.182955
0.520488 0.611798 0.853869 0.867387
0.958533 0.969545 -0.284982 -0.269451
-0.978963 -0.977859 -0.204036 -0.193433
-0.905485 -0.877378 0.424378 0.436468
-0.438091 -0.431923 -0.898930 -0.885023
-0.985005 -0.986799 0.172525 0.187228
-0.400110 -0.335487 -0.916467 -0.900038
-0.248896 -0.262504 0.968530 0.979315
0.443793 0.505794 -0.896129 -0.895520
-0.182260 -0.181806 -0.983250 -0.978943
0.322289 0.332658 0.946641 0.947475
-0.291678 -0.259068 0.956516 0.947655
0.102971 0.131531 -0.994684 -0.988836
-0.779955 -0.766913 0.625836 0.646452
0.906765 0.914458 0.421636 0.311500
0.877212 0.894052 0.480104 0.521253
-0.653435 -0.650393 0.756982 0.782935
-0.727792 -0.727932 0.685798 0.689874
-0.085042 -0.070495 -0.996377 -0.988692
0.124316 0.127393 0.992243 0.992283
-0.996733 -0.990036 -0.080773 -0.043595
-0.472917 -0.456472 -0.881107 -0.884031
-0.594023 -0.583654 0.804448 0.834040
0.980313 0.967607 0.197448 0.239313
-0.590370 -0.562181 0.807133 0.805516
-0.289410 -0.270091 -0.957205 -0.961424
-0.555648 -0.559167 0.831418 0.860375
0.642193 0.643480 0.766543 0.785267
0.486129 0.552793 -0.873887 -0.881352
-0.491824 -0.466695 0.870694 0.875552
-0.974766 -0.975019 -0.223228 -0.213617