在矩阵的每一行上应用函数的最快方法是什么？

Question

到目前为止我已经搜索过，我知道有几种方法（1，2，3，4）到目前为止我使用了以下代码：

Fv_calc(:,2) = arrayfun(@(n) MaxPositiveRoot2DegreePolynomial(QuadraticCoefficients(n,:)), 1:size(QuadraticCoefficients,1));  

其中 MaxPositiveRoot2DegreePolynomial 是以下函数：

function root = MaxPositiveRoot2DegreePolynomial(c)
    d = c./c(1);
    root = eig([0 -d(3);1 -d(2)]);
    condition = root == abs(root);
    root = root(condition);
    if isempty(root)
        root = 0;
    end
    root = max(root);
end  

其中 QuadraticCoefficients 是一个 62271x3 矩阵，每行包含一般二次方程的系数 a、b、c。 ax^2+bx+c 关于我要解决的问题，所有的根都是真实的，因此我使用了一个修改后的函数来查找根，以免浪费时间检查根是否真实。
通过分析代码，我发现函数 MaxPositiveRoot2DegreePolynomial 的每次运行大约需要 0.000047 秒，对于 62271 次运行大约需要 2.914925371 秒。但是代码

Fv_calc(:,2) = arrayfun(@(n) MaxPositiveRoot2DegreePolynomial(QuadraticCoefficients(n,:)), 1:size(QuadraticCoefficients,1));  

需要 3.198597803 才能运行。所以大约 0.283672431 仅由 arrayfun 占用。我想看看有什么办法可以减少这个时间？

Answer 1

这是一个完美的例子，它显示了你应该在哪里继续使用你的循环。我预期的矢量化解决方案最终比您的循环方法慢。

但这可能取决于您的 Matlab 版本，因为我已经在使用 Matlab R2015b 和 optimzed JIT-Compiler，在旧版本中，矢量化方法可能更快。

[n,~] = size(C);
%// division
d = bsxfun(@rdivide,C,C(:,1));
%// create blocks for eigen matrix block procession
X = ( [zeros(n,1) ones(n,1) -d(:,3) -d(:,2)] ).'; %'
Y = reshape(X,2,[]);
%// use block processing to get eigenvalues of blocks
rt = blockproc(Y,[2 2],@(x) eig(x.data) ).'; %'
%// check if eigen values are real positives
condition = rt == abs(rt);
%// output filter
out = rt(condition);

基准测试

function [t] = bench()
    %// load your data
    DATA = load('matlab');
    C = DATA.QuadraticCoefficients;

    % functions to compare
    fcns = {
        @() thewaywewalk(C);
        @() sepideh(C);
    };

    % timeit
    t = zeros(2,1);
    for ii = 1:10;
        t = t + cellfun(@timeit, fcns);
    end
end

function out = thewaywewalk(C)  %thewaywewalk
    [n,~] = size(C);
    d = bsxfun(@rdivide,C,C(:,1));
    X = ( [zeros(n,1) ones(n,1) -d(:,3) -d(:,2)] ).'; %'
    Y = reshape(X,2,[]);
    rt = blockproc(Y,[2 2],@(x) eig(x.data) ).'; %'
    condition = rt == abs(rt);
    out = rt(condition);
end
function out = sepideh(C)  %sepideh
    [n,~] = size(C);
    out = zeros(n,1);
    for ii = 1:n
        c = C(ii,:);
        d = c./c(1);
        rt = eig([0 -d(3);1 -d(2)]);
        condition = rt == abs(rt);
        rt = rt(condition);
        if isempty(rt)
            rt = 0;
        end
        rt = max(rt);
        out(ii) = rt;
    end
end

---

ans =

    12.2086  %// thewaywewalk
     9.2176  %// sepideh

---

我仍然没有看到 if 条件的意义。