如何使用Thrust进行外积降维操作

Question

我是一位有用的助手，可以为您翻译文本。

我有一个二维数组$a$，存储在具有维度$N$和$m$的索引$(p,i)$的设备向量中。

我想计算：

$$s_{ij} = \sum \limits_{p=1}^{N} a_{p,i} a_{p,j}$$

for $i,j=1,...,m.$ 

有没有一种简单的方式使用Thrust来实现这个？

上面的代码是LaTeX。在C ++中，它会类似于：

Matrix A(N,m); // filled with data

Matrix S(m,m);

for (int i=0; i <m;++i)
  for (int j=0; j <m;++j) 
   { 
     S(i,j)=0;
     for (int p=0; p < N; ++p)
        S(i,j) += A(p,i)*A(p,j);

   }

- Mark Joshi

---

你想要什么？我无法通过这些符号完全理解。 - dreamcrash

看着这个问题，我感觉自己刚刚在智商测试中失败了。你能提供串行代码吗？ - Roger Dahl

我认为在thrust中实现外积或Kronecker积并没有直接的方法。 - talonmies

我认为没有一个好的方法来解决这个问题。你可以尝试使用一些高级迭代表达式生成M X N的乘积，然后再通过某种方式调用reduce_by_key函数。但是性能会很差，因为它无法捕获外部两个嵌套for循环中的重复利用。更直接的方法是将众多专门的N-body CUDA内核之一改编为此问题的解决方案。 - Jared Hoberock

Answer 1

如果我没理解错的话，似乎是这样的：

$$s_{ij} = \sum_{p=1}^{N} a_{p,i} a_{p,j} = \sum_{p=1}^{N} a^T_{i,p} a_{p,j}$$

因此

$$S = A^T A$$

通常的矩阵乘积。我认为你不能使用Thrust来完成这个任务。但是你可以使用CUBLAS或Arrayfire轻松地进行矩阵乘法（实际上我认为Arrayfire在内部使用CUBLAS）。但请记住，这些库将矩阵按列主序存储（类似于Fortran）。

- user1545642

---

是的，一个人可以这样做，但我的主要担心是我需要两份矩阵副本。想想N非常大，M相当小。 - Mark Joshi