"Reshaping" matrix

9 views (last 30 days)
Fredrik P on 1 Nov 2015
Edited: Philip Borghesani on 17 Nov 2015
I have a matrix P_ind that is AxB (where each element is an integer in the range 1:C) and a matrix xgrid that is CxDxB. Is there a smart way to create a new matrix xgrid2 that is AxDxB without resorting to the for loop solution below?
xgrid2=NaN(A,D,B);
for i=1:A
for j=1:B
xgrid2(i,:,j)=xgrid(P_ind(i,j),:,j);
end
end
Edit: Since there were comments on the clarity of my question, let me try again. What I want to do is create the array xgrid2 in the code below, but in a more succint and most of all faster way than the nested for loops in my approach.
A=10000; % no. of individuals
B=65; % maximum age
C=30; % gridpoints for P
D=50; % gridpoints for S
P_ind=randi(C,A,B);
xgrid=randi(5000,A,D,B);
xgrid2=NaN(A,D,B);
for i=1:A
for j=1:B
xgrid2(i,:,j)=xgrid(P_ind(i,j),:,j);
end
end
Fredrik P on 17 Nov 2015
Please excuse my late reply. My initial code was just a sketch for you to get a feeling of what I wanted to do. That was stupid. Now I have added complete, working code.

Philip Borghesani on 17 Nov 2015
Edited: Philip Borghesani on 17 Nov 2015
Simpler and faster solution get rid of the outer loop:
xgrid3=NaN(A,D,B);
tic
for j=1:B
xgrid3(:,:,j)=xgrid(P_ind(:,j),:,j);
end
t2=toc
For me the time went from 3.1 second to 0.13 seconds.
Some performance concepts:
• Not all for loops are bad
• Try to not loop over the first indices of an array.
• If you must loop over the first indices then make that loop the inner loop to traverse memory in the natural order.
• Any solution that uses sub2ind is probably not optimal because of the extra work it must do and because it is rather slow.
Using sub2ind is a programming design trade-off because it may make the code easier to understand and debug.

Kelly Kearney on 17 Nov 2015
My rule of thumb for this type of problem is to permute and reshape the matrices in such a way that you can access the dimension(s) of interest via linear indexing rather than subscripts.
The data:
A=10000; % no. of individuals
B=65; % maximum age
C=30; % gridpoints for P
D=50; % gridpoints for S
P_ind=randi(C,A,B);
xgrid=randi(5000,A,D,B);
tic;
xgrid2=NaN(A,D,B);
for i=1:A
for j=1:B
xgrid2(i,:,j)=xgrid(P_ind(i,j),:,j);
end
end
t(1) = toc;
New way (permute, reshape, index, unreshape, unpermute)
tic;
xgrdtmp = reshape(permute(xgrid, [1 3 2]), [], D);
[ii,jj] = ndgrid(1:A,1:B);
pidx = sub2ind(size(P_ind), ii, jj);
idx = sub2ind([A,B], P_ind(pidx), jj);
xgrid3 = permute(reshape(xgrdtmp(idx(:),:), [A B D]), [1 3 2]);
t(2) = toc;
Check
check = isequal(xgrid2, xgrid3);
fprintf('Time (old): %f\nTime (new): %f\nEqual: %d\n', t,check);
...
Time (old): 2.325150
Time (new): 0.441705
Equal: 1

Nitin Khola on 3 Nov 2015
My understanding is that you want to reshape an array that is of size "CxDxB", xgrid to another array of size "AxDxB". So that you can use "reshape" you need the total number of elements to remain the same i.e. A=C in your case. This implies you will have values repeating (I am assuming A>C). You need to use "repmat" to do it in that particular dimension. Refer to the following documentation for details: http://www.mathworks.com/help/matlab/ref/repmat.html http://www.mathworks.com/help/matlab/ref/reshape.html
Fredrik P on 17 Nov 2015
Please excuse my late reply. I think that you understand what I want to do. I can't figure out how to do it in any other way than I already do. I have added working code (as opposed to my initial "sketch") in my edit to the question, if you would care to give it a try. Many thanks!