rgb2yuv_matrix = np.array([[0.299, 0.587, 0.114], [-0.1473, -0.28886, 0.436],[0.615, -0.51499, 0.10001]])
for i in range(n_train):
    for j in range(32):
        for k in range(32):
            # Note here that X_train[i, j, k] is a 3 element array
            X_train[i,j,k] = np.dot(rgb2yuv_matrix, X_train[i,j,k])

import numpy as np
from scipy.stats import pearsonr

x = np.array([[[1,2,3,4],
               [5,6,7,8]],
              [[11,22,23,24],
               [25,26,27,28]]])


i,j,k = x.shape

y = np.array([[[31,32,33,34],
               [35,36,37,38]],
              [[41,42,43,44],
               [45,46,47,48]]])



xx = np.row_stack(np.dstack(x))
yy = np.row_stack(np.dstack(y))

results = []

for a, b in zip(xx,yy):
    r_sq, p_val = pearsonr(a, b)
    results.append(r_sq)

results = np.array(results).reshape(j,k)

print results

[[ 1.  1.  1.  1.]
 [ 1.  1.  1.  1.]]