SGD learning in CNN gets stuck when using a max pooling layer (x-post from DataScience) I'm working on a CNN library for a university project and I'm having some trouble implementing the backpropagation through the max pooling layer.

Please note that the whole thing was built from scratch following the
  neuralnetworksanddeeplearning book and other resources, so I'm aware there are ready-to-use libraries out there, as well as more optimized algorithms to perform certain operations (like the FFT convolution algorithm), but the point of this thing was to understand how a network works and to implement everything without using 3rd party libraries.

I'm quite sure the fully connected layers work fine, as I've tried training a network with the MNIST dataset and it quickly gets to over 96% in less than 60 epochs, just like the one in the book mentioned above.
I went on and implemented the convolutional layers, and they seem to be working fine too, in fact if I try to have a network like this one (sigmoid activation on all the layers):


*

*Convolution: 28*28*1 inputs, 20 5*5*1 kernels

*Fully Connected: 24*24*20 inputs, 100 outputs

*Fully Connected: 100 inputs, 10 outputs, cross-entropy cost


And here's what I got: 
As you can see, the network quickly reaches over 96% accuracy in just 3 epochs.
But, as soon as I try to put a 2*2 pooling layer between the convolutional and fully connected layer (just like Michael Nielsen did in the 6th chapter of his book), here's what I get:

The network reaches 30% and then just gets stuck there.
Update - here's a plot of the cross-entropy cost using the pooling layer:

The cost gets stuck there at around 2.52, and the accuracy stays at around 30%, as stated before.
I am aware that using a max pooling layer causes only the weights that contributed to the activated neurons to be updated, but the thing is that here I'm exactly replicating the network structure used by Michael Nielsen in his book, and he states that this network easily gets to over 98% of accuracy in less than 60 epochs (using the same training parameters).
Here instead it's clear that my network will never even get close to 98%, even after a 100 epochs.
Update #2 - here's another cost plot over 200 epochs:

I've checked and everything seems to be working fine (but of course I can be wrong, and I probably am, since this thing is not working as it should), the forward pooling layer takes the maximum value in each 2*2 window across all the input images (doing so for each depth layer), and during the backpropagation through the pooling layer I just upscale the input delta using the previous outputs from the convolutional layer, so that each delta goes to the pixel that had the maximum value in the convolutional output, while all the other pixels get 0.
Where should I start looking for the error? I mean, clearly something happens with the pooling layer, but the code looks fine to me, both in theory and with the actual Unit tests I've added to the project. What's the right approach to take here to investigate the issue? Are there probable explanations for this issue?
Here's the code I'm using, if anyone has .NET Core 2.x installed and wants to try the library out (The code to download/parse the MNIST dataset is already in the library):
((float[,] X, float[,] Y) training, (float[,] X, float[,] Y) test) = DataParser.LoadDatasets();
INeuralNetwork network = NetworkTrainer.NewNetwork(
    NetworkLayers.Convolutional((28, 28, 1), (5, 5), 20, ActivationFunctionType.Identity),
    NetworkLayers.Pooling((24, 24, 20), ActivationFunctionType.Sigmoid),
    NetworkLayers.FullyConnected(12 * 12 * 20, 100, ActivationFunctionType.Sigmoid),
    NetworkLayers.FullyConnected(100, 10, ActivationFunctionType.Sigmoid, CostFunctionType.CrossEntropy));
await NetworkTrainer.TrainNetworkAsync(network, (training.X, training.Y), 60, 10, null,
    new TestParameters(test, new Progress<BackpropagationProgressEventArgs>(p =>
    {
        Printf($"Epoch {p.Iteration}, cost: {p.Cost}, accuracy: {p.Accuracy}");
    })));

Thanks for your help!
 A: Build a unit test for your max-pooling layer.
Here's an example: Forwards: https://github.com/hughperkins/DeepCL/blob/master/test/testpoolingforward.cpp#L20
PoolingForward *poolingForward = PoolingForward::instanceForTest( cl, false, numPlanes, imageSize, poolingSize );
float data[] = { 1, 2, 5, 3,
                 3, 8, 4, 1,
                 3, 33, 14,23,
                 -1, -3.5f,37.4f,5
};
int outputNumElements = poolingForward->getOutputNumElements( batchSize );
int *selectors = new int[outputNumElements];
float *output = new float[outputNumElements];

poolingForward->forward( batchSize, data, selectors, output );

EXPECT_EQ( selectors[0], 3 );
EXPECT_EQ( selectors[1], 0 );
EXPECT_EQ( selectors[2], 1 );
EXPECT_EQ( selectors[3], 2 );

EXPECT_EQ( output[0], 8 );
EXPECT_EQ( output[1], 5 );
EXPECT_EQ( output[2], 33 );
EXPECT_EQ( output[3], 37.4f );

Backwards: https://github.com/hughperkins/DeepCL/blob/master/test/testpoolingbackward.cpp#L19
PoolingBackward *poolingBackprop = PoolingBackward::instanceForTest( cl, false, numPlanes, imageSize, poolingSize );
float errors[] = {
    3, 5,
    2, 9
};
int selectors[] = {
    2, 1,
    0, 3
};
float *errorsForUpstream = new float[ poolingBackprop->getInputNumElements( batchSize ) ];

poolingBackprop->backward( batchSize, errors, selectors, errorsForUpstream );

float expectedErrorsForUpstream[] = {
    0,0,0,5,
    3,0,0,0,
    2,0,0,0,
    0,0,0,9,
};
for( int i = 0; i < 16; i++ ) {
    ASSERT_EQ( expectedErrorsForUpstream[i], errorsForUpstream[i] );
}

