Matlab - 创建两个数组，考虑两个对象的所有可能值答案

【问题标题】：Matlab - Creating two arrays that consider all possible values for two objectsMatlab - 创建两个数组，考虑两个对象的所有可能值
【发布时间】：2017-06-06 19:56:14
【问题描述】：

我有两个物理“对象”。我用两个不同的数组代表他们的位置。

• 对象 1 仅在 xy 平面中移动

• 对象 2 在所有三个物理维度中移动

目标：矢量化四个for 循环而不扭曲数据。此外，目的是对要与对象 2 比较的对象 1 的所有可能值执行此操作。

这里是sn-p的代码：

Npos = 21;
Nsam = 200;

% dummy initialisation    
AX = rand(1, Npos);
AY = zeros(1, Npos);
AZ = rand(1, Npos);
Bx = rand(Nsam);
By = rand(Nsam);
Bz = rand(Nsam);

for qx = 1 : Npos
    for yx = 1 : Npos
        for zx = 1 : Nsam
            for cx = 1 : Nsam
                Tx2Array( qx, yx, zx, cx ) = sqrt( ( AX( qx ) - Bx( zx, cx ) ).^2 + ( AY( yx ) - By( zx, cx ) ).^2 + ( AZ( yx ) - Bz( zx, cx ) ).^2 );
            end
        end
    end
end
% Result is a 21 x 21 x 200 x 200 matrix filled with all real numbers

传奇

AX、AY、AZ 是 1 x 21 的数组，代表对象 1 的 (x,y=0,z)

AY 全部为零，但为了便于阅读，仍然包含在内（因此没有第五个循环！）

Bx、By、Bz 都是 200 x 200 的数组，代表对象 2 的 (x,y,z)

Npos = 21; Nsam = 200;

上面用到的公式是：

sqrt( (a1-b1)^2 + (a2-b2)^2 + (a3-b3)^2 )

【问题讨论】：

提供minimal reproducible example 可能很有用，即初始化所有变量（使用随机值）
您要做的第一件事是预先分配Tx2Array。即在循环之前写：Tx2Array = zeros(Npos,Npos,Nsam,Nsam)。
可以，但是向量运算不需要预先分配。
我不了解您输入的尺寸，但在我看来，您可以使用 ndgrid 和 pdist2 完成所有操作。
@beaker，我以前从未使用过这些功能。你能创建一个答案来扩展它们吗？

标签： matlab matrix matrix-multiplication

【解决方案1】：

如果您有可用的统计工具箱，您可以使用pdist2 计算对象 1 的每个坐标与对象 2 的每个坐标之间的距离：

[X1, Z1] = ndgrid(AX(:), AZ(:));   % X1 and Z1 will be 21x21
D = pdist2([X1(:), zeros(size(X1(:))), Z1(:)], [Bx(:), By(:), Bz(:)]);

在这种情况下，输出将是一个 441 x 40,000 数组，其中 D(i, j) 为您提供对象 1 的点 i 和对象 2 的点 j 之间的距离，两者都使用线性索引。

【讨论】：

是不是假设计算的数量类似于 21 x 21 x 200 x 200 = 17,640,000 个元素的总数？在这种情况下，只有 21 x 40,000 = 840,000 个元素。
对象 1 只有 21 分，对吧？对象 2 有 40,000 分？这就是我对你的数学感到困惑的原因，因为在我看来，你要么重复计算对象 1 的分数，要么只计算对象 2 的三分之一。
我同意烧杯的观点，这种计算似乎更合乎逻辑。如果您有 21 个 A 点和 40 000 个 B 点，则可以计算 21 x 40 000 个距离。您正在计算 A 的所有可能排列的距离，这可能是您想要的，但通常并非如此。
对象 1 有 21 个“x 坐标”和 21 个“z 坐标”。如上所述，它仅存在于 xz 平面中。这意味着对象有时会沿 x 轴存在，但也可能沿 z 轴存在并介于两者之间。
好吧，你为什么不使用你的代码计算Tx2Array，然后用这个和reshape(D, 21, 21, 200, 200)比较它们呢？

【解决方案2】：

您可以通过将zx 和cx 替换为: 来避免内部循环，如下所示：

Tx2Array = zeros(Npos, Npos, Nsam, Nsam); % preallocate memory
for qx = 1 : Npos
    for yx = 1 : Npos
        Tx2Array( qx, yx, :, : ) = sqrt( ( AX( qx ) - Bx( :, : ) ).^2 + ( AY( yx ) - By( :, : ) ).^2 + ( AZ( yx ) - Bz( :, : ) ).^2 );
    end
end

通过这种方式，最大的维度被矢量化。所以，最大的改进已经完成。

通过将您的 B* 转换为 4D 并为您的 A* 矩阵生成网格，您甚至可以删除所有 for 循环，如下所示：

[AX_, AZ_] = meshgrid(AX, AZ);
AX_ = AX_';
AZ_ = AZ_';
AY_ = zeros(Npos);

Bx_(1, 1, :, :) = Bx;
By_(1, 1, :, :) = By;
Bz_(1, 1, :, :) = Bz;

Tx2Array2 = sqrt( ( AX_ - Bx_ ).^2 + ( AY_ - By_ ).^2 + ( AZ_ - Bz_ ).^2 );

您可以使用以下检查检查结果是否相同：

max(max(max(max(abs(Tx2Array - Tx2Array2))))) < eps

【讨论】：

您可以从Bs 中删除(:,:)
我知道，我只是想明确说明我所做的事情，即将zx 和cx 替换为:
谢谢！这确实将其减少到两个循环。
我认为您可能会在最后一行进行浮点比较时遇到麻烦。您可能希望将其更改为 < eps。
@beaker 这不是问题（在我的电脑上），但我更改了它，因为在比较浮点数时考虑浮点精度总是一个好习惯

【解决方案3】：

如果数组正确初始化，您的任务将非常简单：

用正确的维度初始化数组

AX = rand( Npos,1);
AY = zeros(1, Npos);
AZ = rand(1, Npos);
Bx = rand(1,1,Nsam,Nsam);
By = rand(1,1,Nsam,Nsam);
Bz = rand(1,1,Nsam,Nsam);

然后在 MATLAB r2016b / Octave 中您可以简单地编写：

Tx2Array = sqrt( ( AX - Bx ).^2 + ( AY - By ).^2 + ( AZ - Bz ).^2 );

在 r2016b 之前，您可以使用 bsxfun:

Tx2Array = sqrt(bsxfun(@plus,bsxfun(@plus,bsxfun(@minus,AX , Bx).^2,bsxfun(@minus,AY , By).^2),bsxfun(@minus,AZ , Bz).^2));

【讨论】：