CODEWITHSUNDEEP
matlabvectorization
Ka Wa Yip
Ka Wa Yip

Reputation: 3009

Vectorization: matrix array multiplication element wise one by one

I have a matrix:

R = [0 -1;1 0]; 

array = 1:1:10;

Also x0 = [2;1]

How can I obtain another array in the most efficient way without loop?

array2 = [expm(1*R) expm(2*R) expm(3*R) .... expm(10*R)];

Then I want to obtain array3 of dimension 2 by 10 such that:

array3 = [expm(1*R)*x0 expm(2*R)*x0 expm(3*R)*x0 .... expm(10*R)*x0];

Upvotes: 3

Views: 168

Answers (3)

ammportal
ammportal

Reputation: 991

Well I see that the matrix R that you have is 2x2. In case it is always 2x2, then you can use the following function (Wikipedia) to calculate the exponential:

function output = expm2d(A)
% Assuming t = 1 from Evaluation by Laurent series (https://en.wikipedia.org/wiki/Matrix_exponential#Evaluation_by_Laurent_series)
s = trace(A) / 2;
q = sqrt(-det(A - s*eye(size(A))));
output = exp(s) * ((cosh(q) - s * sinh(q) / q) * eye(size(A)) + (sinh(q) * A / q));
end

Using the excellent comparison function provided by thewaywewalk, I got the following results:

When using expm:

>> bench
ans =
    0.0181 %// rahnema
    0.1075 %// thewaywewalk arrayfun
    0.1139 %// thewaywewalk accumarray

When using expm2d:

>> bench
ans =
    0.0048 %// rahnema
    0.0161 %// thewaywewalk arrayfun
    0.0222 %// thewaywewalk accumarray

As you can see, using the function for 2d matrices leads to a 10x decrease in the runtime. Of course, this cannot be used when R is not 2x2.

Edit: When using expm2d for A = 1:100:

>> bench
ans =
    0.1379 %// rahnema
    0.1415 %// thewaywewalk arrayfun
    0.1756 %// thewaywewalk accumarray

Upvotes: 2

Robert Seifert
Robert Seifert

Reputation: 25232

I still don't know if I got your question right. Here are two solutions which are not fully vectorized, but fairly fast:

R = [0 -1;1 0];
A = 1:1:10;
x0 = [2;1];

%// option 1
temp = arrayfun(@(x) (expm(R*x)*x0).', A, 'uni', 0);
array3 = vertcat( temp{:} )

%// option 2
temp = accumarray( (1:numel(A)).', A(:), [], @(x) {(expm(R*x)*x0).'})
array3 = vertcat( temp{:} )

Benchmark

I haven't considered Leander's Answer as it doesn't calculate array3:

function [t] = bench()
    R = [0 -1;1 0];
    A = 1:1:10;
    x0 = [2;1];

    % functions to compare
    fcns = {
        @() compare1(A,R,x0);
        @() compare2(A,R,x0);
        @() compare3(A,R,x0);
    };

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

function array3 = compare1(A,R,x0)  %rahnema1
    n = numel(A);
    array3 = reshape(expm(kron(diag(A),R))*repmat(x0,n,1),2,[])
end
function array3 = compare2(A,R,x0)  %thewaywewalk 1
    temp = arrayfun(@(x) (expm(R*x)*x0).', A, 'uni', 0);
    array3 = vertcat( temp{:} )
end
function array3 = compare3(A,R,x0)  %thewaywewalk 2
    temp = accumarray( (1:numel(A)).', A(:), [], @(x) {(expm(R*x)*x0).'});
    array3 = vertcat( temp{:} )
end

Results:

for A = 1:1:10;

0.1006   %// rahnema
0.2831   %// thewaywewalk arrayfun
0.3103   %// thewaywewalk accumarray

As kron gets really slow for large arrays, the benchmark results change for A = 1:1:100;:

4.0068   %// rahnema
1.8045   %// thewaywewalk arrayfun
2.4257   %// thewaywewalk accumarray

Upvotes: 1

rahnema1
rahnema1

Reputation: 15867

From wikipedia:

If a matrix is diagonal its exponential can be obtained by exponentiating each entry on the main diagonal.

Given that a block diagonal matrix can be created from {1*R, 2*R,...} then its exponential can be obtained and reshaped to a [2 * n] and it can be multiplied by x0. However its performance may be worse than for loop.

R = [0 -1;1 0];
array = 1:1:10;
x0 = [2;1]
n = numel(array);
result = reshape(expm(kron(spdiags(array.',0,n,n),R))*repmat(x0,n,1),2,[]);

For array of small size (less than 70 elements) full matrix is more efficient:

result = reshape(expm(kron(diag(array),R))*repmat(x0,n,1),2,[]);

Upvotes: 2

Related Questions