Using rbf_kernel with two 2D numpy arrays causes Python 3 kernel to die

I am trying to calculate the maximum mean discrepancy between two datasets, X, Y, where the entries of X, Y are of type numpy.float64. Part of this calculation involves computing all pairwise dot products between X and X using some kernel function. In this case, I want to use sklearn.metrics.pairwise.rbf_kernel from sklearn. Using Jupyter, the following code causes the Python 3 kernel to die after ~30 seconds, then prompts me to restart it (I show two different methods that causes the kernel to die):

import sklearn as sklearn
from sklearn.metrics.pairwise import pairwise_kernels

def mmd(X, Y):

# X.shape -> (922315, 14)
# type(X) -> <class 'numpy.ndarray'>

# Kernel dies here
Kxx = sklearn.metrics.pairwise.rbf_kernel(X, X)

# Kernel also dies using this
Kxx = pairwise_kernels(X, X, metric='rbf')
.
.
.


After more debugging, I found that Python is Unable to allocate 6.18 TiB for an array with shape (921819, 921819) and data type float64. How does one get around this issue? Do I need to reduce the number of data points in X? If so, what is the proper way to do so?