从 3 个向量及其原始坐标生成三元组对的所有可能组合

Question

给定 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

Answer 1

填充矩阵D 的最后一个嵌套for 循环可以更有效。 OP 在他们的声明中是正确的：

“例如，在我的代码中，我没有预先分配矩阵 D，因为我不知道如何获取每对三元组的索引，这样效率不高。”

我们可以通过注意到 OP 在他们的 cmets 中提到的关于矩阵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 在 cmets 中给出了非常出色的答案：

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 的更新解决方案不匹配。我可能会补充说，它的订单速度更快且相当优雅。