我刚刚在 Octave/Matlab 中创建了一个代码片段，旨在创建一个显示初始概率向量准确性的图从转移概率矩阵导出使用其幂的马尔可夫链 (）。
n = 1:10000;
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^10000,2);
for i=1:n(end)
y(i) = (tProbs^n(i)*[1;1]-iProbs)(1);
end
有没有办法在提到的行上使用矢量化来删除主循环并希望提高运行时间？我无法将矩阵提升到向量的幂.
PD
这并不能解决我的问题。我正在寻找一种矢量化技术。以防万一有任何混淆。
我也希望这是发布此问题的正确 Stack Exchange。如果它完全不合适，如果是 Mod，我将不胜感激。可以将这个问题迁移到更合适的地方。
甚至在讨论任何矢量化之前，您提供的代码都存在严重的低效率。仅通过使用基本技术，我就获得了 21 倍的加速：
n = 1000000;
y = zeros(1,n);
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);
iProbsSqr = iProbs.^2;
V0 = [1;1];
for i=1:n
V0 = tProbs*V0;
y(i) = sqrt(dot(V0.^2,iProbsSqr));
end
如果以前的版本在您的机器上需要 35 秒，那么这应该需要大约 1.5 秒。如果不显着增加内存使用量，就不可能进一步改进。此外，它看起来非常丑陋和神秘（参见 Roedy Green 的 How to Write Unmaintainable Code - 这是一个链接：https ://www.se.rit.edu/~tabeec/RIT_441/Resources_files/How%20To%20Write% 20Unmaintainable%20Code.pdf )。如果您仍想做矢量化，请告诉我；我认为使用repmat一个可以以矢量化方式预先计算所有矩阵幂，然后进行张量缩减。这可能会起作用。
编辑：有一个小的修改可以提高 Octave 的效率（不影响 MATLAB）。Octave 在函数调用和仅调用dot而不是使用方面有巨大的开销'*。sum并且sqrt似乎很好。
n = 1000000;
y = zeros(1,n);
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);
iProbsSqr = iProbs.^2;
V0 = [1;1];
for i=1:n
V0 = tProbs*V0;
y(i) = sqrt((V0.^2)'*iProbsSqr);
end
这段代码的向量化非常困难。它需要计算所有矩阵的幂，n而使用内置函数无法有效地做到这一点（好吧，整个想法只是冗余计算，但可以通过它看到过去）。如果您可以预先计算所有必要的矩阵幂，这将成为矩阵向量和向量向量乘法的练习。
您可以尝试使用 oct 文件预编译 for 循环（请参阅https://octave.org/doc/v6.2.0/Getting-Started-with-Oct_002dFiles.html）
编辑 2：交付的不可维护的代码。
n = 1:1000000;
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);
[V,D] = eig(tProbs);
d = diag(D);
dp = d.^n;
aux1 = inv(V)*[1;1];
aux1 = dp.*aux1;
aux1 = V*aux1;
aux1 = aux1.^2;
aux1 = aux1'*(iProbs.^2);
aux1 = sqrt(aux1);
如果我不解释为什么这段代码无法维护，我会感觉很糟糕。所以我开始了：我们必须同时计算所有矩阵的幂，但这很愚蠢，而且是很多冗余计算。但是，我们在这里没有使用 C。冗余计算是获得加速的好工具。
但是你问，我们如何一次计算所有的矩阵幂？我们假设矩阵是可对角化的（唷，幸运的是这次是。希望下一个人足够聪明来检查），因为如果矩阵是可对角化的，即如果存在对角矩阵和一个可逆矩阵， 然后.
但这不只是在推动问题吗？并不真地！我们从矩阵中取出对角线部分，将其放入一个向量中，我们可以对该向量进行元素幂运算，n以获得一个 2 乘以多少矩阵dp。（耶！代码变得越来越难以遵循！）。怎么办？
只需使用一系列操作计算结果：乘和，将结果与向量dp相乘，将结果与，...希望你跟着。
生成的代码在 Octave 中运行时间为 0.1 秒（接近我使用 MATLAB 更简单、更易于维护的代码所得到的），并且相同的代码在 MATLAB 中运行时间为 0.01 秒。我的偏好是在 MATLAB 上使用更简单的版本，但 MATLAB 是专有的，所以请随意使用复杂/不可维护的版本 :) 不要忘记评论它。至少如果你必须自己回去阅读它，你可以理解它。
编辑 3：修改以适应对问题的编辑。这是代码：
n = 1:1000000;
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);
[V,D] = eig(tProbs);
d = diag(D);
dp = d.^n;
aux1 = inv(V)*[1;1];
aux1 = dp.*aux1;
aux1 = V*aux1;
aux1 = aux1-iProbs;
aux1 = aux1(1,:);
这是此代码不可维护的另一个原因。一个简单的操作更改需要大量的代码重组。