6

我正在寻找一种快速计算的方法

(1:N)'*(1:N)

对于相当大的N。我觉得问题的对称性使得实际进行乘法和加法是浪费的。

4

4 回答 4

14

你为什么要这样做的问题真的很重要。

从理论上讲，其他答案中建议的三角形方法将节省您的操作。@jgmao 的答案在减少乘数方面特别有趣。

在实际意义上，CPU 操作的数量不再是编写快速代码时最小化的指标。当 CPU 操作如此之少时，内存带宽占主导地位，因此调整缓存感知访问模式是如何使其快速运行的方法。矩阵乘法代码的实现非常高效，因为它是一种常见的操作，并且每个值得一提的 BLAS 数值库的实现都将使用优化的访问模式以及 SIMD 计算。

即使您直接编写 C 并将运算次数减少到理论上的最小值，您可能仍然无法击败完整的矩阵乘法。这归结为找到与您的操作最匹配的数字原语。

综上所述，有一个比 DGEMM（矩阵乘法）更接近的 BLAS 操作。它被称为 DSYRK，即 rank-k 更新，它可以用于精确的A'*A. 我很久以前为此编写的 MEX 函数就在这里。我已经很久没有弄乱它了，但是当我第一次写它的时候它确实起作用了，而且实际上跑得比直接的快A'*A。

/* xtrx.c: calculates x'*x taking advantage of the symmetry.
Peter Boettcher <email removed>
Last modified: <Thu Jan 23 13:53:02 2003> */

#include "mex.h"

const double one = 1;
const double zero = 0;

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  double *x, *z;
  int i, j, mrows, ncols;

  if(nrhs!=1) mexErrMsgTxt("One input required.");

  x = mxGetPr(prhs[0]);
  mrows = mxGetM(prhs[0]);
  ncols = mxGetN(prhs[0]);

  plhs[0] = mxCreateDoubleMatrix(ncols,ncols, mxREAL);
  z = mxGetPr(plhs[0]);

  /* Call the FORTRAN BLAS routine for rank k update */
  dsyrk_("U", "T", &ncols, &mrows, &one, x, &mrows, &zero, z, &ncols);

  /* Result is in the upper triangle.  Copy it down the lower part */
  for(i=0; i<ncols; i++)
      for(j=i+1; j<ncols; j++)
          z[i*ncols + j] = z[j*ncols + i];
}
于 2013-10-01T20:22:43.660 回答
6

MATLAB 的矩阵乘法通常非常快，但这里有几种方法可以只得到上三角矩阵。它们比天真的计算v'*v（或使用在 BLAS 中调用更合适的对称秩 k 更新函数的 MEX 包装器，不足为奇！）要慢。无论如何，这里有一些仅限 MATLAB 的解决方案：

第一个使用线性索引：

% test vector
N = 1e3;
v = 1:N;

% compute upper triangle of product
[ii, jj] = find(triu(ones(N)));
upperMask = false(N,N);
upperMask(ii + N*(jj-1)) = true;
Mu = zeros(N);
Mu(upperMask) = v(ii).*v(jj); % other lines always the same computation

% validate
M = v'*v;
isequal(triu(M),Mu)

下一种方法也不会比天真的方法更快，但这是另一种计算下三角形的bsxfun解决方案：

Ml = bsxfun(@(x,y) [zeros(y-1,1); x(y:end)*y],v',v);

对于上三角形：

Mu = bsxfun(@(x,y) [x(1:y)*y; zeros(numel(x)-y,1)],v',v);
isequal(triu(M),Mu)

使用这种特殊情况（其中）的整个矩阵的另一种解决方案。这个速度实际上很接近。cumsumv=1:N

M = cumsum(repmat(v,[N 1]));

也许这些可以成为更好的起点。

于 2013-10-01T18:11:11.240 回答
5

这比 (1:N).'*(1:N) 快 3 倍，前提是int32结果是可以接受的（如果数字足够小，可以使用int16而不是 ，它会更快int32）：

N = 1000;
aux = int32(1:N);
result = bsxfun(@times,aux.',aux);

基准测试：

>> N = 1000; aux = int32(1:N); tic, for count = 1:1e2, bsxfun(@times,aux.',aux); end, toc
Elapsed time is 0.734992 seconds.

>> N = 1000; aux = 1:N; tic, for count = 1:1e2, aux.'*aux; end, toc
Elapsed time is 2.281784 seconds.

请注意，aux.'*aux不能用于aux = int32(1:N).

正如@DanielE.Shub 所指出的，如果需要将结果作为double矩阵，则必须进行最终转换，在这种情况下增益非常小：

>> N = 1000; aux = int32(1:N); tic, for count = 1:1e2, double(bsxfun(@times,aux.',aux)); end, toc
Elapsed time is 2.173059 seconds.
于 2013-10-03T16:00:15.927 回答
3

由于输入的特殊有序结构，考虑N=4的情况

(1:4)'*(1:4) = [1 2 3 4
                2 4 6 8
                3 6 9 12
                4 8 12 16]

您会发现第一行只是（1：N），从第二（j = 2）行开始，该行的值是前一行（j = 1）加上（1：N）。所以1.你不要做很多乘法。相反，您可以通过 N*N 加法生成它。2. 由于输出是对称的，所以只需要计算输出矩阵的一半。所以总计算量是 (N-1)+(N-2)+...+1 = N^2 / 2 次加法。

于 2013-10-01T19:31:46.847 回答