MATLAB: Map a function over specific dimensions of matrix

functionMATLABmatrixoptimization

Consider an N dimensional matrix M. Consider a function f that eats k-dimensional matrices.
I want to build a function F such that the following code :
B=F(M,@f,k);
Returns a matrix B having N-k dimensions. The coefficient B(i_1,…,i_{N-k}) is equal to f(squeeze(M(i_1,…,i_{N-k},:))
To do so I created the following function :
function B = mapping(M, f, k)
% MAPPING map a function f that must have a k dimension matrix as input to 
% the last k dimensions of a matrix M
% B = mapping(M,@f,k) returns a matrix M of N-k dimensions where N is the
% dimension of M such that B(i1,..,i_{N-k})=f(M(i1...i_{N-k},:))
size_M=size(M);
% We start by creating an intermediate matrix A that will be required in
% the function arrayfun. This matrix has the same size as the returned B
% matrix but verifies A(k)=k where k is the 1 dimensional indice
if(numel(size_M)-k>1) 
    % general case where A (and B) will have at least 2 dimensions
    A=[1:prod(size_M(1:end-k))];
    A=reshape(A, size_M(1:end-k));
elseif(numel(size_M)-k==1)
    A=[1:size_M(1)];
else 
    error("The function to map maps all the dimensions of the matrix or more.");
end
      function fct = F(Ai)
          % We first take the subscript associated to the number Ai
          b_dim=numel(size_M(1:end-k)); 
          subscripts=cell(b_dim,1);
          % We take the indices corresponding to Ai
          [subscripts{:}]=ind2sub(size_M(1:end-k),Ai);
          subscripts=subscripts.';
          % We then apply the function f to the good element of matrix M
          fct=f(squeeze(M(subscripts{:},:)));
      end
B=arrayfun(@F, A);
end
My question are :
Above all : does it actually already exist a matlab function doing the purpose I want ? Or am I forced to create it ?
It it doesn't already exist such a function I would like to know if there are tricks to speed up what I want, because here it is too slow for my purpose.
If you want to make a try :
M=rand(8,10,2)
% The matrix B will have 2 dimensions (B is a 8x10 matrix) and B(i,j)=norm(squeeze(M(i,j,:)))
B=mapping(M,@norm,1);

Best Answer

You would have to modify your function F like this to achieve what I think you want:
      function fct = F(Ai)
          % We first take the subscript associated to the number Ai
          b_dim=numel(size_M(1:end-k)); 
          subscripts=cell(1, b_dim);
          % We take the indices corresponding to Ai
          [subscripts{:}]=ind2sub(size_M(1:end-k),Ai);
          % replicate individual indices to the size of the k-dimensional submatrices
          size_k = size_M(end-k+1:end);
          subscripts = cellfun(@(idx) repmat(idx, size_k), subscripts, 'UniformOutput', false);
          % emulate colon for the remaining dimension by explicitly listing all indices
          k_subscripts = arrayfun(@(sz) 1:sz, size_k, 'UniformOutput', false);
          [k_subscripts{:}] = ndgrid(k_subscripts{:});
          % convert subscripts to linear indices and index
          M_k = M(sub2ind(size(M), subscripts{:}, k_subscripts{:}));
          % We then apply the function f to the good element of matrix M
          fct=f(M_k);
      end
That is you have to translate your subscripts back to linear indices and explicitly list all the indices for the remaining dimensions instead of an unknown number of colons.
However, yes the whole thing is way overly complicated. You only need to do:
function B = mapping(M, f, k)
   B = cellfun(@(M_k) f(squeeze(M_k)), num2cell(M, k+1:ndims(M)));
end
Alternatively, you could permute the dimensions which would avoid the squeeze in the cellfun (and hence the cost of the anonymous function) but require a squeeze at the end. As permute is expensive, you may not gain anything:
function B = mapping(M, f, k)
   B = squeeze(cellfun(f, num2cell(permute(M, [k+1:ndims(M), 1:k]), 1:k)));
end
Related Question