如何更有效地计算二项式和?

问题描述

我必须按如下方式计算方程:

I must calculate an equation as follows:

给出k1,k2的位置.我正在使用MATLAB计算P.我认为我对上述方程式有一个正确的实现.但是，我的执行太慢了.我认为问题出在二项式系数上.从等式中，我可以有一种有效的方法来加快时间吗?谢谢大家

where k1,k2 are given. I am using MATLAB to compute P. I think I have a correct implementation for the above equation. However, my implementation is so slow. I think the issue is from binomial coefficient. From the equation, could I have an efficient way to speed up the time? Thank all.

对于k1=150; k2=150; D=200;，它需要 11.6秒

function main
warning ('off');
  function test_binom()
      k1=150; k2=150; D=200; P=0;
      for i=0:D-1
          for j=0:i
              if (i-j>k2||j>k1)
                  continue;
              end
              P=P+nchoosek(k1,j)*nchoosek(k2,i-j)/nchoosek((k1+k2),i);
          end
      end
  end
f = @()test_binom();
timeit(f)
end

更新:对于测量时间，我发现nchoosek是计算时间较长的原因.因此，我将函数重写如下

Update: For measure time, I found that nchoosek is the reason for large computational time. Hence, I rewrite the function as follows

function re=choose(n, k)
    if (k == 0)
        re=1;
    else
        re=(n * choose(n - 1, k - 1)) / k;
    end
end

现在，计算时间减少为0.25秒.有更好的办法吗?

Now, the computational time is reduced as 0.25 second. Is has any better way?

推荐答案

您可以对所有过程进行矢量化处理，这将使它变得超快，而无需使用mex.

You can vectorize all the process and it will make it super-fast without any need for mex.

首先是nchoosek函数:

function C = nCk(n,k)
% use smaller k if available
k(k>n/2) = n-k(k>n/2);
k = k(:);
kmat = ones(numel(k),1)*(1:max(n-k));
kmat = kmat.*bsxfun(@le,kmat,(n-k));
pw = bsxfun(@power,kmat,-1./(n-k));
pw(kmat==0) = 1;
kleft = ones(numel(k),1)*(min(k):n);
kleft = kleft.*bsxfun(@gt,kleft,k);
t = bsxfun(@times,kleft,prod(pw,2));
t (kleft==0) = 1;
C = prod(t,2);
end

然后beta和P计算:

function P = binomial_coefficient(k1,k2,D)
warning ('off','MATLAB:nchoosek:LargeCoefficient');
i_ind = nonzeros(triu(ones(D,1)*(1:D)))-1;
j_ind = nonzeros(tril(ones(D,1)*(1:D+1)).')-1;
valid = ~(i_ind-j_ind>=k2 | j_ind>=k1);
i_ind = i_ind(valid);
j_ind = j_ind(valid);
beta = @(ii,jj) nCk(k1,jj).*nCk(k2,ii-jj)./nCk((k1+k2),ii);
b = beta(i_ind,j_ind);
P = sum(b(:));
end

，执行时间从 10.674 减少到 0.49696 秒.

and execution time drops from 10.674 to 0.49696 seconds.

采用@ rahnema1的想法，我设法使用一个表来进行所有唯一的nCk计算，从而使这一过程变得更快，因此，一次完成的操作都不多.使用上面相同的nCk函数，这就是新的binomial_coefficient函数的外观:

Taking the idea of @rahnema1, I managed to make this even faster, using a table for all unique nCk computations, so none of them will be done more than once. Using the same nCk function from above, this is how the new binomial_coefficient function will look:

function P = binomial_coefficient(k1,k2,D)
warning ('off','MATLAB:nchoosek:LargeCoefficient');
i_ind = nonzeros(triu(ones(D,1)*(1:D)))-1;
j_ind = nonzeros(tril(ones(D,1)*(1:D+1)).')-1;
valid = ~(i_ind-j_ind>=k2 | j_ind>=k1);
i_ind = i_ind(valid);
j_ind = j_ind(valid);
ni = numel(i_ind);
all_n = repelem([k1; k2; k1+k2],ni); % all n's to be used in thier order
all_k = [j_ind; i_ind-j_ind; i_ind]; % all k's to be used in thier order
% get all unique sets of 'n' and 'k':
sets_tbl = unique([all_n all_k],'rows');
uq_n = unique(sets_tbl(:,1));
nCk_tbl = zeros([max(all_n) max(all_k)+1]);
% compute all the needed values of nCk:
for s = 1:numel(uq_n)
    curret_n = uq_n(s);
    curret_k = sets_tbl(sets_tbl(:,1)==curret_n,2);
    nCk_tbl(curret_n,curret_k+1) = nCk(curret_n,curret_k).';
end
beta = @(ii,jj) nCk_tbl(k1,jj+1).*nCk_tbl(k2,ii-jj+1)./nCk_tbl((k1+k2),ii+1);
b = beta(i_ind,j_ind);
P = sum(b(:));
end

现在，当仅需要 0.01212 秒运行时，它不仅是超快速的代码，而且是飞行通话超快代码！

and now, when it takes only 0.01212 second to run, it's not just super-fast code, it's a flying-talking-super-fast code!

这篇关于如何更有效地计算二项式和?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！