MATLAB：高效生成多索引的大整数矩阵
MATLAB: efficient generation of a large integer matrix of multi-indices
令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).
如果你对这个奇怪矩阵的起源感到好奇，请告诉我！
我会这样解决：
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)，此商告诉您创建的矩阵中有多少实际用于解决方案。如果这对您的输入来说很低，可能存在更好的解决方案。
解决方案使用 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 解决方案会更快。
旧方法：
这是 Dennis 的 dec2base 方法的替代方法，后者对最大基数有限制。它仍然比上面的解决方案慢很多。
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,:);
这是一个非常简单的方法：
L = dec2base(0:4^3-1,4);
idx=sum(num2str(L)-'0',2)<=3;
L(idx,:)
我认为第一行可以非常节省时间来创建候选人列表，但不幸的是我不知道之后如何有效地减少列表。
所以第二行有效，但可以明智地使用改进性能。
令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).
如果你对这个奇怪矩阵的起源感到好奇，请告诉我！
我会这样解决：
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)，此商告诉您创建的矩阵中有多少实际用于解决方案。如果这对您的输入来说很低，可能存在更好的解决方案。
解决方案使用 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 解决方案会更快。
旧方法：
这是 Dennis 的 dec2base 方法的替代方法，后者对最大基数有限制。它仍然比上面的解决方案慢很多。
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,:);
这是一个非常简单的方法：
L = dec2base(0:4^3-1,4);
idx=sum(num2str(L)-'0',2)<=3;
L(idx,:)
我认为第一行可以非常节省时间来创建候选人列表，但不幸的是我不知道之后如何有效地减少列表。
所以第二行有效，但可以明智地使用改进性能。