Question: I want to perform Matthews correlation coefficient (MCC) in sklearn to find the correlation between different features (boolean vectors) in a 2D numpyarray. What I have done so far is to loop through each column and find correlation value between features one by one. Here is my code: from sklearn.metrics import matthews_corrcoef import numpy as np X = np.array([[1, 0, 0, 0, 0], [1, 0, 0, 1, 0], [1, 0, 0, 0, 1], [1, 1, 0, 0, 0], [1, 1, 0, 1, 0], [1, 1, 0, 0, 1], [1, 0, 1, 0, 0], [1, 0, 1, 1, 0], [1, 0, 1, 0, 1], [1, 0, 0, 0, 0]]) n_sample, n_feature = X.shape rff_all = [] for i in range(n_feature): for j in range(i + 1, n_feature): coeff_f_f = abs(matthews_corrcoef(X[:, i], X[:, j])) rff_all.append(coeff_f_f) rff = np.mean(rff_all)  As I have a huge dimension of 2D numpyarray, it seems to be really slow and impractical. What is the most efficient way to perform this kind of operation simultaneously without using the loops? Edit: I then come up with this idea but it is still pretty slow. from more_itertools import distinct_combinations all_c = [] for item in distinct_combinations(np.arange(X.shape[1]), r=2): c = matthews_corrcoef(X[:, item][:, 0], X[:, item][:, 1]) all_c.append(abs(c)) 

Answer: You can use numba to speed up the computing, e.g: import numba import numpy as np @numba.njit def _fill_cm(m, c1, c2): m[:] = 0 for a, b in zip(c1, c2): m[a, b] += 1 @numba.njit def mcc(confusion_matrix): # https://stackoverflow.com/a/56875660/992687 tp = confusion_matrix[0, 0] tn = confusion_matrix[1, 1] fp = confusion_matrix[1, 0] fn = confusion_matrix[0, 1] x = (tp + fp) * (tp + fn) * (tn + fp) * (tn + fn) return ((tp * tn) - (fp * fn)) / sqrt(x + 1e-6) @numba.njit def get_all_mcc_numba(X): rows, columns = X.shape confusion_matrix = np.zeros((2, 2), dtype="float32") out = [] for i in range(columns): c1 = X[:, i] for j in range(i + 1, columns): # make confusion matrix c2 = X[:, j] _fill_cm(confusion_matrix, c1, c2) out.append(abs(mcc(confusion_matrix))) return out  Benchmark: from timeit import timeit from math import sqrt import numba import numpy as np from sklearn.metrics import matthews_corrcoef def get_all_mcc_normal(X): n_sample, n_feature = X.shape rff_all = [] for i in range(n_feature): for j in range(i + 1, n_feature): coeff_f_f = abs(matthews_corrcoef(X[:, i], X[:, j])) rff_all.append(coeff_f_f) return rff_all @numba.njit def _fill_cm(m, c1, c2): m[:] = 0 for a, b in zip(c1, c2): m[a, b] += 1 @numba.njit def mcc(confusion_matrix): # https://stackoverflow.com/a/56875660/992687 tp = confusion_matrix[0, 0] tn = confusion_matrix[1, 1] fp = confusion_matrix[1, 0] fn = confusion_matrix[0, 1] x = (tp + fp) * (tp + fn) * (tn + fp) * (tn + fn) return ((tp * tn) - (fp * fn)) / sqrt(x + 1e-6) @numba.njit def get_all_mcc_numba(X): rows, columns = X.shape confusion_matrix = np.zeros((2, 2), dtype="float32") out = [] for i in range(columns): c1 = X[:, i] for j in range(i + 1, columns): # make confusion matrix c2 = X[:, j] _fill_cm(confusion_matrix, c1, c2) out.append(abs(mcc(confusion_matrix))) return out # make 2000 x 100 0-1 matrix: np.random.seed(42) X = np.random.randint(low=0, high=2, size=(2000, 100), dtype="uint8") assert np.allclose(get_all_mcc_normal(X), get_all_mcc_numba(X)) t_normal = timeit("get_all_mcc_normal(X)", number=1, globals=globals()) t_numba = timeit("get_all_mcc_numba(X)", number=1, globals=globals()) print(f"{t_normal=}") print(f"{t_numba=}")  Prints on my computer (AMD 5700x): t_normal=4.352220230008243 t_numba=0.008588693017372862   You can further speed up the computing using numba multiprocessing: @numba.njit(parallel=True) def get_all_mcc_numba_parallel(X): rows, columns = X.shape num_cpus = numba.get_num_threads() # for each thread, allocate array for confusion matrix cms = [] for i in range(num_cpus): cms.append(np.empty((2, 2), dtype="float32")) out = np.empty(shape=(columns, columns), dtype="float32") # make indexes for each thread thread_column_idxs = np.array_split(np.arange(columns), num_cpus) for i in range(columns): c1 = X[:, i] for thread_idx in numba.prange(num_cpus): for j in thread_column_idxs[thread_idx]: if j < i + 1: continue c2 = X[:, j] _fill_cm(cms[thread_idx], c1, c2) out[i, j] = abs(mcc(cms[thread_idx])) out2 = [] for i in range(columns): for j in range(i + 1, columns): out2.append(out[i, j]) return out2  Benchmark (using only numba/numba+parallel versions): import perfplot np.random.seed(42) X = np.random.randint(low=0, high=2, size=(2000, 100), dtype="uint8") assert np.allclose(get_all_mcc_normal(X), get_all_mcc_numba(X)) assert np.allclose(get_all_mcc_numba(X), get_all_mcc_numba_parallel(X)) perfplot.show( setup=lambda n: np.random.randint(low=0, high=2, size=(2000, n), dtype="uint8"), kernels=[ get_all_mcc_numba, get_all_mcc_numba_parallel, ], labels=["single", "parallel"], n_range=[100, 250, 500, 750, 1000], xlabel="2000 rows x N columns", logx=True, logy=True, equality_check=np.allclose, )  Creates this graph:   For matrix with 10_000 rows and 10_000 columns: t_numba=782.5051075649972 (~13 minutes) t_numba_parallel=215.43611432801117 (~3.5 minutes)