如何加速并扩展处理20000×20000指示矩阵的Python函数？

我有一个处理形状为(20000×20000)指示矩阵的函数，需要执行4亿次（20000×20000）。注意传入函数的indicator_Matrix必须是Pandas DataFrame格式，实际场景中该DataFrame包含时间索引和整数列（此处已简化问题）。

Pandas实现 #

indicator_Matrix = pd.DataFrame(np.random.randint(0,2,[20000,20000]))
def operations(indicator_Matrix):
   s = indicator_Matrix.sum(axis=1)
   d = indicator_Matrix.div(s,axis=0)
   res = d[d>0].mean(axis=0)
   return res.iloc[-1]

我尝试用Numpy改写优化，但运行耗时仍极长；也试过concurrent.futures.ThreadPoolExecutor，但相比列表推导式提升不大。

Numpy实现 #

indicator_Matrix = pd.DataFrame(np.random.randint(0,2,[20000,20000]))
def operations(indicator_Matrix):
   s = indicator_Matrix.to_numpy().sum(axis=1)
   d = (indicator_Matrix.to_numpy().T / s).T
   d = pd.DataFrame(d, index = indicator_Matrix.index, columns = indicator_Matrix.columns)
   res = d[d>0].mean(axis=0)
   return res.iloc[-1]

output = [operations(indicator_Matrix) for i in range(0,20000**2)]

注：将d重新转为DataFrame是因为需要获取列均值并保留最后一列的均值（使用.iloc[-1]）。d[d>0].mean(axis=0)返回列均值，示例如下：

2478    1.0
0       1.0

更新：目前仍受困于此问题，想了解在本地桌面使用cudf、CuPy等GPU库是否能带来性能改善。

内容的提问来源于stack exchange，提问作者user1769197