Fast dataset manipulation in MATLAB

Problem definition

I have a dataset composed by m 2-dimensional vectors, stored in a 2xm matrix D. Given a rotation angle theta (in rads) and given a 2-dimensional translation vector p, I have to compute a (inverse) roto-translation of each vectors in the dataset. I'm wondering what is the fastest solution to perform this task because I've noticed that in my code the bottleneck is, indeed, caused by the computation of a roto-translation (99.8% of the total computation time).

To fix the notation, I denote the generic roto-translation as the following affine transformation

R    = [cos(theta) -sin(theta); sin(theta) cos(theta)];
bary = R' * (y - p);

where y is the original vector, bary is the roto-translated vector, R is the rotation matrix and p is the translation vector.

Naive solution

The simplest solution is to apply to each single vector in the dataset D the previous transformation

for j = 1 : m
    barD(:, j) = R' * (D(:, j) - p);
end

where barD is the rototranslated dataset. Hence, here we perform m "small" affine transformations.

A (hopefully) better solution

I'm not an expert in programming, but I understand, if I'm not wrong, that we should avoid to use for cycles as much as possible because MATLAB is optimized to manipulate matrices and most of the times we can express for cycles as matrix manipulations. If I'm not wrong, this concept is called code vectorization. In this case we can replace the previous for cycle with a single "big" affine transformation. This objective can be achieved by stacking the dataset vectors in one single "big" 2*m x 1 vector vecD.

vecD = reshape(D, [2*m 1]);

Then, we can compute the m rototranslations via the following "diagonal" trasformation

invRR   = kron(eye(m), R');
pp      = kron(ones(m, 1), p);
vecbarD = invRR * (vecD - pp);

where vecbarD is the vectorial representation of the rototranslated dataset barD. To conclude the task, vecbarD must be converted back in the 2 x m matrix representation barD. So far I'm using the following code, but probably it would be better to use reshape once again.

xidx = 1 : 2 : 2 * m - 1; 
yidx = 2 : 2 : 2 * m;
barD = [vecbarD(xidx)'; vecbarD(yidx)'];

Questions

1. I would like to know if what I've written makes sense and where I can find some resource about code vectorization
2. I would like to know if and how my solutions can be improved in terms of the computational time
3. Since each vector in the dataset is treated independently, I would like to know how to parallelize my code. A trivial solution is to replace the for of my naive solution with a parfor, but I don't know if this strategy is better than a vectorized one.

Answer 1

[Matt J]

The implementation that you wrote first,

R    = [cos(theta) -sin(theta); sin(theta) cos(theta)];
bary = R' * (y - p);

is already vectorized and internally parallelized and should have the best performance, at least on the CPU. Replication of R using kron() will just add more overhead and is unnecessary.

If you want further speed-up, you could consider converting R,y,p to gpuArrays assuming you have a high-performance graphics card.

[Matteo Tesori]

Thanks for your advice Matt. The two line of codes mentioned by you work only for a single vector y.

In order to deal with multiple vectors, i.e. the dataset D, I run

for j = 1 : m
    barD(:, j) = R' * (D(:, j) - p);
end

Is this the best way to manipulate D?

[Matt J]

The two line of codes mentioned by you work only for a single vector y.

I'm not sure why you think so. Have you tried it? It doesn't sound like you have.

D=rand(2,1e5);
p=rand(2,1); 
theta=pi/4;

tic
R    = [cos(theta) -sin(theta); sin(theta) cos(theta)];
barD = R' * (D - p);
toc

Elapsed time is 0.004413 seconds.

[Steven Lord]

Depending on how often you're doing this, manually transposing R as you define it and skipping transposing it when you perform the multiplication might save some time. This assumes sin(theta) is real and scalar.

D=rand(2,1e5);
p=[rand(2,1)]; 
theta=pi/4;

tic
R    = [cos(theta) -sin(theta); sin(theta) cos(theta)];
barD1 = R' * (D - p);
toc

Elapsed time is 0.025855 seconds.

tic
R    = [cos(theta) sin(theta); -sin(theta) cos(theta)];
barD2 = R * (D - p);
toc

Elapsed time is 0.004135 seconds.

norm(barD1-barD2)

ans = 0

[Matteo Tesori]

@Matt J: I didn't realize that

barD = R' * (D - p);

is equivalent to

for j = 1 : m
    barD(:, j) = R' (D(:, j) - p);
end

thanks for tip.

@Steven Lord: thanks, your suggestion further reduces the computational time