快速计算方法(1:N)'*(1:N)

Question

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

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

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

Answer 1

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

在理论意义上,其他答案中建议的三角形方法将为您节省操作.@ jgmao的答案在减少乘法时特别有趣.

在实际意义上,CPU操作的数量不再是编写快速代码时最小化的度量.当您拥有如此少的CPU操作时,内存带宽占主导地位,因此调整缓存感知访问模式是如何快速实现这一目标的.矩阵乘法代码非常有效地实现,因为它是如此常见的操作,并且值得使用的BLAS数值库的每个实现都将使用优化的访问模式和SIMD计算.

即使你写了直接C并将你的操作数减少到理论最小值,你可能仍然没有击败全矩阵乘法.这归结为找到与您的操作最匹配的数字原语.

总而言之,BLAS操作比DGEMM(矩阵乘法)更接近.它被称为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];
}

Answer 2

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]));

也许这些可以成为更好的事情的起点.

Answer 3

这比(1:N)快3倍.'*(1:N)提供的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.