def _pearson_corr(mat_X, mat_Y):
    """Pearson's correlation between every columns in mat_X and mat_Y
    
    Args
    ----
        mat_X (N,M1): np.ndarray
        
        mat_Y (N,M2): np.ndarray

    Returns
    -------
        mat_corr: (M1,M2): np.ndarray
            Correlation matrix
    """
    
    # Reshape 
    if len(mat_X.shape)==1:
        mat_X = mat_X.reshape([-1,1])
    if len(mat_Y.shape)==1:
        mat_Y = mat_Y.reshape([-1,1])
        
    mat_X = (mat_X-mat_X.mean(axis=0))/mat_X.std(axis=0).clip(min=1e-8)
    mat_Y = (mat_Y-mat_Y.mean(axis=0))/mat_Y.std(axis=0).clip(min=1e-8)
    mat_corr = mat_X.T.dot(mat_Y)/mat_X.shape[0]
    
    if mat_corr.shape[1]==1:
        return mat_corr.reshape([-1])
    else:
        return mat_corr