Compute 1D IDFTs of a subset of 2D slices in a 5D Matrix

I want to compute 1D IDFTs of a subset of 2D slices within a 5D matrix yc. The dimensions of the 5D matrix yc are N1xN2xN3xN4xL. [a,b,c,d] is a permutation of [1,2,3,4] and [Na,Nb,Nc,Nd] is the corresponding permutation of [N1,N3,N3,N4]. V_arr is a sparse N1xN2xN3xN4 logical matrix that selects a sparse subset of elements in each N1xN2xN3xN4 slice of yc. The 1D IDFTs are computed along dimension d, so that each 1D IDFT is of length Nd. nzd = nnz(Ed), where Ed = (sum(V_arr,d) > 0) determines the subset of 2D slices within yc whose 1D IDFTs along dimension d need to be computed.

The code below achieves this in two different ways. The first method only computes nzd*L 1D IDFTs along dimension d. The second method comutes Na*Nb*Nc*L 1D IDFTs along dimension d. In the first method, is there a more efficient way to perform the "extract" and "unpack" steps, which are achieved with for loops in the implementation below?

N1 = 64; N2 = 32; N3 = 64; N4 = 16; L = 10;
%N1 = 128; N2 = 128; N3 = 256; N4 = 64; L = 4;
%N1 = 32; N2 = 4; N3 = 36; N4 = 12; L = 4;
%N1 = 4; N2 = 3; N3 = 5; N4 = 2; L = 4;
N = N1*N2*N3*N4;
N_vec = [N1,N2,N3,N4];
% Construct V_arr.
C = rand(N1,N2,N3,N4);
V_arr = (C > .999);
nnz_V_arr = nnz(V_arr);
fprintf('Number of nonzeros in V_arr: %d. Number of elements in V_arr: %d. Sparsity of A: %1.2f%%.\n',nnz_V_arr,N,(N-nnz_V_arr)/N*100);
dim_perm = randperm(4); % Random permutation of [1,2,3,4].
fprintf('The permutation of the 4 dimensions is: %d %d %d %d.\n',dim_perm(1),dim_perm(2),dim_perm(3),dim_perm(4));
N_perm = N_vec(dim_perm);
a = dim_perm(1); b = dim_perm(2); c = dim_perm(3); d = dim_perm(4);
Na = N_perm(1); Nb = N_perm(2); Nc = N_perm(3); Nd = N_perm(4);
Ed = (sum(V_arr,d) > 0); % Nonempty dimension d vectors in the virtual array.
nzd = nnz(Ed);
fsEd = find(squeeze(Ed));
switch d
    case 1
        [i1,i2,i3] = ind2sub([N2 N3 N4],fsEd);
    case 2
        [i1,i2,i3] = ind2sub([N1 N3 N4],fsEd);
    case 3
        [i1,i2,i3] = ind2sub([N1 N2 N4],fsEd);
    case 4
        [i1,i2,i3] = ind2sub([N1 N2 N3],fsEd);
end
Nabc = Na*Nb*Nc;
fprintf('The number of method 1 vs 2 IDFTs of size %d is: %d vs %d. The ratio 2/1 is: %1.2f. The percentage 1/2*100 is: %1.3f%%.\n',Nd,nzd*L,Nabc*L,Nabc/nzd,nzd/Nabc*100);
rV_arr = repmat(V_arr,1,1,1,1,L); % rV_arr is a binary-valued N1xN2xN3xN4xL matrix.
yc = rand(N1,N2,N3,N4,L)+1j*rand(N1,N2,N3,N4,L); % yc is a complex-valued N1xN2xN3xN4xL matrix.
y = zeros(N1,N2,N3,N4,L);
w = zeros(nzd,Nd,L);
J = 10; % Repeat J iterations of each method, soley for timing the two different methods.
%% Fast method.
ST1 = tic;
for j = 1:J
    % Extract each 2D slice from yc.
    for k = 1:nzd
        switch d
            case 1
                w(k,:,:) = yc(:,i1(k),i2(k),i3(k),:);
            case 2
                w(k,:,:) = yc(i1(k),:,i2(k),i3(k),:);
            case 3
                w(k,:,:) = yc(i1(k),i2(k),:,i3(k),:);
            case 4
                w(k,:,:) = yc(i1(k),i2(k),i3(k),:,:);
        end
    end
    
    z = ifft(w,Nd,2); % Dimension d IDFTs of size Ndx1. Only perform nzd*L dimension d IDFTs, each of size Ndx1.
    
    % Repack each 2D slice in z into y.
    for k = 1:nzd
        switch d
            case 1
                y(:,i1(k),i2(k),i3(k),:) = z(k,:,:);
            case 2
                y(i1(k),:,i2(k),i3(k),:) = z(k,:,:);
            case 3
                y(i1(k),i2(k),:,i3(k),:) = z(k,:,:);
            case 4
                y(i1(k),i2(k),i3(k),:,:) = z(k,:,:);
        end
    end
    
    out = y(rV_arr); % Downselect.
end
ET1 = toc(ST1);
%% Slow method.
ST2 = tic;
for j = 1:J
    y2 = ifft(yc,Nd,d); % Dimension d IDFTs of size Ndx1. Perform Na*Nb*Nc*L dimension d IDFTs, each of size Ndx1.
    out2 = y2(rV_arr); % Downselect.
end
ET2 = toc(ST2);
err = norm(out-out2);
spdup = ET2/ET1;
fprintf('The error between methods 1 and 2 is: %1.3e. The speedup 2/1 is: %1.2f.\n',err,spdup);

Answer 1

[Shishir Reddy]

Hi Stuart

As per my understanding, you would like to compute 1D inverse DFTs (IDFTs) only along a specific dimension 'd' for the selected slices, i.e., those that contain non-zero entries in V_arr.

The efficiency of the "extract" and "repack" steps can be improved by avoiding explicit 'for' loops using 'permute' and 'reshape'. The key idea is to bring the target dimension d (along which you compute the IDFT) to a fixed position (e.g. dimension 2), and then extract slices using vectorized indexing.

1. Permute yc so that dimension d becomes dimension 2

% Suppose we want to bring dim `d` to position 2
perm = [a, b, c, d, 5]; % Original dim order
[~, inv_perm] = sort(perm);

yc_perm = permute(yc, perm); % yc_perm has shape [Na, Nb, Nc, Nd, L]

2. Reshape into a 3D matrix for vectorized IDFT

yc_reshaped = reshape(yc_perm, [], Nd, L); % size: [Nabc, Nd, L]

3. Select only the required slices (nzd)

fsEd = find(squeeze(Ed)); % indices of non-zero slices
selected_yc = yc_reshaped(fsEd,:,:); % size: [nzd, Nd, L]

4. Compute IDFT along dimension 2 (now fixed)

z = ifft(selected_yc, Nd, 2); % size: [nzd, Nd, L]

For more information regarding 'permute' function in MATLAB, kindly refer the following documentation -

https://www.mathworks.com/help/control/ref/inputoutputmodel.permute.html

I hope this helps.

[Stuart Rogers]

Your code doesn't show how to repack z into y. Your code doesn't use inv_perm.

[Stuart Rogers]

I was able to get your code to work with my repacking for loop, by making the following modification to the variable perm in your implementation: perm = [sort([a, b, c]), d, 5];