MATLAB：为蒙特卡洛模拟索引一个大矩阵答案

【问题标题】：MATLAB: Indexing a large matrix for Monte Carlo SimulationMATLAB：为蒙特卡洛模拟索引一个大矩阵
【发布时间】：2012-12-15 17:50:54
【问题描述】：

我试图在 MATLAB 中索引一个大型矩阵，该矩阵包含跨行和跨列单调递增的数字，即，如果矩阵被称为A，则对于每个(i,j)、A(i+1,j) > A(i,j) 和A(i,j+1) > A(i,j)。

我需要创建一个随机数 n 并将其与矩阵 A 的值进行比较，以查看该随机数应放置在矩阵 A 中的哪个位置。换句话说，n 的值可能不等于任何矩阵的内容，但它可能位于任意两行和任意两列之间，并确定了一个“bin”，该“bin”标识了它在 A 中的位置。一旦找到这个位置，我就会在新矩阵中增加相应的索引和 A 一样大。

问题是我想这样做 1,000,000 次。我需要创建一个随机数一百万次，并对这些数字中的每一个进行索引检查。这是从一个点降落在屏幕上的一百万个光子的蒙特卡罗模拟；矩阵A由球坐标中的角度组成，随机数为每个入射光子的立体角。

到目前为止，我的代码是这样的（我没有在这里复制粘贴，因为细节并不重要）：

for k = 1:1000000  
    n = rand(1,1)*pi;  
    for i = length(A(:,1))  
        for j = length(A(1,:))  
            if (n > A(i-1,j)) && (n < A(i+1,j)) && (n > A(i,j-1)) && (n < A(i,j+1))  

                new_img(i,j) = new_img(i,j) + 1;   % new_img defined previously as zeros

            end
        end
    end
end

“if”语句只是检查 A 中构成 n 边界的索引。

这工作得很好，但它需要很长的时间，特别是因为我的矩阵 A 是尺寸为 11856 x 11000 的图像。有没有更快/更聪明/更简单的方法？

提前致谢。

【问题讨论】：

代码原样将出错，索引越界到A。
你是对的 - 谢谢。我的实际代码考虑到了这一点；我只是想了解它的要点。
另外，你不应该用“点击”的数量来衡量增量吗？否则，角度越大，光子就越多。
如果您发布surf 渲染A 可能是个好主意。
每个“命中”对应于落在该“bin”中的光子，因此靠近中心的“bin”将有更多“命中”，从而导致中心的传播最强烈。

标签： matlab indexing

【解决方案1】：

您可以通过同时对A 的所有元素执行计算来摆脱内部循环。此外，您可以一次创建所有随机数，而不是一次创建一个。请注意，new_img 的最外层像素永远不能不为零。

randomNumbers = rand(1,1000000)*pi;
new_img = zeros(size(A));
tmp_img = zeros(size(A)-2);

for r = randomNumbers
    tmp_img = tmp_img + A(:,1:end-2)<r & A(:,3:end)>r & A(1:end-1,:)<r & A(3:end,:)>r;
end

new_img(2:end-1,2:end-1) = tmp_img;

/aside：如果数组更小，我会使用bsxfun 进行比较，但使用 OP 中的数组大小，该方法会耗尽内存。

【讨论】：

谢谢。有没有办法在避免通过随机数进行迭代的同时做到这一点？这就是一直在消耗的东西。
不遍历随机数将需要大量内存 - 基本上你会创建一个 imageSize-by-#randomNumbers 数组。我认为在你的情况下这是不可行的。
@shash：但是，如果您能告诉我们更多关于A 的结构的信息，那么可以有更聪明的解决方案。你如何创建A？
实际上有两组“A”矩阵，一组由“theta”角组成；即每列与矩阵中心点之间的角度，距矩阵一定垂直距离。另一个“A”矩阵由相应的行角组成。这很难解释，但它只是“胶片”（“A”矩阵）上每个像素相对于距胶片一定距离的点源的球坐标表示。为简单起见，我将 theta 和 phi 矩阵分开；处理立体角很痛苦。
所以实际上应该有两组随机数：一组用于thetas，一组用于phis，从而定义了每个光子采取的随机方向。

【解决方案2】：

A bin 中的值是边缘吗？即A 是否指定了一个网格？如果是这种情况，那么您可以使用 hist3 快速填充 A。

这是一个例子： numRand = 1e n = 兰迪（100,1e6,1）； nMatrix = [floor(data./10), mod(data,10)];

edges = {0:1:9, 0:10:99};

A = hist3(dataMat, edges);

如果您的A 没有指定网格，那么您应该创建一次所有随机值并对其进行排序。然后遍历这些值。

因为您知道n(i) >= n(i-1)，所以您不必检查对于n(i-1) 来说太小的垃圾箱。这是优化大部分冗余检查的一种非常简单的方法。

【讨论】：

谢谢，这也很有用。

【解决方案3】：

这是一个对内部循环有很大帮助的 sn-p，它会找到小于您的 value 的最大点的位置。

idx1 = A<value
idx2 = A(idx1) == max(A(idx1))

如果你想找到确切的位置，你可以用find 包裹它。

【讨论】：

@Jonas 我不明白为什么它不适用于 2D。如果您在 idx2 中找到一个点，则新值可以在它的正下方或下方。