"Reshaping" matrix

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

[Triveni]

reshape(xgrid2',i,j,[])

tell your question clearly.

your script is not running.

[Fredrik P]

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.

Answer 1

[Philip Borghesani]

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.

Answer 2

[Kelly Kearney]

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

Your way (loops)

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

Answer 3

[Nitin Khola]

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]

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!