How can I efficiently extract a max along two dimensions, from a 4 dimensional array, without writing a double loop?

I have the following loop. It picks values from a matrix A(nb,nd). For each vaue of this matrix A it constructs a matrix B = F(A(i,j),Z), where Z is some matrix of a given dimension Z(n,m), and F is some operator that imposes constraints which vary with the index pair (i,j). Then I find the maximum value of these matrices B, and preserve 3 things: The maximum value of B, BMAX = bmax(i,j), the column index where each max is found, C = C(i,j), and the row index, R=R(i,j). QUESTION: The question is the following: what is the most efficient way of eliminating the double loop? I assume there are 2 steps. the first to build an array BARRAY(nb,nd,n,m) directly without a loop, (it is of course obvious how to build the array with a loop) and the second to maximize BARRAY along the appropriate dimensions and collect the values and indices. Here is the code below: PS: eliminating the double max WITHIN the loop by using reshape and colon, and then using ind2sub actually makes the double loop slower, for the dimensions I am working with.

         BMAX = ones(nb,nd);
         C = ones(nb,nd);
         R = ones(nb,nd);
         for i = 1:nb
         for j = 1:nd;     
               % We write B = F(A(i,j), Z), so that
              [v1 ix1] = max(B);
              [BMAX(i,j) ix2] = max(v1); % ix2 is column index.
              C(i,j) = ix2;
              R(i,j) = ix1(ix2);                    
         end
         end

Answer 1

[Matt J]

In each iteration of the loop, you generate a matrix,

   B_ij = F(A(i,j),Z)

but its 2D nature is irrelevant, as is the 2D nature of A(i,j).

         dims=[nb,nd];
         BMAX = nan(dims);
         idx=nan(dims);

         for k=1:numel(A)

              Bij=F(A(k),Z);

              [BMAX(k),idx(k)]= max(Bij(:));

         end

         [C,R]=ind2sub([n,m],idx);

[Matt J]

A further question is whether you really need the loop over A(k) at all. If each Bij(:) is the column of some larger matrix Btotal, which can be computed in some vectorized fashion, you can just do

   [BMAX,idx]=max(Btotal,[],1);

    BMAX=reshape(BMAX,dims);
    [C,R]=ind2sub([n,m],reshape(idx,dims));

[Joao Ejarque]

Your first answer is sure to improve on what I have. I will think about your second comment. I will post something here sometime this weekend. Regards. J.

[Joao Ejarque]

The first answer improves matters but only slightly. It turns out that stacking the initial matrix and therefore getting rid of one of the indices reduces computational time from about 4 seconds to 3.8 seconds. It is a more compact code so it is of course better but only just. I am working on the second suggestiom. I am happy to send you the m files later if you are curious about the problem. I will be in touch.

[Matt J]

I'm guessing that you're not seeing a significant speed-up because the principal computational effort is the evaluation of Bij=F(A(k),Z). There must be heavy operations inside F().