What is the fastest way to apply a function on each row of a matrix?

Question

I have searched so far and I know that there are several ways (1, 2, 3, 4) so far I have used the following code:

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

Where MaxPositiveRoot2DegreePolynomial is the following function:

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

Where QuadraticCoefficients is a 62271x3 matrix each row containing the coefficients a, b, c of a general quadratic equation. ax^2+bx+c Regarding the problem that I'm solving, all the roots will be real and so I've used a modified function for finding roots in order to not waste time to check if a root is real or not.
By profiling the code, I've found that each run of the function MaxPositiveRoot2DegreePolynomial takes about 0.000047 seconds which will be about 2.914925371 seconds for 62271 runs. But the code

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

takes 3.198597803 to run. So about 0.283672431 is taken just by the arrayfun. I want to see if there is any way to lessen this time?

Have you tried a loop? While `arrayfun` (and its related `*fun` family) is nice for writing compact code, they add error checking overhead that is often significant. Runtime of a basic `for` loop through `QuadraticCoefficients` is about half the `arrayfun` implementation in R2015b. — sco1, May 10 '16 at 16:31
If you spend more than 10 minute trying to optimize something that takes 3 seconds, there's a good chance you'll spend more time optimizing than you will ever spend running it, unless perhaps you're writing something that has to run hyper-optimized in real time, in which case, you probably shouldn't be using Matlab anyway. — Matthew Gunn, May 10 '16 at 16:45
Can't you just use a vectorized version of the quadratic formula? `sqrt`'s gotta have less overhead than `eig`. — TroyHaskin, May 10 '16 at 16:50
@MatthewGunn very good idea then I am going to write a question to help me decide whether to use matlab or not. — Sepideh Abadpour, May 10 '16 at 17:00
@TroyHaskin in fact I want to MEX the [gsl_poly_solve_quadratic](http://www.gnu.org/software/gsl/manual/html_node/Quadratic-Equations.html#Quadratic-Equations) and use that instead of `MaxPositiveRoot2DegreePolynomial`. I'm asking this question to just find the fastest way to apply a general function on each row of a matrix — Sepideh Abadpour, May 10 '16 at 18:48

Robert Seifert · Accepted Answer · 2016-05-11T07:27:58.527

This is a perfect example, which shows where you should just keep using your loop. My intended vectorized solution is slower than your loop approach finally.

But it could depend you on your Matlab version, as I'm already using Matlab R2015b with the optimzed JIT-Compiler, in older versions the vectoried approach is probably faster.

[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);

Benchmark

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

I still don't see the point of the if-condition.

I've uploaded that array for you [link](https://drive.google.com/file/d/0Bx_9vCAvXd5kcXZOeHl2ZFpFRE0/view?usp=sharing) — Sepideh Abadpour, May 10 '16 at 18:38
In fact I'm doing some more experiments for my special problem. And I haven't yet used your code yet. But in general if your answer is faster than mine. I accept it for the time being — Sepideh Abadpour, May 29 '16 at 22:40

What is the fastest way to apply a function on each row of a matrix?

1 Answers1

Benchmark