从 3 个向量及其原始坐标生成所有可能的三元组对组合
Generate all possible combinations of pairs of triplets from 3 vectors together with their original coordinates
给定 Matlab 中的 3 个行向量 A,B,C
,我想生成矩阵 D
报告来自 A,B,C
的三元组对的所有可能组合及其在 A,B,C
。
我已经编写了一个代码来完成我想要的。由于我正在尝试尽可能多地优化我的代码(代码应该重复数百万次),我想知道您是否可以想出更有效的解决方案。例如,在我的代码中,我没有预先分配矩阵 D
,因为我不知道如何获取每对三元组的索引,这样效率不高。
下面的代码会更好地解释它:
clear
A=[1 2];
B=[-4 -2 5];
C=[8 9 -3 0];
sA=size(A,2);
sB=size(B,2);
sC=size(C,2);
sT=sA*sB*sC;
%Generate the matrix D of dimension [sT*(sT-1)/2]x[12]
%reporting all the possible combinations of pairs of triplets from A,B,C
%together with their original coordinates in A,B,C
[ca, cb, cc] = ndgrid(A, B, C);
T = [ca(:), cb(:), cc(:)]; %matrix of dimension sTx3 reporting all the possible triplets
%from A,B,C
[ca, cb, cc] = ndgrid(1:sA, 1:sB, 1:sC);
Tcoord = [ca(:), cb(:), cc(:)]; %matrix of dimension sTx3 reporting the coordinates of all
%the possible triplets from A,B,C
D=[];
for w=1:sA*sB*sC
for r=w+1:sA*sB*sC
D=[D; T(w,:) T(r,:) Tcoord(w,:) Tcoord(r,:)];
end
end
填充矩阵 D
的最后一个嵌套 for 循环效率更高。 OP 在他们的声明中是正确的:
"For example, in my code I am not pre-allocating the matrix D because I don't know how to get the index of each pair of triplets and this is not efficient."
我们可以对这些循环中的大部分工作进行向量化,方法是注意到 OP 在他们关于矩阵最终大小的评论中提到了一种模式 D
(即 Generate the matrix D of dimension [sT*(sT-1)/2]x[12]
) .任何熟悉系列和序列的人都会熟悉第一个维度。就是Triangle Numbers.
的公式
考虑到这一点,我们可以看到在最终结果中,前 3 列(以及第 7 列到第 9 列)重复了 23 次,然后是 22 次,依此类推,而第 4 列到第 6 列(以及第 10 列)到 12) 是 T/Tcoord
的最后 23 行,T/Tcoord
的最后 22 行等等。在代码中我们有:
D1 = zeros(sT * (sT - 1) / 2, 12);
s = 1;
e = sT - 1;
for w = 1:(sT - 1)
D1(s:e,[1:3,7:9]) = repmat([T(w,:),Tcoord(w,:)], sT - w, 1);
D1(s:e,[4:6,10:12]) = [T((w+1):sT,:),Tcoord((w+1):sT,:)];
s = e + 1;
e = e + (sT - (w + 1));
end
和运行每个方法使用tic
和toc
200次,我们发现我们的效率几乎提高了35%
。
% OP's setup code goes here
tic
for i=1:200
D=[];
for w=1:sA*sB*sC
for r=w+1:sA*sB*sC
D=[D; T(w,:) T(r,:) Tcoord(w,:) Tcoord(r,:)];
end
end
end
toc
tic
for i = 1:200
D1 = zeros(sT * (sT - 1) / 2, 12);
s = 1;
e = sT - 1;
for w = 1:(sT - 1)
D1(s:e,[1:3,7:9]) = repmat([T(w,:),Tcoord(w,:)], sT - w, 1);
D1(s:e,[4:6,10:12]) = [T((w+1):sT,:),Tcoord((w+1):sT,:)];
s = e + 1;
e = e + (sT - (w + 1));
end
end
toc
% Gives same result
isequal(D, D1)
% Timing for 200 runs on 24 total combinations
Elapsed time is 2.09613 seconds.
Elapsed time is 1.35988 seconds.
ans = 1
如果我们使输入向量更大,我们会看到效率有更大的提高。以下是 运行 在以下向量上对每种方法进行 15 次的结果:
A=[1 2 3 4 23];
B=[-4 -2 5 74];
C=[8 9 -3 0];
% Timing for 15 run on 80 total combinations
Elapsed time is 4.00448 seconds.
Elapsed time is 0.379919 seconds.
ans = 1
快了 10 倍多。随着输入向量的大小变大,差距呈指数增长。
A=[1 2 3 4 23];
B=[-4 -2 5 74 28];
C=[8 9 -3 0 -100 -5];
% Timing for 1 run on 150 total combinations
Elapsed time is 3.63065 seconds.
Elapsed time is 0.0481789 seconds.
ans = 1
大约快了 75 倍!!!
更新
OP 在评论中给出了更好的答案:
indices=nchoosek((1:1:sT),2);
D=[T(indices(:,1),:) T(indices(:,2),:) Tcoord(indices(:,1),:) Tcoord(indices(:,2),:)];
这是我用来做基准测试的代码:
clear
A=[1 2 3 4 23 24 25 26];
B=[-4 -2 5 74 28 10 11 12 13];
C=[8 9 -3 0 -100 -5 60 120];
sA=size(A,2);
sB=size(B,2);
sC=size(C,2);
sT=sA*sB*sC;
tic
for i = 1:10
[ca, cb, cc] = ndgrid(A, B, C);
T = [ca(:), cb(:), cc(:)];
[ca, cb, cc] = ndgrid(1:sA, 1:sB, 1:sC);
Tcoord = [ca(:), cb(:), cc(:)];
D1 = zeros(sT * (sT - 1) / 2, 12);
s = 1;
e = sT - 1;
for w = 1:(sT - 1)
D1(s:e,[1:3,7:9]) = repmat([T(w,:),Tcoord(w,:)], sT - w, 1);
D1(s:e,[4:6,10:12]) = [T((w+1):sT,:),Tcoord((w+1):sT,:)];
s = e + 1;
e = e + (sT - (w + 1));
end
end
toc
tic
for i = 1:10
indices=nchoosek((1:1:sT),2);
D=[T(indices(:,1),:) T(indices(:,2),:) Tcoord(indices(:,1),:) Tcoord(indices(:,2),:)];
end
toc
isequal(D, D1)
结果如下:
% Timing for 10 runs on 576 total combinations
Elapsed time is 1.9834 seconds.
Elapsed time is 0.13818 seconds.
ans = 1
我提供的改进解决方案比原来的解决方案好很多,但与 OP 更新的解决方案不匹配。它的订单速度更快,非常优雅,我可能会补充。
给定 Matlab 中的 3 个行向量 A,B,C
,我想生成矩阵 D
报告来自 A,B,C
的三元组对的所有可能组合及其在 A,B,C
。
我已经编写了一个代码来完成我想要的。由于我正在尝试尽可能多地优化我的代码(代码应该重复数百万次),我想知道您是否可以想出更有效的解决方案。例如,在我的代码中,我没有预先分配矩阵 D
,因为我不知道如何获取每对三元组的索引,这样效率不高。
下面的代码会更好地解释它:
clear
A=[1 2];
B=[-4 -2 5];
C=[8 9 -3 0];
sA=size(A,2);
sB=size(B,2);
sC=size(C,2);
sT=sA*sB*sC;
%Generate the matrix D of dimension [sT*(sT-1)/2]x[12]
%reporting all the possible combinations of pairs of triplets from A,B,C
%together with their original coordinates in A,B,C
[ca, cb, cc] = ndgrid(A, B, C);
T = [ca(:), cb(:), cc(:)]; %matrix of dimension sTx3 reporting all the possible triplets
%from A,B,C
[ca, cb, cc] = ndgrid(1:sA, 1:sB, 1:sC);
Tcoord = [ca(:), cb(:), cc(:)]; %matrix of dimension sTx3 reporting the coordinates of all
%the possible triplets from A,B,C
D=[];
for w=1:sA*sB*sC
for r=w+1:sA*sB*sC
D=[D; T(w,:) T(r,:) Tcoord(w,:) Tcoord(r,:)];
end
end
填充矩阵 D
的最后一个嵌套 for 循环效率更高。 OP 在他们的声明中是正确的:
"For example, in my code I am not pre-allocating the matrix D because I don't know how to get the index of each pair of triplets and this is not efficient."
我们可以对这些循环中的大部分工作进行向量化,方法是注意到 OP 在他们关于矩阵最终大小的评论中提到了一种模式 D
(即 Generate the matrix D of dimension [sT*(sT-1)/2]x[12]
) .任何熟悉系列和序列的人都会熟悉第一个维度。就是Triangle Numbers.
考虑到这一点,我们可以看到在最终结果中,前 3 列(以及第 7 列到第 9 列)重复了 23 次,然后是 22 次,依此类推,而第 4 列到第 6 列(以及第 10 列)到 12) 是 T/Tcoord
的最后 23 行,T/Tcoord
的最后 22 行等等。在代码中我们有:
D1 = zeros(sT * (sT - 1) / 2, 12);
s = 1;
e = sT - 1;
for w = 1:(sT - 1)
D1(s:e,[1:3,7:9]) = repmat([T(w,:),Tcoord(w,:)], sT - w, 1);
D1(s:e,[4:6,10:12]) = [T((w+1):sT,:),Tcoord((w+1):sT,:)];
s = e + 1;
e = e + (sT - (w + 1));
end
和运行每个方法使用tic
和toc
200次,我们发现我们的效率几乎提高了35%
。
% OP's setup code goes here
tic
for i=1:200
D=[];
for w=1:sA*sB*sC
for r=w+1:sA*sB*sC
D=[D; T(w,:) T(r,:) Tcoord(w,:) Tcoord(r,:)];
end
end
end
toc
tic
for i = 1:200
D1 = zeros(sT * (sT - 1) / 2, 12);
s = 1;
e = sT - 1;
for w = 1:(sT - 1)
D1(s:e,[1:3,7:9]) = repmat([T(w,:),Tcoord(w,:)], sT - w, 1);
D1(s:e,[4:6,10:12]) = [T((w+1):sT,:),Tcoord((w+1):sT,:)];
s = e + 1;
e = e + (sT - (w + 1));
end
end
toc
% Gives same result
isequal(D, D1)
% Timing for 200 runs on 24 total combinations
Elapsed time is 2.09613 seconds.
Elapsed time is 1.35988 seconds.
ans = 1
如果我们使输入向量更大,我们会看到效率有更大的提高。以下是 运行 在以下向量上对每种方法进行 15 次的结果:
A=[1 2 3 4 23];
B=[-4 -2 5 74];
C=[8 9 -3 0];
% Timing for 15 run on 80 total combinations
Elapsed time is 4.00448 seconds.
Elapsed time is 0.379919 seconds.
ans = 1
快了 10 倍多。随着输入向量的大小变大,差距呈指数增长。
A=[1 2 3 4 23];
B=[-4 -2 5 74 28];
C=[8 9 -3 0 -100 -5];
% Timing for 1 run on 150 total combinations
Elapsed time is 3.63065 seconds.
Elapsed time is 0.0481789 seconds.
ans = 1
大约快了 75 倍!!!
更新
OP 在评论中给出了更好的答案:
indices=nchoosek((1:1:sT),2);
D=[T(indices(:,1),:) T(indices(:,2),:) Tcoord(indices(:,1),:) Tcoord(indices(:,2),:)];
这是我用来做基准测试的代码:
clear
A=[1 2 3 4 23 24 25 26];
B=[-4 -2 5 74 28 10 11 12 13];
C=[8 9 -3 0 -100 -5 60 120];
sA=size(A,2);
sB=size(B,2);
sC=size(C,2);
sT=sA*sB*sC;
tic
for i = 1:10
[ca, cb, cc] = ndgrid(A, B, C);
T = [ca(:), cb(:), cc(:)];
[ca, cb, cc] = ndgrid(1:sA, 1:sB, 1:sC);
Tcoord = [ca(:), cb(:), cc(:)];
D1 = zeros(sT * (sT - 1) / 2, 12);
s = 1;
e = sT - 1;
for w = 1:(sT - 1)
D1(s:e,[1:3,7:9]) = repmat([T(w,:),Tcoord(w,:)], sT - w, 1);
D1(s:e,[4:6,10:12]) = [T((w+1):sT,:),Tcoord((w+1):sT,:)];
s = e + 1;
e = e + (sT - (w + 1));
end
end
toc
tic
for i = 1:10
indices=nchoosek((1:1:sT),2);
D=[T(indices(:,1),:) T(indices(:,2),:) Tcoord(indices(:,1),:) Tcoord(indices(:,2),:)];
end
toc
isequal(D, D1)
结果如下:
% Timing for 10 runs on 576 total combinations
Elapsed time is 1.9834 seconds.
Elapsed time is 0.13818 seconds.
ans = 1
我提供的改进解决方案比原来的解决方案好很多,但与 OP 更新的解决方案不匹配。它的订单速度更快,非常优雅,我可能会补充。