Why BatchNormalization fails in Keras
I have played around with it and managed to get a working configuration.
Necessary changes:


axis=1 for BatchNormalization. Normally, this should be used when data_format="channels_first", but it worked for me despite TensorFlow backend
Standardize the input:     

X_train = X_train * 2.0 - 1.0; X_test = X_test * 2.0 - 1.0
Avoid batchnorm after the last convolution in res_unit.


The working solution looks like this:

# build a model                                                                                                                                                                                                                                                       
def res_unit(x, prefix=''):                                                                                                                                                                                                                                           
    x_shortcut = x                                                                                                                                                                                                                                                    

    x = Conv2D(16, (1, 1), padding='same', name='{}/{}'.format(prefix, 'conv1'))(x)                                                                                                                                                                                   
    x = BatchNormalization(name='{}/{}/bn'.format(prefix, 'conv1'), axis=1)(x)                                                                                                                                                                                        
    x = Activation(activation='relu', name='{}/{}/relu'.format(prefix, 'conv1'))(x)                                                                                                                                                                                   

    x = Conv2D(16, (3, 3), padding='same', name='{}/{}'.format(prefix, 'conv2'))(x)                                                                                                                                                                                   
    x = BatchNormalization(name='{}/{}/bn'.format(prefix, 'conv2'), axis=1)(x)                                                                                                                                                                                        
    x = Activation(activation='relu', name='{}/{}/relu'.format(prefix, 'conv2'))(x)                                                                                                                                                                                   

    x = Conv2D(32, (1, 1), padding='same', name='{}/{}'.format(prefix, 'conv3'))(x)                                                                                                                                                                                   
    x = Add(name='{}/add'.format(prefix))([x_shortcut, x])                                                                                                                                                                                                            
    x = Activation('relu', name='{}/relu'.format(prefix))(x)                                                                                                                                                                                                          

    return x                                                                                                                                                                                                                                                          


def main():                                                                                                                                                                                                                                                           
    (X_train, y_train), (X_test, y_test) = load_data()                                                                                                                                                                                                                
    X_train = X_train.astype('float32') / 255                                                                                                                                                                                                                         
    X_test = X_test.astype('float32') / 255                                                                                                                                                                                                                           
    X_train = X_train * 2.0 - 1.0                                                                                                                                                                                                                                     
    X_test = X_test * 2.0 - 1.0                                                                                                                                                                                                                                       
    y_train = keras_utils.to_categorical(y_train)                                                                                                                                                                                                                     
    y_test = keras_utils.to_categorical(y_test)                                                                                                                                                                                                                       

    print(X_train.min())                                                                                                                                                                                                                                              
    print(X_train.max())                                                                                                                                                                                                                                              

    X_input = Input((X_train.shape[1:]))                                                                                                                                                                                                                              

    X = Conv2D(32, (3, 3), padding='same', name='conv0')(X_input)                                                                                                                                                                                                     
    X = BatchNormalization(name='conv0/bn', axis=1)(X)                                                                                                                                                                                                                
    X = Activation('relu', name='conv0/relu')(X)                                                                                                                                                                                                                      

    X = res_unit(X, prefix='block1')                                                                                                                                                                                                                                  
    X = res_unit(X, prefix='block2')                                                                                                                                                                                                                                  
    X = res_unit(X, prefix='block3')                                                                                                                                                                                                                                  

    X = Flatten()(X)                                                                                                                                                                                                                                                  
    X = Dense(32)(X)                                                                                                                                                                                                                                                  
    X = Activation('relu')(X)                                                                                                                                                                                                                                         
    X = Dense(10)(X)                                                                                                                                                                                                                                                  

    X = Activation('softmax', name='softmax')(X)                                                                                                                                                                                                                      

    model = Model(inputs=X_input, outputs=X)                                                                                                                                                                                                                          
    model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])


I've got 59.07% validation accuracy with it. However, sometimes training does not converge, it feels like it is very sensitive to the random init.
Still don't know why making batchnorm work here is so hard...