正如您所指出的，您的问题有两个方面。一方面，你有记忆问题，因为你需要做很多试验。另一方面，您有性能问题，因为您必须处理所有这些试验。
每个问题的解决方案通常会对另一个问题产生负面影响。恕我直言，最好的方法是找到妥协。
只有在您摆脱矢量化所需的那些 庞大 数组并使用不同的策略来执行循环时，才有可能进行更多试验。我将优先考虑使用更多试验的可能性，可能以最佳性能为代价。
当我在Matlab profiler 中按原样执行您的代码时，它立即显示所有变量的初始内存分配需要大量时间。它还表明plot 和hold all 命令是它们中最耗时的行。更多的反复试验表明，在OUT OF MEMORY 错误开始出现之前，您可以执行的trials 的最大值令人失望地低。
如果您了解一些有关其在 Matlab 中的局限性的信息，则可以极大地加速循环。在旧版本的 Matlab 中，应该完全避免循环以支持“矢量化”代码，这是事实。在最近的版本中（我相信是 R2008a 及更高版本），Mathworks 引入了一种称为 JIT 加速器（即时编译器）的技术，它在执行期间将 M 代码即时翻译成 机器 语言.简而言之，JIT 加速器允许您的代码绕过 Matlab 的解释器并更直接地与底层硬件对话，这可以节省大量时间。
你会听到很多关于在 Matlab 中应该避免循环的建议，现在不再是一般正确的了。尽管矢量化仍然有其价值，但仅使用矢量化代码实现的任何相当复杂的过程通常都难以辨认、难以理解、难以更改和难以维护。使用循环的相同过程的实现通常没有这些缺点，而且它通常会更快并且需要更少的内存。
不幸的是，JIT 加速器有一些令人讨厌的（恕我直言，不必要的）限制，您必须了解这些限制。
这样的事情之一是plot；通常最好让循环什么都不做，除了收集和操作数据，并将任何绘图命令等延迟到循环之后。
另一个这样的事情是hold； hold 函数不是 Matlab 内置函数，也就是说，它是用 M 语言实现的。 Matlab 的 JIT 加速器在循环中使用时无法加速非内置函数，这意味着您的整个循环将以 Matlab 的解释速度运行，而不是机器语言的速度！因此，还要将此命令延迟到 循环之后 :)
现在，如果您想知道，最后一步可能会产生巨大的不同——我知道有一种情况是将函数体复制粘贴到上层循环带来了 1200 倍的性能提升。 天的执行时间已减少到分钟！）。
实际上，您的循环中还有另一个次要 问题（它非常小，而且相当不方便，我会立即同意）—— 循环变量的名称不应该是i。 i这个名字是Matlab中虚数单位的名字，名字解析也会在每次迭代中不必要地消耗时间。它很小，但不容忽视。
现在，考虑到这一切，我得出了以下实现：
function [hot, cold, h] = MuonTracks(tracks)
% NOTE: no variables larger than 1x1 are initialized
width = 1e-4;
height = 2e-4;
% constant used for Laplacian noise distribution
bL = 15 / sqrt(2);
% Loop through all tracks
X = [];
hot = 0;
ii = 0;
while ii <= tracks
ii = ii + 1;
% Note that I've inlined (== copy-pasted) the original laprnd()
% function call. This was necessary to work around limitations
% in loops in Matlab, and prevent the nececessity of those HUGE
% variables.
%
% Of course, you can still easily generalize all of this:
% the new data
u = rand-0.5;
Ystart = 15;
Xstart = 800*rand-400;
Xend = Xstart - bL*sign(u)*log(1-2*abs(u));
b = (Ystart*Xend)/(Xend-Xstart);
% the test
if ((b < height && b > 0)) ||...
(Xend < width/2 && Xend > -width/2)
hot = hot+1;
% growing an array is perfectly fine when the chances of it
% happening are so slim
X = [X [Xstart; Xend]]; %#ok
end
end
% This is trivial to do here, and prevents an 'else' in the loop
cold = tracks - hot;
% Now plot the chosen ones
h = figure;
hold all
Y = repmat([15;0], 1, size(X,2));
plot(X, Y, 'r');
end
有了这个实现，我可以这样做：
>> tic, MuonTracks(1e8); toc
Elapsed time is 24.738725 seconds.
内存占用完全可以忽略不计。
分析器现在还显示了代码中工作量的均匀分布；没有因为内存使用或性能而真正脱颖而出的行。
这可能不是最快的实现（如果有人看到明显的改进，请随时编辑它们）。但是，如果你愿意等待，你将能够做到MuonTracks(1e23)（或更高：）
我还用C做了一个实现，可以编译成Matlab MEX文件：
/* DoMuonCounting.c */
#include <math.h>
#include <matrix.h>
#include <mex.h>
#include <time.h>
#include <stdlib.h>
void CountMuons(
unsigned long long tracks,
unsigned long long *hot, unsigned long long *cold, double *Xout);
/* simple little helper functions */
double sign(double x) { return (x>0)-(x<0); }
double rand_double() { return (double)rand()/(double)RAND_MAX; }
/* the gateway function */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int
dims[] = {1,1};
const mxArray
/* Output arguments */
*hot_out = plhs[0] = mxCreateNumericArray(2,dims, mxUINT64_CLASS,0),
*cold_out = plhs[1] = mxCreateNumericArray(2,dims, mxUINT64_CLASS,0),
*X_out = plhs[2] = mxCreateDoubleMatrix(2,10000, mxREAL);
const unsigned long long
tracks = (const unsigned long long)mxGetPr(prhs[0])[0];
unsigned long long
*hot = (unsigned long long*)mxGetPr(hot_out),
*cold = (unsigned long long*)mxGetPr(cold_out);
double
*Xout = mxGetPr(X_out);
/* call the actual function, and return */
CountMuons(tracks, hot,cold, Xout);
}
// The actual muon counting
void CountMuons(
unsigned long long tracks,
unsigned long long *hot, unsigned long long *cold, double *Xout)
{
const double
width = 1.0e-4,
height = 2.0e-4,
bL = 15.0/sqrt(2.0),
Ystart = 15.0;
double
Xstart,
Xend,
u,
b;
unsigned long long
i = 0ul;
*hot = 0ul;
*cold = tracks;
/* seed the RNG */
srand((unsigned)time(NULL));
/* aaaand start! */
while (i++ < tracks)
{
u = rand_double() - 0.5;
Xstart = 800.0*rand_double() - 400.0;
Xend = Xstart - bL*sign(u)*log(1.0-2.0*fabs(u));
b = (Ystart*Xend)/(Xend-Xstart);
if ((b < height && b > 0.0) || (Xend < width/2.0 && Xend > -width/2.0))
{
Xout[0 + *hot*2] = Xstart;
Xout[1 + *hot*2] = Xend;
++(*hot);
--(*cold);
}
}
}
在 Matlab 中编译
mex DoMuonCounting.c
(在运行mex setup :) 之后 然后将它与像这样的小型 M-wrapper 一起使用：
function [hot,cold, h] = MuonTrack2(tracks)
% call the MEX function
[hot,cold, Xtmp] = DoMuonCounting(tracks);
% process outputs, and generate plots
hot = uint32(hot); % circumvents limitations in 32-bit matlab
X = Xtmp(:,1:hot);
clear Xtmp
h = NaN;
if ~isempty(X)
h = figure;
hold all
Y = repmat([15;0], 1, hot);
plot(X, Y, 'r');
end
end
这让我可以做
>> tic, MuonTrack2(1e8); toc
Elapsed time is 14.496355 seconds.
请注意，MEX 版本的内存占用略大，但我认为这没什么好担心的。
我看到的唯一缺陷是 Muon 计数的最大固定数（硬编码为 10000 作为 Xout 的初始数组大小；需要，因为标准 C 中没有动态增长的数组）......如果你'再次担心这个限制可能会被打破，只需增加它，将其更改为等于tracks 的一小部分，或者做一些更聪明（但更痛苦）的动态数组增长技巧。