当A为稀疏矩阵时，如何加速矩阵乘法XAX^T的计算？

假设X为r行c列的稠密嵌套数组，A为c×c的稀疏矩阵（以{(i,j): value}形式存储）。当A的非零元素总数z较小时，以下Python伪代码效率足够；但当c较大且z大于c时，即便想近似加速XAX^T的三方乘积也缺乏可行思路。

from itertools import product

def naive_threeway_matmul(X, A):
  """
    X: (r,c) dense nested arrays
    A: (c,c) sparse matrix, stored as {(i,j): value}
    
    rtn: (r,r) dense nested arrays: X @ A @ X.transpose()
  """
  r = len(X)
  c = len(A)
  rtn = [[0]*r for _ in range(r)]
  for (u,v) in product(range(r), repeat=2):
    rtn[u][v] = sum(
      X[u][i] * value * X[v][j]
      for (i,j), value in A.items()
    )
  return rtn

目前已有复杂度为*O(r²z)的解决方案，通过先计算稀疏-稠密乘积再计算稠密-稠密乘积可得到O(r²c)*的方案（当z≈c时），但二者均为类立方复杂度。

请问当c较大且z≈c时，是否存在近似算法可将运行时间降至接近O(r²)？暂不考虑numpy、jax等加速器，后续可再优化代码。

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