MATLAB：通过重复向量的元素创建一个大矩阵，每列的步幅增加

Question

在 MATLAB 中，我有一个长度为 n 的向量 x，其中 n 通常为 O(10)，我想构建一个大小为 [n^m,m] 的高矩阵 A，其中 m 再次为 0(10 ）。矩阵有一个特殊的形式：如果 n=4 和 m=6，让

x=[x1; x2; x3; x4]

那么A是

   x1 x1 x1 x1 x1 x1
   x2 x1 x1 x1 x1 x1
   x3 x1 x1 x1 x1 x1
   x4 x1 x1 x1 x1 x1
   x1 x2 x1 x1 x1 x1
   x2 x2 x1 x1 x1 x1
   x3 x2 x1 x1 x1 x1
   x4 x2 x1 x1 x1 x1
   x1 x3 x2 x1 x1 x1
   .               .         
   .               . 
   .               .
   x4 x4 x4 x4 x4 x4

在实践中，每一列都是通过重复 x 的元素获得的，每列的步幅增加。我怎样才能做到这一点？如果可能，我更喜欢高效（矢量化）的解决方案，因为如您所见，A 的行数随着 m 呈指数增长。

编辑：哎呀，对不起！我忘了我还需要基于具有相同长度 x 的向量 w 构建一个大小为 [n^m,1] 的向量 V

w=[w1; w2; w3; w4]

V 是

   w1^6
   w2*w1^5
   w3*w1^5
     .
     .
     .
   w4^6
     

希望糟糕的图形足够清晰。无论如何，V 是一个长度为 n^m 的列向量。猜猜我可以从 w 创建一个矩阵 B，就像从 x 创建一个矩阵 A，然后使用 prod(B,2)？

Answer 1

使用 allcomb tool from MATLAB file-exchange 生成索引[1 2 3 4] 的可能组合，然后使用它们来索引x -

v = repmat({1:numel(x)},1,m);
A = x(fliplr(allcomb(v{:})));

另外，您似乎可以使用 -allcomb(v{:},'matlab') 代替 fliplr。

对于问题的已编辑部分，您可以使用它的修改版本-

V = prod(x(allcomb(v{:})),2)

基准测试

请注意，这些是针对此处发布的可运行解决方案。

基准代码

%// Parameters and input x
n = 10; m = 6;num_runs = 20; x =  randi(9,n,1);

disp('-------- With allcomb')
tic
for runs = 1:num_runs
    v = repmat({1:numel(x)},1,m);
    A = x(fliplr(allcomb(v{:})));
end
toc,    clear v A

disp('-------- With bsxfun')
tic
for runs = 1:num_runs
    A = x(floor(mod(bsxfun(@rdivide, (0:n^m-1).', n.^[0:m-1] ), n)+1)); %//'
end
toc,    clear A

disp('-------- With ttable')
tic
for runs = 1:num_runs
    I = ttable(n*ones(1,m));
    A = x(I);
end
toc,    clear I A

disp('-------- With arrayfun')
tic
for runs = 1:num_runs
    A = cell2mat(arrayfun(@(i)...
        (repmat(reshape(repmat(x',n^(i-1),1),[],1),n^(m-i),1)),1:m,'uni',0));
end
toc

结果

-------- With allcomb
Elapsed time is 6.544981 seconds.
-------- With bsxfun
Elapsed time is 11.547062 seconds.
-------- With ttable
Elapsed time is 15.729932 seconds.
-------- With arrayfun
Elapsed time is 4.319048 seconds.

Answer 2

仅基于内置函数的单行代码（即mod 和非常强大的bsxfun）：

result = x(floor(mod(bsxfun(@rdivide, (0:n^m-1).', n.^[0:m-1] ), n)+1));

Answer 3

试试这个：（都是内置函数）

A = cell2mat(arrayfun(@(i)(repmat(reshape(repmat(x',n^(i-1),1),[],1),n^(m-i),1)),1:m,'UniformOutput',0))

解释：

n = 2;
m = 4;
x = (1:n)';
A = [];
for i = 1:m
%// temp1 is (n^(i-1)) x n matrix with each row equal to x' 
    temp1 = repmat(x',n^(i-1),1); 
%// temp2 is (n^(i-1))*n x 1 column vector with corresponding elements of temp1
    temp2 = reshape(temp1,[],1);
%// temp3 is a (n^(m-i))*(n^(i-1))*n x 1, i.e n^m x 1 column vector with elements of temp2 repeated n^(m-i) times
    temp3 = repmat(temp2,n^(m-i),1);
%// A is appending temp3 into its ith column
    A = cat(2,A,temp3);
end

对于编辑部分：

你可以按照你说的做，即prod(B,2)，其中 B 是使用上述代码计算的矩阵

Answer 4

我认为文件交换中的generalized truth table 功能会对您有所帮助

尝试（未测试）：

I = ttable(n*ones(1,m));
x(I);

Answer 5

好的，对此有一些棘手的单行解决方案。但是，如果顺序无关紧要，则有一个关于 matlab 文件交换的文件提供您正在查看的解决方案。见这里combn。它基本上使用了与其他答案提出的单线相同的技巧。 如果顺序很重要，您可能需要在后面进行一些排序或直接调整源代码。