Looping through a matrix of unknown dimensions

Hi everyone, I know similar questions have been asked before but I can't find anything that would be relevant for my particular situation. I need to write a function that would retrieve a vector from a matrix of unknown size, operate on it, and then return the results into an output matrix of size equal to the starting matrix. Here is an example of what I need:

start_matrix = rand(5,6,7,8,1000);
s = size(start_matrix);
output_matrix = zeros(s);

for a = 1:s(1)
    for b = 1:s(2)
        for c = 1:s(3)
            for d = 1:s(4)
                operate_data = start_matrix(a,b,c,d,:); % retrieving the vector of length=1000
                
                outputmatrix(a,b,c,d,:) = fft(operate_data); % using fast fourier transform as an example, there are a few more calculations here
            end
        end
    end
end

Now the number of dimensions can vary and the size of them can vary too, meaning that the number of loops needed will be variable. So my question is, how do I make this work with e.g. size = [10,15,20,500], where I loop through the first three dimensions, doing calculations on vectors of length = 500?

I've been trying to implement something like in this example by Jan:

https://uk.mathworks.com/matlabcentral/answers/345551-function-with-varying-number-of-for-loops

But I don't know how to make it work when I am not operating on each element, but rather a vector, each time.

Help is greatly appreciated! Thanks.

[KSSV]

Did you check this ?

https://in.mathworks.com/matlabcentral/answers/386638-how-to-loop-through-an-unknown-number-of-matrix-dimensions

[Greta]

Yes, I saw this but I'm completely lost with Jan's answer to that one... Everything from what the variables refer to (e.g. XX?) to how permutation is relevant here, or why the dimension to operate on needs to be moved to the start. Can anyone help me decipher this?

Answer 1

[Jan]

What about reshaping the array?

start_matrix = rand(5,6,7,8,1000);
s            = size(start_matrix);
M            = reshape(start_matrix, [], s(end));
Mfft         = fft(M, [], 2);
Result       = reshape(Mfft, s);

isequal(Result, output_matrix)
>> logical(1)

This is much faster than a loop. But if you really need a set of loop with distinct dimensions, this is the way to go:

X   = rand(5,6,7,8,1000);
Out = zeros(size(X));  % Pre-allocate the output

% Prepare index vector as cell, e.g. {1,2,3,4,':'} :
nv   = ndims(X) - 1;  % Exclude last dimension
v    = [repmat({1}, 1, nv), {':'}];
vLim = size(X);       % Loop until size (last element not used)

ready = false;
while ~ready
  % The actual operation:
  Out(v{:}) = fft(X(v{:}));
  
  % Update the index vector:
  ready = true;        % Assume that the WHILE loop is ready
  for k = 1:nv         % Loop through dimensions
    v{k} = v{k} + 1;   % Increase current index by 1
    if v{k} <= vLim(k) % Does it exceed the length of current dim?
       ready = false;  % No, WHILE loop is not ready now
       break;          % v(k) increased successfully, leave "for k" loop
    end
    v{k} = 1;          % Reset v{k}, proceed to next k
  end
end

The idea is to create a cell vector, which contains the indices, e.g. {1,2,3,4,':'}. Then this is equivalent:

index = {1, 2, 3, 4, ':'}
A(1,2,3,4,:)
A(index{:})

Now tis index vector is used and updated in the for loop: The first element is increased by 1. If it does not exceed size(X, 1), the next iteration is performed. If it exceeds size(X, 1), the element is reset to 1 and the next one is increased.

[Guillaume]

Oh, yes much simpler!

I'm leaving my answer as it shows some interesting manipulation technique, but use Jan's method.

[Greta]

Brilliant, that's exactly what I wanted. I will still need to include a loop in there (I think) to perform other calculations, but I can do that looping over M. Thank you so much!

[Greta]

Thanks, I will implement the second thing also.

Answer 2

[Guillaume]

I can see two ways of solving your problem.

The first way, would be to permute the dimensions, so that your : is move to the first dimension, This way, you can linear indices. For m a matrix of size MxNxPx..., m(:, 1, 1, 1, ...) is the same as m(1:M), m(:, 2, 1, 1, ...) is m((1:M) + M), etc. so:

m = rand(5,6,7,8,1000);
%move the colon dimension to the first one:
m = permute(m, [ndims(m), 1:ndims(m) - 1]);
output = zeros(size(m));
%use linear indexing
sz = size(m);
indices = 1:sz(1); 
for stride = 1:prod(sz(2:end))
    operate_indices = indices + (stride-1)*sz(1);
    output(operate_indices) = fft(m(operate_indices));
end
%at the end you can restore the original ordering if you want:
output = permute(output, [2:ndims(m), 1]);

The other way is to use the expansion of cell arrays into comma-separated list and store the indices into a cell array:

m = rand(5,6,7,8,1000);
%calculate all indices excluding last dimension
sz = size(m)
indices = arrayfun(@(s) 1:s, sz(1:end-1), 'UniformOutput', false);
[indices{:}] = ndgrid(indices{:});  %cartesian product of all indices of all dimensions. Using expansion of cell array into csl
%convert into a 2D matrix
indices = reshape(cat(numel(indices)+1, indices{:}), numel(indices{1}), []; 
% at this point each row of indices, is a list of index of each dimension. Split back into cell array
indices = num2cell(indices);
%iterate over the indices combination
output = zeros(size(m));
for row = 1:size(indices, 1)
    output(indices{row, :}, :) = fft(m(indices{row, :}, :));
end

[Greta]

Thanks for this, I will be using Jan's method (as you suggested as well) but hopefully one day this too will make sense to me :)

Answer 3

[Andrei Bobrov]

For your case:

start_matrix = rand(5,6,7,8,1000);
nd = 1:ndims(start_matrix);
a = permute(start_matrix,nd([end,1:end-1]));
out_matrix = permute(fft(a),nd([2:end,1]));