Vectorization: element of multiplying an array matrix one by one
I have a matrix:
R = [0 -1;1 0];
array = 1:1:10;
Also x0 = [2;1]
How can I get another array in the most efficient way without looping?
array2 = [expm(1*R) expm(2*R) expm(3*R) .... expm(10*R)];
Then I want to get the
array3
dimensions 2 by 10, so:
array3 = [expm(1*R)*x0 expm(2*R)*x0 expm(3*R)*x0 .... expm(10*R)*x0];
source to share
I still don't know if I got your question right. Here are two solutions that 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 compute 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
it gets very slow for large arrays, the test results vary for A = 1:1:100;
:
4.0068 %// rahnema
1.8045 %// thewaywewalk arrayfun
2.4257 %// thewaywewalk accumarray
source to share
From wikipedia :
If the diagonal of the matrix is a diagonal, its exponent can be obtained by highlighting each entry on the main diagonal.
Given that a block diagonal matrix can be created from {1*R, 2*R,...}
, then its exponent can be obtained and changed by a [2 * n]
, and it can be multiplied by x0
. However, its performance may be worse than for a 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
small size (less than 70 elements) the full matrix is more efficient:
result = reshape(expm(kron(diag(array),R))*repmat(x0,n,1),2,[]);
source to share
Ok I can see that the matrix R
has 2x2. If it is always 2x2, you can use the following function ( Wikipedia ) to calculate the exponent:
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 compare 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 results in a 10x faster execution time. Of course, this cannot be used when R is not 2x2.
Edit: When used expm2d
for A = 1:100
:
>> bench
ans =
0.1379 %// rahnema
0.1415 %// thewaywewalk arrayfun
0.1756 %// thewaywewalk accumarray
source to share