3

如果我有一个颂歌并以两种方式编写它，例如这里：

function re=rabdab()
x=linspace(0,2000,2000)';
tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;



function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

时间上几乎没有区别。一个需要和另一个一样长。但我认为那fun是fun2. 还是我在这里弄错了？目的是加快我的代码速度。该示例取自 matlab 网页。我想我还没有真正理解“矢量化”是什么意思。如果这已经是矢量化的，那么非矢量化代码会是什么样子？

4

2 回答 2

4

矢量化是一个与函数式编程密切相关的概念。在 MATLAB 中，这意味着在不隐式编写循环的情况下对数组（向量或矩阵）执行操作。

例如，如果你要计算1 到 100 之间f(x) = 2x的每个整数的函数，你可以写：x

for x = 1:100
    f(x) = 2 * x;
end

这不是矢量化代码。矢量化版本是：

x = 1:100; %// Declare a vector of integer values from 1 to 100
f = 2 * x; %// Vectorized operation "*"

甚至更短：

f = 2 * (1:100);

MATLAB 是一种解释型语言，因此很明显，解释器会将其翻译成某种“幕后”循环，但它已经过优化，通常比实际解释循环要快得多（阅读此问题以供参考）。嗯，有点——直到最近的 MATLAB 版本都集成了JIT 加速（在这里阅读）。

现在回到您的代码：这里有两个向量化版本的代码：一个垂直连接三个值，另一个直接将这些值分配到列向量中。这基本上是一样的。如果您使用显式的 for 循环来完成，这将不会被“矢量化”。关于“矢量化”循环（即将 for 循环转换为矢量化操作）所获得的实际性能增益，这首先取决于由于 JIT 加速而导致的 for 循环的实际速度。

似乎没有太多工作可以加快您的代码速度。您的功能非常基本，因此归结为ode45您无法修改的内部实现。

如果您有兴趣进一步阅读有关矢量化和编写更快的 MATLAB 代码的一般性知识，这里有一篇有趣的文章：Loren 谈 MATLAB 的艺术：“加速MATLAB 应用程序”。

快乐编码！

于 2013-05-26T22:42:33.417 回答
1

虽然前面的答案在一般意义上是正确的，但在 ODE 的上下文中向量化意味着更具体的东西。简而言之，一个函数f(t,y)是向量化的，当且仅当f(t,[y1 y2 ...])返回[f(t,y1) f(t,y2) ...], 对于y1,y2列向量。根据文档 [1]，“这允许求解器减少计算雅可比矩阵的所有列所需的函数评估次数。”

在 ODE 意义上，函数fun3及fun4以下函数已正确矢量化：

function re=rabdab()
x=linspace(0,20000,20000)';
opts=odeset('Vectorized','on');

tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
toc;

tic;
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
toc;


function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

function dy = fun3(t,y)
dy = zeros(size(y)); % a matrix with arbitrarily many columns, rather than a column vector
dy = [y(2,:) .* y(3,:);...
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

function dy = fun4(t,y)
dy = [y(2,:) .* y(3,:);... % same as fun3()
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

（作为旁注：省略不必要的内存分配，zeros让fun4运行速度略快于fun3.）

  • 问：关于第一个参数的向量化怎么样？
  • 答：对于 ODE 求解器，ODE 函数仅针对第二个参数进行矢量化。然而，边值问题求解器bvp4c确实需要对第一个和第二个参数进行矢量化。[1]

官方文档 [1] 提供了有关 ODE 特定矢量化的更多详细信息（请参阅“雅可比属性描述”部分）。

[1] https://www.mathworks.com/help/releases/R2015b/matlab/ref/odeset.html

于 2021-01-31T06:53:44.513 回答