Thanks for the feedback. I have moved relu after max_pool2d and indeed was applying softmax on the output twice (didn’t realise CrossEntropyLoss() includes it). The updated script is here 1 however the script still gets 0.72 accuracy and takes around 6 minutes, so weirdly not much has chained.
I was purposefully only max-pooling twice because my feature-maps are already quite small (and do the same thing with Chainer).
Ross, the data is the same for all the frameworks (it comes from a common function, instead of using the library’s versions of it). If it’s easier here are the relevant code extracts:
def create_symbol():
    class SymbolModule(nn.Module):
        def __init__(self):
            super(SymbolModule, self).__init__()
            self.conv1 = nn.Conv2d(3, 50, kernel_size=(3, 3), padding=(1, 1))
            self.conv2 = nn.Conv2d(50, 50, kernel_size=(3, 3), padding=(1, 1))
            self.conv3 = nn.Conv2d(50, 100, kernel_size=(3, 3), padding=(1, 1))
            self.conv4 = nn.Conv2d(100, 100, kernel_size=(3, 3), padding=(1, 1))
            # feature map size is 8*8 by pooling
            self.fc1 = nn.Linear(100*8*8, 512)
            self.fc2 = nn.Linear(512, N_CLASSES)

        def forward(self, x):
            x = self.conv2(F.relu(self.conv1(x)))
            # Apply relu after max-pool
            x = F.relu(F.max_pool2d(x, kernel_size=(2, 2), stride=(2, 2)))
            x = F.dropout(x, 0.25)
            
            x = self.conv4(F.relu(self.conv3(x)))
            # Apply relu after max-pool
            x = F.relu(F.max_pool2d(x, kernel_size=(2, 2), stride=(2, 2)))
            x = F.dropout(x, 0.25)
                       
            x = x.view(-1, 100*8*8)   # reshape Variable
            x = F.dropout(F.relu(self.fc1(x)), 0.5)
            return self.fc2(x)
    return SymbolModule()

def init_model(m):
    # Implementation of momentum:
    # v = \rho * v + g \\
    # p = p - lr * v
    opt = optim.SGD(m.parameters(), lr=LR, momentum=MOMENTUM)
    # Combines softmax output with negative log-likelihood
    criterion = nn.CrossEntropyLoss()
    return opt, criterion

x_train, x_test, y_train, y_test = cifar_for_library(channel_first=True)
y_train = y_train.astype(np.int64)
y_test = y_test.astype(np.int64)

sym = create_symbol()
sym.cuda() # CUDA!
optimizer, criterion = init_model(sym)
sym.train()  
for j in range(EPOCHS):
    for data, target in yield_mb(x_train, y_train, BATCHSIZE, shuffle=True):
        # Get samples
        data = Variable(torch.FloatTensor(data).cuda())
        target = Variable(torch.LongTensor(target).cuda())
        # Init
        optimizer.zero_grad()
        # Forwards
        output = sym(data)
        # Loss
        loss = criterion(output, target)
        # Back-prop
        loss.backward()
        optimizer.step()