我有下面的代码，现在太慢了，我正在寻找更好的方法来编写它。
这个代码是一个更长代码的一部分，包括for循环t=1:t，因此下面的代码对每个t运行。
我有F家公司，他们生产Yd。在每个时期，他们购买一定数量的K机床，其生产率由a给出；因此，在t中购买的机床的生产能力等于AK=A.*K。企业的总生产能力等于总和(AK(。由于总和(AK(可能大于他们想要生产的Yd，公司必须确定他们将利用每台机器的哪个分数(ω(；v〃；他们拥有，考虑到他们首先开始使用生产率最高的机床(A最高的机床(。因此，从最具生产力的机器开始，如果某台机器的生产能力低于Yd，则该机器的ω将被设置为1，如果生产能力已经达到Yd，那么该机器的Ω将被设置成0，否则ω将被设为所需的分数。我该如何编码？
这就是我所做的，但肯定太慢了：
T=100;
F=50;
A=rand(T,F); %initialized randomly (not as in my real code)
K=rand(T,F); %idem
AK=A.*K; %productivity of machine tools
Yd=rand(1,F); %initialized randomly
omega=zeros(T,F); %needs to be computed
OAK(:,:)=zeros(T,F); %needs to be computed
for f=1:F
[A_sorted_value,A_sorted_index]=sort(A(:,f)); %sort machines according to their productivity level
while sum(A_sorted_value)>0 && sum(OAK(1:t,f))<Yd(f) %apply rule as long as there are machines available and desired production has not been achieved
v=A_sorted_index(end); %pick the best machine of capital, i.e. highest productivity A
if AK(v,f)+sum(OAK(1:t,f))<=Yd(f) %assign omega=1 if existing capacity is below desired production
omega(v,f)=1;
OAK(v,f)=omega(v,f).*AK(v,f); %update existing capacity
elseif AK(v,f)+sum(OAK(1:t,f))>Yd(f) && Yd(f)>sum(OAK(1:t,f)) %assign omega so that machine v can fill desired production
omega(v,f)=(Yd(f)-sum(OAK(1:t,f)))/AK(v,f);
OAK(v,f)=omega(v,f).*AK(v,f);
end
A_sorted_index(end)=[]; %remove machine v from the list and pick the next one if condition "while" is satisfied
A_sorted_value(end)=[]; %otherwise go to the next firm
end
end
我希望它是明确的，并提前感谢！
您似乎正在运行某种模拟。这意味着(如果结果取决于以前的迭代(，则无法优化t和f循环。因此，唯一可用的优化目标将是内部while循环。我建议你考虑一下你是否可以预先计算或向量化任何东西(在纸上检查你的数学(。
但既然我真的帮不了你，这里有一些可能会有所帮助的东西。
我可以立即提供的一个低挂结果是，您可以使用sort函数在某个集合方向上对整个矩阵进行排序，这样您就可以将其移动到f循环之外。
% Sort the entire A array, column-wise.
[AVals, AIdxs] = sort(A,1);
for f = 1:F
A_sorted_value = AVal(:,f);
A_sorted_index = AIdx(:,f);
...
end
接下来，我建议您使用for循环，而不是while循环。
由于您已经知道将在整个已排序的A列表中进行迭代，或者直到满足另一个条件(sum(OAK(1:t,f))<Yd(f)(，因此可以使用for循环来完成此操作。结合上面的建议，你可能会得到这样的东西：
% Sort A in descending order instead, so first element in the array is now highest in value
[AVal, AIdx] = sort(A,1,'descend');
for f = 1:F
A_sorted_value = AVal(:,f); % Appears to be unused?
A_sorted_index = AIdx(:,f);
% Iterate through entire A_sorted_index array, eliminates rewriting of A_sorted_index
for i = 1:length(A_sorted_index)
v = A_sorted_index(i);
% Loop Exit Conditions
if (sum(OAK(1:t,f))<Yd(f)) || (sum(A_sorted_value(1:i)) <= 0)
break
end
% Your code below
%assign omega=1 if existing capacity is below desired production
if AK(v,f)+sum(OAK(1:t,f)) <= Yd(f)
omega(v,f)=1;
%update existing capacity
OAK(v,f)=omega(v,f).*AK(v,f);
%assign omega so that machine v can fill desired production
elseif AK(v,f)+sum(OAK(1:t,f))>Yd(f) && Yd(f)>sum(OAK(1:t,f))
omega(v,f) = (Yd(f)-sum(OAK(1:t,f)))/AK(v,f);
OAK(v,f) = omega(v,f).*AK(v,f);
end
end
end
现在，A_sorted_value似乎未使用。我不能假设A中元素的值是非零的正数。所以我必须包含检查条件sum(A_sorted_value(1:i)) <= 0。但是，如果只有非零的正A元素，则可以删除该片段以进一步加快循环速度。
我将此代码与旧代码进行了比较，在T=200，F=100，在循环中运行了50次。
new = 6.835 secs
old = 65.634 secs
这大约减少了90%的时间。
减少的大量时间来自于删除这些代码行，这是相当浪费的。
A_sorted_index(end)=[];
A_sorted_value(end)=[];
我相信还有更多的优化空间，但这对我来说已经足够了。
请验证我提供的代码是否正确。
这是一个速度更快的版本，它存储了多次使用的一些结果。除了momocha的建议之外，您应该能够实施这些更改。
为了大幅提高速度，您应该转置所有矩阵，使f按列而非行访问数组，或者交换F和Tfor循环的顺序。一次访问矩阵的单列比处理单行更高效，所以总是先遍历列，然后遍历该列中的行。
% precompute sorted array
[A_sorted,A_sorted_I]=sort(A);
for f = 1:F
% get sorted vectors from percomputed array
A_sorted_value = A_sorted(:,f);
A_sorted_index = A_sorted_I(:,f);
% use nnz() instead of sum() (this may not make a difference)
while nnz(A_sorted_value)>0 && sum(OAK(1:t,f))<Yd(f)
% precompute a few things that are used several times
v = A_sorted_index(end);
ak_ = AK(v,f);
sOAK = sum(OAK(1:t,f));
if ak_+sOAK<=Yd(f)
omega(v,f)=1;
OAK(v,f)=ak_; % don't multiply by omega(v,f), we know it's one, this saves memory lookup
else % doesn't need to be an else if, condition will always be true
ex = Yd(f)-sOAK;
omega(v,f)=ex/ak_;
OAK(v,f)=ex;
end
% this is a faster way to de-allocate
A_sorted_index=A_sorted_index(1:end-1);
A_sorted_value=A_sorted_value(1:end-1);
end
end