MATLAB：高效生成多索引的大整数矩阵

Question

令 d 和 p 为两个整数。我需要生成一个大的整数矩阵 A，有 d 列和 N=nchoosek(d+p,p) 行。请注意，nchoosek(d+p,p) 会随着 d 和 p 快速增加，因此我可以快速生成 A 非常重要。 A 的行都是分量从 0 到 p 的多元指数，使得分量之和小于等于 p。这意味着，如果 d=3 且 p=3，则 A 是一个 [N=nchoosek(3+3,3)=20x3] 矩阵，其结构如下：

A=[0 0 0;
   1 0 0;
   0 1 0;
   0 0 1;
   2 0 0;
   1 1 0;
   1 0 1;
   0 2 0;
   0 1 1;
   0 0 2;
   3 0 0;
   2 1 0;
   2 0 1;
   1 2 0;
   1 1 1;       
   1 0 2;   
   0 3 0;      
   0 2 1;
   0 1 2;
   0 0 3]       

严格遵循我使用的行排序并不是必不可少的，尽管它会让我的生活更轻松（对于那些感兴趣的人，它被称为分级词典排序，如下所述： http://en.wikipedia.org/wiki/Monomial_order)。

如果您对这个奇怪矩阵的来源感到好奇，请告诉我！

Answer 1

使用nchoosek 和diff 的解决方案

以下解决方案基于 Mark Dickinson 的 this clever answer。

function degrees = monomialDegrees(numVars, maxDegree)
if numVars==1
    degrees = (0:maxDegree).';
    return;
end
degrees = cell(maxDegree+1,1);
k = numVars;
for n = 0:maxDegree
    dividers = flipud(nchoosek(1:(n+k-1), k-1));
    degrees{n+1} = [dividers(:,1), diff(dividers,1,2), (n+k)-dividers(:,end)]-1;
end
degrees = cell2mat(degrees);

您可以通过调用monomialDegrees(d,p)获取您的矩阵。

使用nchoosek 和accumarray/histc 的解决方案

这种方法基于以下思想：所有 k-multicombinations 与我们正在寻找的矩阵之间存在双射。多重组合给出了应添加条目的位置。例如，多重组合[1,1,1,1,3] 将映射到[4,0,1]，因为有四个1，一个3。这可以使用accumarray 或histc 进行转换。这是accumarray-方法：

function degrees = monomialDegrees(numVars, maxDegree)
if numVars==1
    degrees = (0:maxDegree).';
    return;
end
degrees = cell(maxDegree+1,1);
degrees{1} = zeros(1,numVars);
for n = 1:maxDegree
    pos = nmultichoosek(1:numVars, n);
    degrees{n+1} = accumarray([reshape((1:size(pos,1)).'*ones(1,n),[],1),pos(:)],1);
end
degrees = cell2mat(degrees);

这里是使用histc的替代方法：

function degrees = monomialDegrees(numVars, maxDegree)
if numVars==1
    degrees = (0:maxDegree).';
    return;
end
degrees = cell(maxDegree+1,1);
degrees(1:2) = {zeros(1,numVars); eye(numVars);};
for n = 2:maxDegree
    pos = nmultichoosek(1:numVars, n);
    degrees{n+1} = histc(pos.',1:numVars).';
end
degrees = cell2mat(degrees(1:maxDegree+1));

两者都使用以下函数来生成多重组合：

function combs = nmultichoosek(values, k)
if numel(values)==1
    n = values;
    combs = nchoosek(n+k-1,k);
else
    n = numel(values);
    combs = bsxfun(@minus, nchoosek(1:n+k-1,k), 0:k-1);
    combs = reshape(values(combs),[],k);
end

基准测试：

如果您的numVars 较低且maxDegree 较高，则对上述代码进行基准测试可以得出diff 解决方案更快。如果numVars 高于maxDegree，则histc 解决方案会更快。

旧方法：

这是dec2base 的Dennis' 方法的替代方法，该方法对最大基数有限制。它仍然比上述解决方案慢很多。

function degrees = monomialDegrees(numVars, maxDegree)
Cs = cell(1,numVars);
[Cs{:}] = ndgrid(0:maxDegree);
degrees = reshape(cat(maxDegree+1, Cs{:}),(maxDegree+1)^numVars,[]);
degrees = degrees(sum(degrees,2)<=maxDegree,:);

Answer 2

我会这样解决：

ncols=d;
colsum=p;
base=(0:colsum)';
v=@(dm)permute(base,[dm:-1:1]);

M=bsxfun(@plus,base,v(2));
for idx=3:ncols
    M=bsxfun(@plus,M,v(idx));
end
L=M<=colsum;
A=cell(1,ncols);
[A{:}]=ind2sub(size(L),find(L));
a=cell2mat(A);
%subtract 1 because 1 based indexing but base starts at 0
a=a-1+min(base);

它建立了一个包含总和的 p 维矩阵。这段代码的效率取决于sum(L(:))/numel(L)，这个商告诉你有多少创建的矩阵实际用于解决方案。如果这对您的输入来说太低了，那么可能会有更好的解决方案。

Answer 3

这是一个非常简单的方法：

L = dec2base(0:4^3-1,4);
idx=sum(num2str(L)-'0',2)<=3;
L(idx,:)

我认为第一行可以非常有效地创建候选人列表，但不幸的是我不知道如何以有效的方式减少列表。

所以第二行有效，但可以明智地使用改进性能。