Matlab中三点相关函数的矢量化计算？

Question

我想通过向量元素的适当直方图（长度 system_size 的 num_samples 个样本）和相应的聚类函数 T2、T3 来计算向量样本的 2 点和 3 点相关函数 R2、R3 .为简单起见，我正在考虑跨统一箱进行直方图。

什么是矢量化和/或加速以下代码的好方法？

n = length(mesh);
R2 = zeros(n, n);
R3 = zeros(n, n, n);
for sample_id=1:num_samples 
    s = samples(:, sample_id);
    d = mesh(2) - mesh(1);
    % Which bin does the ith sample s belong to?
    bins = ceil((s - mesh(1))/d);

    % Compute two-point correlation function
    for i = 1:system_size
        for j = 1:system_size
            if i ~= j
                R2(bins(i), bins(j))=R2(bins(i), bins(j))+1;
            end
        end
    end

    % Compute three-point correlation function
    for i = 1:system_size
        for j = 1:system_size
            if i ~= j
                for k = 1:system_size
                    if k ~= j && k ~= i
                        R3(bins(i), bins(j), bins(k))=R3(bins(i), bins(j), bins(k))+1;
                        T3(x1, x2, x3) = R3(x1,x2,x3)-R1(x1)*R2(x2,x3)-R1(x2)*R2(x1,x3)...
                             -R1(x3)*R2(x1,x2)+2*R1(x1)*R1(x2)*R1(x3);
                    end
                end
            end
        end
    end
end
R2 = R2/sum(R2(:));
R3 = R3/sum(R3(:));

T3 = zeros(n, n, n);
% Compute three-point cluster function
for i = 1:n
    for j = 1:n
        if i ~= j
            for k = 1:n
                if k ~= j && k ~= i
                    T3(x1, x2, x3) = R3(x1,x2,x3)-R1(x1)*R2(x2,x3)-R1(x2)*R2(x1,x3)...
                         -R1(x3)*R2(x1,x2)+2*R1(x1)*R1(x2)*R1(x3);
                end
            end
        end
    end
end

我天真地认为 hist3(bins, bins...) 或 crosstab(bins, bins) 几乎可以做我想做的事，即寻找向量元素的相关出现，但事实并非如此。

---

例子：

如果我在最外层循环中的输入是

s = [1.2 3.1 4.6 4.7 5.1]
mesh = 0:0.5:6

那么量化后的数据应该是

bins = [3 7 10 10 11]

R2 应该是

>> R2

R2 =

     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     0     0     2     1     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0     2     1     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     2     0     0     0     2     0     0     2     2     0
     0     0     1     0     0     0     1     0     0     2     0     0
     0     0     0     0     0     0     0     0     0     0     0     0

Answer 1

'R2andR3`很简单：

R2 = R2 + 1 - diag(ones(size(R2, 1), 1); % you can replace the loop with this
eye3 = zeros(n, n, n);
eye3(linspace(1, numel(eye3), n)) = 1;
R3 = R3 + 1 - eye3; % can move R3 computation outside the loop

对于T3：

temp = repmat(R2, [1 1 n]).*permute(repmat(R1, [n, 1, n]), [1, 3, 2]);
T3 = R3 - temp - permute(temp, [2 3 1]) - permute(temp, [3 1 2]);
temp2 = repmat(R1'*R1, [1 1 n]).*permute(repmat(R1, [n, 1, n]), [1, 3, 2]);
T3 = T3 + temp2;

假设R1 是一个行向量。

您可能需要稍微尝试一下，因为您的代码中仍有一些不清楚的地方，但这应该与您最终需要的非常接近。

编辑澄清后：

对于R2：

ubins = unique(bins);
bincounts = histc(bins, ubins);
for i=1:max(bincounts)
    indices = find(bincounts == i);
    R2(indices, indices) = R2(indices, indices) + i
end

这仅对大型向量和数组有用。实际上，您正在向量化矩阵块的计算，而不是整个矩阵（因为bins 中可能存在重复）。

您可以为R3 编写类似的内容。 T3 应该仍然与我之前的答案相似。